Fork me on GitHub

Changeset 39 in svn


Ignore:
Timestamp:
Nov 18, 2008, 10:30:58 AM (16 years ago)
Author:
severine ovyn
Message:

ok for ETmis

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Resolutions.cpp

    r26 r39  
    5252//------------------------------------------------------------------------------
    5353
    54 
    55 
    56 void PairingJet(TLorentzVector &JETSm, TLorentzVector JET,  vector<fastjet::PseudoJet> sorted_jetsS)
     54// //********************************** PYTHIA INFORMATION*********************************
     55
     56TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
     57 {
     58   TIter it((TCollection*)GEN);
     59   it.Reset();
     60   TRootGenParticle *gen1;
     61   TSimpleArray<TRootGenParticle> array,array2;
     62
     63   while((gen1 = (TRootGenParticle*) it.Next()))
     64      {
     65        array.Add(gen1);
     66      }
     67   it.Reset();
     68   bool tauhad;
     69   while((gen1 = (TRootGenParticle*) it.Next()))
     70    {
     71       tauhad=false;
     72       if(abs(gen1->PID)==15)
     73         {
     74            int d1=gen1->D1;
     75            int d2=gen1->D2;
     76
     77            if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
     78               {
     79                 tauhad=true;
     80                 for(int d=d1; d < d2+1; d++)
     81                    {
     82                      if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
     83                    }
     84               }
     85         }
     86        if(tauhad)array2.Add(gen1);
     87     }
     88   return array2;
     89 }
     90
     91
     92
     93void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET,  vector<fastjet::PseudoJet> sorted_jetsS)
    5794{
    5895  JETSm.SetPxPyPzE(0,0,0,0);
    … …  
    146183 
    147184  TRootGenParticle *particle;
    148   TRootETmis *etmisc;
    149  
    150   RESOLELEC *elementElec;
     185 
     186  RESOLELEC * elementElec;
    151187  RESOLMUON *elementMuon;
    152188  RESOLJET *elementJet;
    … …  
    172208  vector<fastjet::PseudoJet> inclusive_jetsS;
    173209  vector<fastjet::PseudoJet> sorted_jetsS;
     210
     211  vector<fastjet::PseudoJet> input_particlesT;//for FastJet algorithm
     212  vector<fastjet::PseudoJet> inclusive_jetsT;
     213  vector<fastjet::PseudoJet> sorted_jetsT;
     214
    174215  vector<PhysicsTower> towers;
    175216 
    176217  fastjet::JetDefinition jet_def;
    177218  fastjet::JetDefinition jet_defS;
     219  fastjet::JetDefinition jet_defT;
    178220  fastjet::JetDefinition::Plugin * plugins;
    179221  fastjet::JetDefinition::Plugin * pluginsS;
     222  fastjet::JetDefinition::Plugin * pluginsT;
    180223 
    181224  // set up a CDF midpoint jet definition
    … …  
    195238#endif
    196239 
     240 // set up a CDF midpoint jet definition
     241 #ifdef ENABLE_PLUGIN_CDFCONES
     242 pluginsT = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
     243 jet_defT = fastjet::JetDefinition(pluginsT);
     244 #else
     245 pluginsT = NULL;
     246 #endif
     247 
    197248 
    198249  // Loop over all events
    … …  
    212263      TrackCentral.clear();
    213264      TSimpleArray<TRootGenParticle> NFCentralQ;
     265   
     266      sorted_jetsS.clear();
     267      input_particlesS.clear();
     268      inclusive_jetsS.clear();
     269
     270      sorted_jetsT.clear();
     271      input_particlesT.clear();
     272      inclusive_jetsT.clear();
     273
     274      sorted_jets.clear();
    214275      input_particles.clear();
    215276      inclusive_jets.clear();
    216       sorted_jets.clear();
    217       input_particlesS.clear();
    218       inclusive_jetsS.clear();
    219       sorted_jetsS.clear();
    220277      towers.clear();
    221278     
    … …  
    228285          float eta = fabs(particle->Eta);
    229286         
    230           if(particle->Status == 1)
     287          if(particle->Status == 1 && eta < DET->MAX_CALO_FWD)
    231288            {
    232289              input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
    233290            }
    234291         
     292          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmis = PTmis + genMomentum;
    235293          // keeps only final particles, visible by the central detector, including the fiducial volume
    236294          // the ordering of conditions have been optimised for speed : put first the STATUS condition
    … …  
    241299              )
    242300            {
    243               PTmis = PTmis + genMomentum;//ptmis
     301              //if((pid != pNU1) && (pid != pNU2) && (pid != pNU3))PTmis = PTmis + genMomentum;//ptmis
     302
     303              Float_t nonS=0;
    244304              switch(pid) {
    245 
    246305              case pE: // all electrons with eta < DET->MAX_CALO_FWD
    247                 DET->SmearElectron(genMomentum);
     306                nonS = genMomentum.E();
     307                DET->SmearElectron(genMomentum);
     308                if(eta < DET->MAX_TRACKER){
     309                elementElec=(RESOLELEC*) branchelec->NewEntry();
     310                elementElec->NonSmeareE = nonS;
     311                elementElec->SmeareE = genMomentum.E();}
    248312                break; // case pE
    249313              case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
    … …  
    251315                break; // case pGAMMA
    252316              case pMU: // all muons with eta < DET->MAX_MU
     317                elementMuon = (RESOLMUON*) branchmuon->NewEntry();
     318                elementMuon->OneoverNonSmearePT = (1/genMomentum.Pt());
    253319                DET->SmearMu(genMomentum);
     320                elementMuon->OneoverSmearePT = (1/genMomentum.Pt());
    254321                break; // case pMU
    255322              case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
    … …  
    264331              if(pid != pMU)
    265332                {
    266 //                PTmisS = PTmisS + genMomentum;
    267333                  towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
    268                   if(genMomentum.Et()>0.5)input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     334                  input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
    269335                }
    270336             
    … …  
    286352            } // switch
    287353        } // while
    288      
    289      TLorentzVector TowerTLV(0.,0.,0.,0.);
     354   
     355     TLorentzVector Att(0.,0.,0.,0.);
     356     float ScalarEt=0;
    290357     for(unsigned int i=0; i < towers.size(); i++)
    291358       {
    292          TowerTLV.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
    293          if(TowerTLV.Et()>0.5)PTmisS = PTmisS + TowerTLV;
    294       }
    295 
    296  
     359         Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
     360         ScalarEt = ScalarEt + Att.Et();
     361         PTmisS = PTmisS + Att;
     362       }
     363//cout<<"non smeare "<<PTmis.Pt()<<" "<<PTmisS.Pt()<<endl;
     364//cout<<"smeare "<<PTmis.Px()<<" "<<PTmisS.Px()<<endl;
     365//cout<<"**********"<<endl;
     366
    297367      elementEtmis= (ETMIS*) branchetmis->NewEntry();
    298       elementEtmis->Et = (PTmis).Et();
    299       elementEtmis->EtSmeare = (PTmisS).Et()-(PTmis).Et();
     368      elementEtmis->Et = (PTmis).Pt();
     369      elementEtmis->Ex = (-PTmis).Px();
     370      elementEtmis->SEt = ScalarEt;
     371
     372      elementEtmis->EtSmeare = (PTmisS).Pt()-(PTmis).Pt();
     373      elementEtmis->ExSmeare = (-PTmisS).Px()-(PTmis).Px();
    300374     
    301375      //*****************************
    … …  
    316390        }
    317391     
     392      TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
     393
    318394      TLorentzVector JETSm(0,0,0,0);
    319395      for (unsigned int i = 0; i < sorted_jets.size(); i++) {
    … …  
    326402            elementJet->NonSmearePT = JET.Et();
    327403            elementJet->SmearePT = JETSm.Et()/JET.Et();
    328            
    329404          }
     405      }
     406      for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
     407        TLorentzVector JETT(0,0,0,0);
     408        JETT.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
     409        if(fabs(JETT.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS))
     410           {
     411            for(Int_t i=0; i<TausHadr.GetEntries();i++)
     412              {
     413                if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
     414                  {
     415                    elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
     416                    elementTaujet->EnergieCen = (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E());
     417                    elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JETT.Eta(),JETT.Phi());
     418                 }
     419              }
     420           }
     421         
    330422       
    331423      } // for itJet : loop on all jets

  • trunk/interface/FuncDef.h

    r27 r39  
    7979
    8080  string nom = histo.erase(0,histo.find(">>")+2);
    81   TH1F *h = new TH1F(nom.c_str(),"",200,0,3);
     81  TH1F *h = new TH1F(nom.c_str(),"",50,-3,3);
    8282
    8383  string all = min + " && " + max;
    8484  Analyze->Draw(temp.c_str(),all.c_str());
    8585  h->SetMarkerSize(0.6);
    86 //  double MeanFix = ;
    8786  double MeanFix = h->GetMean();
    88   double RangMin = h->GetMean()-h->GetRMS();
    89   double RangMax = h->GetMean()+h->GetRMS();
    90   TF1 *Gauss = new TF1("Gauss","gaus",RangMin,RangMax);
     87  TF1 *Gauss = new TF1("Gauss","gaus",-3,3);
    9188  Gauss->FixParameter(1,MeanFix);
    92   h->Fit("Gauss","R");
    93   h->Fit("Gauss","RI");
    94   h->Fit("Gauss","RI");
     89  h->Fit("Gauss","QR");
     90  h->Fit("Gauss","QRI");
     91  h->Fit("Gauss","QRI");
    9592  Double_t* params = Gauss->GetParameters();
    9693  rms=params[2];
    … …  
    10299}
    103100
    104 void GaussValuesETmis(TTree * Analyze,string histo,double &rms, double &mean, string min,string max)
     101void GaussValuesETmis(TTree * Analyze,string histo,double &rms, string min,string max)
    105102{
    106103  string temp = histo;
    … …  
    109106
    110107  string nom = histo.erase(0,histo.find(">>")+2);
    111   TH1F *h = new TH1F(nom.c_str(),"",50,-300,300);
     108  TH1F *h = new TH1F(nom.c_str(),"",20,-100,100);
    112109
    113110  string all = min + " && " + max;
    … …  
    117114  double RangMax = h->GetMean()+h->GetRMS();
    118115  TF1 *Gauss = new TF1("Gauss","gaus",RangMin,RangMax);
    119   h->Fit("Gauss","R");
    120   h->Fit("Gauss","RI");
    121   h->Fit("Gauss","RI");
     116  //TF1 *Gauss = new TF1("Gauss","gaus");
     117  h->Fit("Gauss","QR");
     118  h->Fit("Gauss","QRI");
     119  h->Fit("Gauss","QRI");
    122120  Double_t* params = Gauss->GetParameters();
    123121  rms=params[2];
    124   mean=params[1];
     122  //rms=h->GetRMS();
     123  //mean=params[1];
    125124  h->Draw("P");
    126125  h->GetXaxis()->SetTitle("E_{T}^{rec}-E_{T}^{MC} [GeV]");

  • trunk/routines/resolutions.C

    r27 r39  
    1111void JetResol()
    1212{
    13 
    14    setTDRStyle();
    15    gROOT->Reset();
    16 
    17    TFile *f1 = new TFile("JETResol.root","read");
    18    TTree *Analyze = (TTree*)f1->Get("Analysis");
    19  
    20    const Int_t numBin=7;
    21    double bins[numBin]={0,20,40,60,80,100,220};
    22    TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
    23 
    24    double rms[numBin];
    25    double mean[numBin];
    26 
    27    TCanvas *c1 = new TCanvas("c1","JET resol",0,0,1000,650);
    28    c1->cd(); int frame=0;
    29    c1->Divide(6,2);
    30 
    31    float x[numBin-1];
    32    float y[numBin-1];
    33    float ex[numBin-1];
    34    float ey[numBin-1];
    35 
    36    float finval=0;//valeur a remplir puis a fitter
    37 
    38    for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
    39       {
    40         TAxis *xaxis = ETSoverET->GetXaxis();
    41         float binCenter = xaxis->GetBinCenter(i+1);
    42         int binMin = (int)xaxis->GetBinLowEdge(i+1);
    43         int binMax = (int)xaxis->GetBinUpEdge(i+1);
    44         cout<<"binMin "<<binMin<<" binMax "<<binMax<<" bin center "<<binCenter<<endl;
    45         char tempMin[500];
    46         if(i==0)binMin=5;
    47         sprintf(tempMin,"JetPTResol.NonSmearePT > %d",binMin);
    48         string mystringMin(tempMin);
    49         char tempMax[500];
    50         sprintf(tempMax,"JetPTResol.NonSmearePT < %d",binMax);
    51         string mystringMax(tempMax);
    52         char tempName[500];
    53         sprintf(tempName,"(JetPTResol.SmearePT)>>hETSoverET%d",i);
    54         string mystringName(tempName);
    55 
    56         c1->cd(++frame);
    57         GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
    58         x[i]=binCenter;
    59         finval=rms[i]/mean[i];
    60         y[i]=(finval*100);
    61         ex[i]=0;
    62         ey[i]=0;
    63         cout<<"finval "<<finval<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl;
    64       }
    65 
    66    TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
    67    c2->cd();
    68 
    69    TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,170);
    70 
    71    TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
    72    gr11->Draw("AP");
    73    gr11->SetTitle("");
    74    gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
    75    gr11->GetYaxis()->SetRangeUser(0,50);
    76    gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})_{fit}/<E_{T}^{rec}/E_{T}^{MC}>_{fit}");
    77 
    78    Double_t* params = fitfun->GetParameters();
    79 
    80    gr11->Fit("user","R");
    81    gr11->Fit("user","RI");
    82    gr11->Fit("user","RI");
    83    gr11->Fit("user","RI");
    84    char tempResol[500];
    85    sprintf(tempResol,"#frac{#sigma(E_{T}^{rec}/E_{T}^{MC})}{<E_{T}^{rec}/E_{T}^{MC}>} = #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus #frac{%f}{E_{T}^{MC}} #oplus %f",params[1]/100,params[2]/100,params[0]/100);
    86 
    87    TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
    88    leg1->Draw();
    89 
    90    delete fitfun;
    91 }
    92 
     13  setTDRStyle();
     14  gROOT->Reset();
     15 
     16  TFile *f1 = new TFile("JET.root","read");
     17  TTree *Analyze = (TTree*)f1->Get("Analysis");
     18 
     19  const Int_t numBin=7;
     20  double bins[numBin]={0,20,40,60,80,100,220};
     21  TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
     22 
     23  double rms[numBin];
     24  double mean[numBin];
     25 
     26  TCanvas *c1 = new TCanvas("c1","JET resol",0,0,1000,650);
     27  c1->cd(); int frame=0;
     28  c1->Divide(6,2);
     29 
     30  float x[numBin-1];
     31  float y[numBin-1];
     32  float ex[numBin-1];
     33  float ey[numBin-1];
     34 
     35  float finval=0;//valeur a remplir puis a fitter
     36 
     37  for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
     38    {
     39      TAxis *xaxis = ETSoverET->GetXaxis();
     40      float binCenter = xaxis->GetBinCenter(i+1);
     41      int binMin = (int)xaxis->GetBinLowEdge(i+1);
     42      int binMax = (int)xaxis->GetBinUpEdge(i+1);
     43      char tempMin[500];
     44      if(i==0)binMin=5;
     45      sprintf(tempMin,"JetPTResol.NonSmearePT > %d",binMin);
     46      string mystringMin(tempMin);
     47      char tempMax[500];
     48      sprintf(tempMax,"JetPTResol.NonSmearePT < %d",binMax);
     49      string mystringMax(tempMax);
     50      char tempName[500];
     51      sprintf(tempName,"(JetPTResol.SmearePT)>>hETSoverET%d",i);
     52      string mystringName(tempName);
     53     
     54      c1->cd(++frame);
     55      GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
     56      x[i]=binCenter;
     57      finval=rms[i]/mean[i];
     58      y[i]=(finval*100);
     59      ex[i]=0;
     60      ey[i]=0;
     61    }
     62 
     63  TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
     64  c2->cd();
     65 
     66  TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,170);
     67 
     68  TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
     69  gr11->Draw("AP");
     70  gr11->SetTitle("");
     71  gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
     72  gr11->GetYaxis()->SetRangeUser(0,50);
     73  gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})_{fit}/<E_{T}^{rec}/E_{T}^{MC}>_{fit}");
     74 
     75  Double_t* params = fitfun->GetParameters();
     76 
     77  gr11->Fit("user","QR");
     78  gr11->Fit("user","QRI");
     79  gr11->Fit("user","QRI");
     80  gr11->Fit("user","QRI");
     81  char tempResol[500];
     82  if(params[2]==0.000000)sprintf(tempResol,"#frac{#sigma(E_{T}^{rec}/E_{T}^{MC})}{<E_{T}^{rec}/E_{T}^{MC}>} = #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus %f",params[1]/100,params[0]/100);
     83  else sprintf(tempResol,"#frac{#sigma(E_{T}^{rec}/E_{T}^{MC})}{<E_{T}^{rec}/E_{T}^{MC}>} = #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus #frac{%f}{E_{T}^{MC}} #oplus %f",params[1]/100,params[2]/100,params[0]/100);
     84 
     85  TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
     86  leg1->Draw();
     87
     88  TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
     89  Delphes->Draw();
     90 
     91 
     92  delete fitfun;
     93}
     94
     95void ElecResol()
     96{
     97 
     98  setTDRStyle();
     99  gROOT->Reset();
     100 
     101  TFile *f1 = new TFile("PTMIS.root","read");
     102  TTree *Analyze = (TTree*)f1->Get("Analysis");
     103 
     104  const Int_t numBin=9;
     105  double bins[numBin]={0,20,40,60,80,100,150,220,400};
     106  TProfile *ETSminusET = new TProfile("ETSminusET","Electron resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
     107 
     108  double rms[numBin];
     109  double mean[numBin];
     110 
     111  TCanvas *c3 = new TCanvas("c3","ELEC resol",0,0,1000,650);
     112  c3->cd(); int frame=0;
     113  c3->Divide(6,2);
     114 
     115  float x[numBin-1];
     116  float y[numBin-1];
     117  float ex[numBin-1];
     118  float ey[numBin-1];
     119 
     120  float finval=0;//valeur a remplir puis a fitter
     121 
     122  for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
     123    {
     124      TAxis *xaxis = ETSminusET->GetXaxis();
     125      float binCenter = xaxis->GetBinCenter(i+1);
     126      int binMin = (int)xaxis->GetBinLowEdge(i+1);
     127      int binMax = (int)xaxis->GetBinUpEdge(i+1);
     128      char tempMin[500];
     129      if(i==0)binMin=5;
     130      sprintf(tempMin,"ElecEResol.NonSmeareE > %d",binMin);
     131      string mystringMin(tempMin);
     132      char tempMax[500];
     133      sprintf(tempMax,"ElecEResol.NonSmeareE < %d",binMax);
     134      string mystringMax(tempMax);
     135      char tempName[500];
     136      sprintf(tempName,"(ElecEResol.NonSmeareE-ElecEResol.SmeareE)>>hETSoverET%d",i);
     137      string mystringName(tempName);
     138     
     139      c3->cd(++frame);
     140      GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
     141      x[i]=binCenter;
     142      finval=rms[i]/binCenter;
     143      y[i]=(finval*100);
     144      ex[i]=0;
     145      ey[i]=0;
     146    }
     147 
     148  TCanvas *c4 = new TCanvas("c4","ELEC resol",100,100,600,450);
     149  c4->cd();
     150 
     151  TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,400);
     152 
     153  TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
     154  gr11->Draw("AP");
     155  gr11->SetTitle("");
     156  gr11->GetXaxis()->SetTitle("E [GeV]");
     157  gr11->GetYaxis()->SetRangeUser(0,5);
     158  gr11->GetYaxis()->SetTitle("#sigma/E");
     159 
     160  Double_t* params = fitfun->GetParameters();
     161 
     162  gr11->Fit("user","QR");
     163  gr11->Fit("user","QRI");
     164  gr11->Fit("user","QRI");
     165  gr11->Fit("user","QRI");
     166  char tempResol[500];
     167  sprintf(tempResol,"#frac{#sigma}{E} = #frac{%f}{#sqrt{E}} #oplus #frac{%f}{E} #oplus %f",params[1]/100,params[2]/100,params[0]/100);
     168 
     169  TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
     170  leg1->Draw();
     171 
     172  TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
     173  Delphes->Draw();
     174 
     175  TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"t#bar{t} #rightarrow l^{+}l^{-}#nu#bar{#nu}b#bar{b}");
     176  events->Draw();
     177 
     178  delete fitfun;
     179}
     180
     181void TauJetInfo()
     182{
     183 
     184  setTDRStyle();
     185  gROOT->Reset();
     186 
     187  TFile *f1 = new TFile("PTMIS.root","read");
     188  TTree *Analyze = (TTree*)f1->Get("Analysis");
     189 
     190  TCanvas *ct = new TCanvas("ct","Tau information",0,0,1000,425);
     191  ct->cd(); int frame=0;
     192  ct->Divide(2,1);
     193 
     194  ct->cd(++frame);
     195  TH1F *tauEnergy =MakeNormTH1F(20,0.8,1,Analyze,"TauJetPTResol.EnergieCen>>tauEnergy",1, 0, 1,2,false);
     196  tauEnergy->Draw();
     197  tauEnergy->SetTitle("");
     198  tauEnergy->GetYaxis()->SetTitle("Fraction of events");
     199  tauEnergy->GetXaxis()->SetTitle("C_{#tau}");
     200 
     201 
     202  TPaveText *Delphes1 = MakeTPave(0.3,0.85,0.45,0.9,"MG/ME + Delphes");
     203  Delphes1->Draw();
     204  TPaveText *ctau = MakeTPave(0.3,0.75,0.45,0.8,"C_{#tau} = #frac{#sum E^{towers} (#DeltaR = 0.15)}{E^{jet}}");
     205  ctau->Draw();
     206  TPaveText *events1 = MakeTPave(0.3,0.65,0.45,0.7,"Events: t#bar{t} #rightarrow l^{+}#nul^{-}#bar{#nu}b#bar{b}");
     207  events1->Draw();
     208 
     209 
     210  ct->cd(++frame);
     211  TH1F *NumTrack =MakeNormTH1F(6,0,6,Analyze,"TauJetPTResol.NumTrack>>NumTrack",1, 0, 1,2,false);
     212  NumTrack->Draw();
     213  NumTrack->SetTitle("");
     214  NumTrack->GetYaxis()->SetTitle("Fraction of events");
     215  NumTrack->GetXaxis()->SetTitle("N^{tracks}");
     216 
     217  TPaveText *Delphes = MakeTPave(0.6,0.85,0.85,0.9,"MG/ME + Delphes");
     218  Delphes->Draw();
     219  TPaveText *numtracks = MakeTPave(0.6,0.75,0.85,0.8,"#DeltaR < 0.4, p_{T}^{track} > 2 GeV");
     220  numtracks->Draw();
     221  TPaveText *events = MakeTPave(0.6,0.65,0.85,0.7,"Events: t#bar{t} #rightarrow l^{+}#nul^{-}#bar{#nu}b#bar{b}");
     222  events->Draw();
     223   
     224}
    93225
    94226void ETmisResol()
    95227{
    96 
    97    setTDRStyle();
    98    gROOT->Reset();
    99 
    100    TFile *f1 = new TFile("PTMIS.root","read");
    101    TTree *Analyze = (TTree*)f1->Get("Analysis");
    102 
    103    const Int_t numBin=6;
    104    double bins[numBin]={10,20,30,50,80,100};
    105    TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-1000,1000);
    106 
    107    double rms[numBin-1];
    108    double mean[numBin-1];
    109 
    110    TCanvas *c5 = new TCanvas("c5","gerrors2",0,0,1000,650);
    111    c5->cd(); int frame=0;
    112    c5->Divide(6,2);
    113 
    114    double x[numBin-1];
    115    double y[numBin-1];
    116 
    117    for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
    118       {
    119         TAxis *xaxis = ETSoverET->GetXaxis();
    120         float binCenter = xaxis->GetBinCenter(i+1);
    121         int binMin = (int)xaxis->GetBinLowEdge(i+1);
    122         int binMax = (int)xaxis->GetBinUpEdge(i+1);
    123         cout<<"binMin "<<binMin<<" binMax "<<binMax<<" bin center "<<binCenter<<endl;
    124         char tempMin[500];
    125         sprintf(tempMin,"ETmisResol.Et>%d",binMin);
    126         string mystringMin(tempMin);
    127         char tempMax[500];
    128         sprintf(tempMax,"ETmisResol.Et<%d",binMax);
    129         string mystringMax(tempMax);
    130 
    131         char tempName[500];
    132         sprintf(tempName,"ETmisResol.EtSmeare>>hETSoverET%d",i);
    133         string mystringName(tempName);
    134 
    135         c5->cd(++frame);
    136         GaussValuesETmis(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
    137         x[i]=binCenter;
    138         y[i]=rms[i];
    139         cout<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl;
    140 
    141       }
    142 
    143    TCanvas *c6 = new TCanvas("c6","ETmis resolution",100,100,600,450);
    144    c6->cd();
    145 
    146    TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",0,200);
    147 
    148    TGraph *gr11 = new TGraph((numBin-1),x,rms);
    149    gr11->Draw("AP");
    150    gr11->GetXaxis()->SetTitle("MC^{MET} [GeV]");
    151    gr11->GetYaxis()->SetTitle("#sigma(REC^{MET}-MC^{MET})");
    152    gr11->Fit("user","RQ");
    153    gr11->Fit("user","RQ");
    154    gr11->Fit("user","RQ");
    155    gr11->Fit("user","R");
    156 
    157    TCanvas *c7 = new TCanvas("c7","ETmis resolution",100,100,600,450);
    158    c7->cd();
    159 
    160    TGraph *gr12 = new TGraph((numBin-1),x,mean);
    161    gr12->Draw("AP");
    162    gr12->GetXaxis()->SetTitle("MC^{MET} [GeV]");
    163    gr12->GetYaxis()->SetTitle("<(REC^{MET}-MC^{MET})>");
    164 
    165 
    166  delete fitfun;
     228 
     229  setTDRStyle();
     230  gROOT->Reset();
     231 
     232  TFile *f1 = new TFile("test2.root","read");
     233  TTree *Analyze = (TTree*)f1->Get("Analysis");
     234 
     235  TF1 *fitfun = new TF1("user","[0]*sqrt(x)",0,600);
     236  const Int_t numBin=7;
     237//  double bins[numBin]={10,100,150,200,250,300,350,400,450,500};
     238  double bins[numBin]={200,250,300,350,400,450,500};
     239//  double bins[numBin]={0,10,20,30,40,50,60,70};
     240  TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-1000,1000);
     241  Analyze->Draw("ETmisResol.ExSmeare:ETmisResol.SEt>>ETSoverET");
     242  Analyze->Fit("user","RQ");
     243  Analyze->Fit("user","RQ");
     244 
     245 
     246  double rms[numBin-1];
     247 
     248  TCanvas *c5 = new TCanvas("c5","PTmis resol",0,0,1000,650);
     249  c5->cd(); int frame=0;
     250  c5->Divide(6,2);
     251 
     252  double x[numBin];
     253  double y[numBin];
     254 
     255  for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
     256    {
     257      TAxis *xaxis = ETSoverET->GetXaxis();
     258      float binCenter = xaxis->GetBinCenter(i+1);
     259      int binMin = (int)xaxis->GetBinLowEdge(i+1);
     260      int binMax = (int)xaxis->GetBinUpEdge(i+1);
     261      char tempMin[500];
     262      sprintf(tempMin,"ETmisResol.SEt>%d",binMin);
     263      string mystringMin(tempMin);
     264      char tempMax[500];
     265      sprintf(tempMax,"ETmisResol.SEt<%d",binMax);
     266      string mystringMax(tempMax);
     267     
     268      char tempName[500];
     269      sprintf(tempName,"ETmisResol.ExSmeare>>hETSoverET%d",i);
     270      string mystringName(tempName);
     271     
     272      c5->cd(++frame);
     273      GaussValuesETmis(Analyze,tempName,rms[i],mystringMin,mystringMax);
     274      x[i]=binCenter;
     275      y[i]=rms[i];
     276     
     277    }
     278 
     279  TCanvas *c6 = new TCanvas("c6","ETmis resolution",100,100,600,450);
     280  c6->cd();
     281 
     282  x[numBin]=0;
     283  y[numBin]=0;
     284  TGraph *gr11 = new TGraph((numBin),x,rms);
     285  gr11->Draw("AP");
     286  gr11->GetXaxis()->SetTitle("Offline sum of E_{T} [GeV]");
     287  gr11->GetYaxis()->SetTitle("Resolution of x-component  of MET [GeV]");
     288  gr11->Fit("user","RQ");
     289  gr11->Fit("user","RQ");
     290  gr11->Fit("user","RQ");
     291  gr11->Fit("user","RQ");
     292  gr11->GetYaxis()->SetRangeUser(0,30);
     293  gr11->GetXaxis()->SetRangeUser(1,600);
     294
     295  Double_t* params = fitfun->GetParameters();
     296 
     297  char tempResol[500];
     298  sprintf(tempResol,"%f * #sqrt{E_{T}}",params[0]);
     299
     300  TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
     301  leg1->Draw();
     302
     303  TPaveText *leg2 = MakeTPave(0.2,0.8,0.8,0.85,"WHq'#rightarrow W#tau#tauq'#rightarrowjjl#tauq', m_{H}=150 GeV ");
     304  leg2->Draw();
     305
     306  TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
     307  Delphes->Draw();
     308 
     309 
     310  delete fitfun;
    167311}
    168312
    169313void General()
    170314{
    171   JetResol();
     315  //JetResol();
     316  //ElecResol();
    172317  ETmisResol();
    173 }
    174 
    175 
     318  //TauJetInfo();
     319 
     320}
     321
     322
Note: See TracChangeset for help on using the changeset viewer.