Fork me on GitHub

Ignore:
Timestamp:
Jul 14, 2009, 12:51:51 PM (15 years ago)
Author:
Xavier Rouby
Message:

update, better&nicer plots

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/routines/resolutions_atlas.C

    r473 r479  
    4949 
    5050  //TFile *f1 = new TFile("JET2_atlas_jetclu.root","read");
     51  //TFile *f1 = new TFile("JET2_atlas_jetclu_new.root","read"); // 40k
    5152  //TFile *f1 = new TFile("JET2_atlas_siscone.root","read");
    5253  //TFile *f1 = new TFile("JET2_atlas_antikt.root","read");
    53   TFile *f1 = new TFile("JET2_atlas_midpoint.root","read");
     54  //TFile *f1 = new TFile("JET2_atlas_midpoint.root","read");
    5455  //TFile *f1 = new TFile("JET2_atlas_midpoint_eflow.root","read");
    5556  //TFile *f1 = new TFile("JET2_atlas_midpoint_newCaloRes.root","read");
     57  //TFile *f1 = new TFile("JET2_atlas_kt_40k.root","read"); // 40k events
     58  TFile *f1 = new TFile("JET2_atlas_kt_60k.root","read"); // 60k events
     59  //TFile *f1 = new TFile("JET2_atlas_jetclu_new_60k.root","read"); // 60k events
     60
    5661  if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
    5762  if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
    5863  TTree *Analyze = (TTree*)f1->Get("Analysis");
    5964 
    60   const Int_t numBin=16;
    61   double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120,140,180,220,300,1200};
    62   TProfile *ESoverE = new TProfile("ESoverE","Jet resolution for ATLAS ",(numBin-1),bins,-10,10);
     65  //const int numBin=16;
     66  //double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120,140,180,220,300,1200};
     67  //const int numBin=25;
     68  //double bins[numBin]={21,26,32,39,48, 59,72,88,107,130, 158,191,230,278,336, 405,485,581,696,834, 1000,1198,1435,1719,2060};
     69  //const int numBin=20;
     70  //double bins[numBin]={21,26,32,39,48, 59,72,88,107,130, 158,191,230,278,336, 405,485,581,696,834 };
     71  const int numBin=14;
     72  double bins[numBin]={21,26,32,39,48, 59,72,88,130, 191,278, 405,581,834 };
     73  TProfile *ESoverE = new TProfile("ESoverE","Jet resolution for ATLAS ",numBin-1,bins,-10,10);
    6374 
    6475  double mean[numBin], mean2[numBin];
    … …  
    6677  TCanvas *c1 = new TCanvas("c1","JET resol: (E_reco - E_gen)/E_gen",0,0,1000,650);
    6778  c1->cd(); int frame=0;
    68   c1->Divide(6,2);
     79  c1->Divide(10,2);
    6980  TCanvas *c1b = new TCanvas("c1b","JET resol: [(E_reco - E_gen)/E_gen]^2 ",0,0,1000,650);
    7081  c1b->cd(); int frame2=0;
    71   c1b->Divide(6,2);
     82  c1b->Divide(10,2);
    7283 
    7384  float x[numBin-1];
    … …  
    7586  float ex[numBin-1];
    7687  float ey[numBin-1];
     88  float erec_over_etruth[numBin-1];
    7789 
    7890  float finval=0;//valeur a remplir puis a fitter
    79  
    80   for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
     91  // available variables:
     92  //  - E  : true Energy
     93  //  - PT : true transverse energy
     94  //  - SmearedPT : ratio ET_rec / ET_gen
     95  //  - dE : ratio (E_rec - E_truth)/E_truth
     96  //  - dE2 :  ( (E_rec - E_truth)/E_truth )^2
     97
     98 
     99  for ( int i=0; i<numBin-1; i++)  // premiÚre bin : i ==1 et pas i == 0
    81100    {
    82101      TAxis *xaxis = ESoverE->GetXaxis();
    … …  
    92111      string mystringMax(tempMax);
    93112      char tempName[500];
    94       sprintf(tempName,"(JetPTResol.dE)>>hdE%d",i);
     113      //sprintf(tempName,"JetPTResol.dE/JetPTResol.SmearedPT>>hdE%d",i);
     114      //sprintf(tempName,"JetPTResol.dE>>hdE%d",i);
     115      sprintf(tempName,"JetPTResol.dE_reco>>hdE%d",i);
    95116      string mystringName(tempName);
    96117      c1->cd(++frame);
    97       GaussValuesAsymmetry(Analyze,tempName,mean[i],mystringMin,mystringMax);
    98 
    99       sprintf(tempName,"(JetPTResol.dE2)>>hdE2%d",i);
     118      //string cut = "abs(JetPTResol.Eta)<2.0 && abs(JetPTResol.Eta)>1.5";
     119      string cut = "abs(JetPTResol.Eta)<0.5";
     120      //string cut = "1==1";
     121      GaussValuesAsymmetry(Analyze,tempName,mean[i],mystringMin,mystringMax,cut);
     122
     123      //sprintf(tempName,"JetPTResol.dE2/JetPTResol.SmearedPT>>hdE2%d",i);
     124      //sprintf(tempName,"JetPTResol.dE2/JetPTResol.SmearedPT>>hdE2%d",i);
     125      sprintf(tempName,"JetPTResol.dE2_reco>>hdE2%d",i);
    100126      string mystringName2(tempName);
    101127      c1b->cd(++frame2);
    102       GaussValuesAsymmetry(Analyze,tempName,mean2[i],mystringMin,mystringMax);
     128      GaussValuesAsymmetry2(Analyze,tempName,mean2[i],mystringMin,mystringMax,cut);
     129
    103130
    104131      x[i]=binCenter;
    105132      finval=sqrt(mean2[i] - mean[i]*mean[i]);
    106       y[i]=(finval*100);
     133      y[i]=finval*100;
    107134      ex[i]=0;         
    108135      ey[i]=0;
     136      erec_over_etruth[i]=mean[i]+1;
    109137
    110138    }
    111139 
     140  TCanvas *c3 = new TCanvas("c3","JET linearity",100,100,600,450);
     141  c3->cd();
     142  //for(int i=0; i<numBin-1; i++) cout << x[i] << "\t\t" << erec_over_etruth[i] << endl;
     143  TGraph *linearity = new TGraph(numBin-1,x,erec_over_etruth);
     144  linearity->Draw("AP");
     145  linearity->SetTitle("");
     146  linearity->GetXaxis()->SetTitle("E^{Truth} [GeV]");
     147  linearity->GetYaxis()->SetTitle("<E^{Reco}/E^{Truth}>");
     148  linearity->GetXaxis()->SetRangeUser(30,1000);
     149  linearity->SetMinimum(0.85);
     150  linearity->SetMaximum(1.20);
     151  gPad->SetLogx();
     152
     153
    112154  TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
    113155  c2->cd();
    114156 
    115   TF1 *fitfun = new TF1("user","sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))",10,800);
    116   TF1 *fitfunATLAS = new TF1("userATLAS","sqrt(pow(1.03*100/sqrt(x),2)+pow(0.026*100,2)+pow(8*100/x,2))",10,800);
    117  
    118   TGraph *gr11 = new TGraph((numBin-1),x,y);
     157  TGraph *gr11 = new TGraph(numBin-1,x,y);
    119158  gr11->Draw("AP");
    120159  gr11->SetTitle("");
    … …  
    123162  gr11->GetYaxis()->SetTitle("#sigma(E)/E");
    124163 
     164  TF1 *fitfun = new TF1("user","sqrt(pow([0]/x,2)+pow([2]/sqrt(x),2)+pow([1],2))",10,800);
    125165  Double_t* params = fitfun->GetParameters();
    126 
    127166  fitfun->SetLineColor(kBlue);
    128167  gr11->Fit("user","QRN");
    … …  
    130169  gr11->Fit("user","QRN");
    131170  gr11->Fit("user","QRN");
     171  gr11->GetXaxis()->SetRangeUser(30,1200);
     172  //gPad->SetLogx();
    132173  char tempResol[100];
    133   sprintf(tempResol,"Delphes resolution: #frac{#sigma(E)}{E} = #sqrt{ <( #frac{E_{reco} - E_{gen}}{E_{gen}} )^{2} > - < #frac{E_{reco}-E_{gen}}{E_{gen}}>^{2} } =\n #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus %f #oplus #frac{%f}{E}",params[1],params[2],params[3]);
     174  //sprintf(tempResol,"Delphes resolution: #frac{#sigma(E)}{E} = #sqrt{ <( #frac{E_{reco} - E_{gen}}{E_{gen}} )^{2} > - < #frac{E_{reco}-E_{gen}}{E_{gen}}>^{2} } =\n #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus %f #oplus #frac{%f}{E}",params[1],params[2],params[3]);
     175  sprintf(tempResol,"Delphes resolution: #frac{#sigma(E)}{E} = #frac{%.1f}{#sqrt{E^{MC}}} #oplus %.1f #oplus #frac{%.1f}{E^{MC}}",params[2],params[1],params[0]);
    134176 char tempResol2[100];
    135   sprintf(tempResol2,"sqrt(pow(%f/sqrt(x),2)+pow(%f,2))",params[1],params[2]);
     177  sprintf(tempResol2,"sqrt(pow(%f/x,2)+pow(%f/sqrt(x),2)+pow(%f,2) )",params[0],params[2],params[1]);
    136178
    137179
    138180  TF1 *fitfunDelphes = new TF1("userDelphes",tempResol2,7,1000);
    139   fitfunDelphes->SetLineColor(kBlue);
     181  fitfunDelphes->SetLineColor(kBlack);
    140182  fitfunDelphes->SetLineStyle(7);
    141183  fitfunDelphes->Draw("same");
    142184
    143   fitfunATLAS->SetLineColor(kRed);
    144   fitfunATLAS->SetLineWidth(2);
     185  TF1 *fitfunATLAS = new TF1("userATLAS","sqrt(pow(0.69*100/sqrt(x),2)+pow(0.03*100,2)+pow(6.3*100/x,2))",10,800);
     186  fitfunATLAS->SetLineColor(kBlue);
    145187  fitfunATLAS->Draw("same");
    146188
    147   TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: pp #rightarrow gg ");
     189  TF1 *fitfunATLAS2= new TF1("userATLAS2","sqrt(pow(1.19*100/sqrt(x),2)+pow(0.01*100,2)+pow(10*100/x,2))",10,800);
     190  fitfunATLAS2->SetLineColor(kBlue);
     191  //fitfunATLAS2->Draw("same");
     192
     193  TPaveText *etacut = MakeTPave(0.58,0.51,0.65,0.55,"|#eta|<0.5 ");
     194  etacut->Draw();
     195
     196  TPaveText *events = MakeTPave(0.62,0.43,0.75,0.47,"Events: pp #rightarrow ggX ");
    148197  events->Draw();
    149198
    150   TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
     199  TPaveText *algo = MakeTPave(0.64,0.36,0.76,0.40,"k_{T} algorithm #Delta R = 0.6");
     200  //TPaveText *algo = MakeTPave(0.64,0.36,0.76,0.40,"Cone algorithm #Delta R = 0.7");
     201  algo->Draw();
     202
     203  TPaveText *Delphes = MakeTPave(0.62,0.29,0.85,0.33,"MG/ME + Pythia + Delphes");
    151204  Delphes->Draw();
    152205
    153   TLegend *legend = new TLegend(0.2,0.6,0.9,0.85,NULL,"NDC");
     206  TLegend *legend = new TLegend(0.18,0.72,0.87,0.87,NULL,"NDC");
    154207  legend->AddEntry(fitfunATLAS,"ATLAS resolution","l");
    155208  legend->AddEntry(fitfunDelphes,tempResol,"l");
    … …  
    157210  legend->SetBorderSize(0);
    158211  legend->Draw();
     212
    159213
    160214  delete fitfun;
Note: See TracChangeset for help on using the changeset viewer.