Index: trunk/routines/resolutions_atlas.C
===================================================================
--- trunk/routines/resolutions_atlas.C	(revision 473)
+++ trunk/routines/resolutions_atlas.C	(revision 479)
@@ -49,16 +49,27 @@
   
   //TFile *f1 = new TFile("JET2_atlas_jetclu.root","read");
+  //TFile *f1 = new TFile("JET2_atlas_jetclu_new.root","read"); // 40k
   //TFile *f1 = new TFile("JET2_atlas_siscone.root","read");
   //TFile *f1 = new TFile("JET2_atlas_antikt.root","read");
-  TFile *f1 = new TFile("JET2_atlas_midpoint.root","read");
+  //TFile *f1 = new TFile("JET2_atlas_midpoint.root","read");
   //TFile *f1 = new TFile("JET2_atlas_midpoint_eflow.root","read");
   //TFile *f1 = new TFile("JET2_atlas_midpoint_newCaloRes.root","read");
+  //TFile *f1 = new TFile("JET2_atlas_kt_40k.root","read"); // 40k events
+  TFile *f1 = new TFile("JET2_atlas_kt_60k.root","read"); // 60k events
+  //TFile *f1 = new TFile("JET2_atlas_jetclu_new_60k.root","read"); // 60k events
+
   if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
   if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
   TTree *Analyze = (TTree*)f1->Get("Analysis");
   
-  const Int_t numBin=16;
-  double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120,140,180,220,300,1200};
-  TProfile *ESoverE = new TProfile("ESoverE","Jet resolution for ATLAS ",(numBin-1),bins,-10,10);
+  //const int numBin=16;
+  //double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120,140,180,220,300,1200};
+  //const int numBin=25;
+  //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};
+  //const int numBin=20;
+  //double bins[numBin]={21,26,32,39,48, 59,72,88,107,130, 158,191,230,278,336, 405,485,581,696,834 };
+  const int numBin=14;
+  double bins[numBin]={21,26,32,39,48, 59,72,88,130, 191,278, 405,581,834 };
+  TProfile *ESoverE = new TProfile("ESoverE","Jet resolution for ATLAS ",numBin-1,bins,-10,10);
   
   double mean[numBin], mean2[numBin];
@@ -66,8 +77,8 @@
   TCanvas *c1 = new TCanvas("c1","JET resol: (E_reco - E_gen)/E_gen",0,0,1000,650);
   c1->cd(); int frame=0;
-  c1->Divide(6,2);
+  c1->Divide(10,2);
   TCanvas *c1b = new TCanvas("c1b","JET resol: [(E_reco - E_gen)/E_gen]^2 ",0,0,1000,650);
   c1b->cd(); int frame2=0;
-  c1b->Divide(6,2);
+  c1b->Divide(10,2);
   
   float x[numBin-1];
@@ -75,8 +86,16 @@
   float ex[numBin-1];
   float ey[numBin-1];
+  float erec_over_etruth[numBin-1];
   
   float finval=0;//valeur a remplir puis a fitter
-  
-  for ( int i=0; i<(numBin-1); i++)  // premiÃšre bin : i ==1 et pas i == 0
+  // available variables:
+  //  - E  : true Energy
+  //  - PT : true transverse energy
+  //  - SmearedPT : ratio ET_rec / ET_gen
+  //  - dE : ratio (E_rec - E_truth)/E_truth
+  //  - dE2 :  ( (E_rec - E_truth)/E_truth )^2
+
+ 
+  for ( int i=0; i<numBin-1; i++)  // premiÃšre bin : i ==1 et pas i == 0
     {
       TAxis *xaxis = ESoverE->GetXaxis();
@@ -92,29 +111,49 @@
       string mystringMax(tempMax);
       char tempName[500];
-      sprintf(tempName,"(JetPTResol.dE)>>hdE%d",i);
+      //sprintf(tempName,"JetPTResol.dE/JetPTResol.SmearedPT>>hdE%d",i);
+      //sprintf(tempName,"JetPTResol.dE>>hdE%d",i);
+      sprintf(tempName,"JetPTResol.dE_reco>>hdE%d",i);
       string mystringName(tempName);
       c1->cd(++frame);
-      GaussValuesAsymmetry(Analyze,tempName,mean[i],mystringMin,mystringMax);
-
-      sprintf(tempName,"(JetPTResol.dE2)>>hdE2%d",i);
+      //string cut = "abs(JetPTResol.Eta)<2.0 && abs(JetPTResol.Eta)>1.5";
+      string cut = "abs(JetPTResol.Eta)<0.5";
+      //string cut = "1==1";
+      GaussValuesAsymmetry(Analyze,tempName,mean[i],mystringMin,mystringMax,cut);
+
+      //sprintf(tempName,"JetPTResol.dE2/JetPTResol.SmearedPT>>hdE2%d",i);
+      //sprintf(tempName,"JetPTResol.dE2/JetPTResol.SmearedPT>>hdE2%d",i);
+      sprintf(tempName,"JetPTResol.dE2_reco>>hdE2%d",i);
       string mystringName2(tempName);
       c1b->cd(++frame2);
-      GaussValuesAsymmetry(Analyze,tempName,mean2[i],mystringMin,mystringMax);
+      GaussValuesAsymmetry2(Analyze,tempName,mean2[i],mystringMin,mystringMax,cut);
+
 
       x[i]=binCenter;
       finval=sqrt(mean2[i] - mean[i]*mean[i]);
-      y[i]=(finval*100);
+      y[i]=finval*100;
       ex[i]=0;         
       ey[i]=0;
+      erec_over_etruth[i]=mean[i]+1;
 
     }
   
+  TCanvas *c3 = new TCanvas("c3","JET linearity",100,100,600,450);
+  c3->cd();
+  //for(int i=0; i<numBin-1; i++) cout << x[i] << "\t\t" << erec_over_etruth[i] << endl;
+  TGraph *linearity = new TGraph(numBin-1,x,erec_over_etruth);
+  linearity->Draw("AP");
+  linearity->SetTitle("");
+  linearity->GetXaxis()->SetTitle("E^{Truth} [GeV]");
+  linearity->GetYaxis()->SetTitle("<E^{Reco}/E^{Truth}>");
+  linearity->GetXaxis()->SetRangeUser(30,1000);
+  linearity->SetMinimum(0.85);
+  linearity->SetMaximum(1.20);
+  gPad->SetLogx();
+
+
   TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
   c2->cd();
   
-  TF1 *fitfun = new TF1("user","sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))",10,800);
-  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);
-  
-  TGraph *gr11 = new TGraph((numBin-1),x,y);
+  TGraph *gr11 = new TGraph(numBin-1,x,y);
   gr11->Draw("AP");
   gr11->SetTitle("");
@@ -123,6 +162,6 @@
   gr11->GetYaxis()->SetTitle("#sigma(E)/E");
   
+  TF1 *fitfun = new TF1("user","sqrt(pow([0]/x,2)+pow([2]/sqrt(x),2)+pow([1],2))",10,800);
   Double_t* params = fitfun->GetParameters();
-
   fitfun->SetLineColor(kBlue);
   gr11->Fit("user","QRN");
@@ -130,26 +169,40 @@
   gr11->Fit("user","QRN");
   gr11->Fit("user","QRN");
+  gr11->GetXaxis()->SetRangeUser(30,1200);
+  //gPad->SetLogx();
   char tempResol[100];
-  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]);
+  //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]);
+  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]);
  char tempResol2[100];
-  sprintf(tempResol2,"sqrt(pow(%f/sqrt(x),2)+pow(%f,2))",params[1],params[2]);
+  sprintf(tempResol2,"sqrt(pow(%f/x,2)+pow(%f/sqrt(x),2)+pow(%f,2) )",params[0],params[2],params[1]);
 
 
   TF1 *fitfunDelphes = new TF1("userDelphes",tempResol2,7,1000);
-  fitfunDelphes->SetLineColor(kBlue);
+  fitfunDelphes->SetLineColor(kBlack);
   fitfunDelphes->SetLineStyle(7);
   fitfunDelphes->Draw("same");
 
-  fitfunATLAS->SetLineColor(kRed);
-  fitfunATLAS->SetLineWidth(2);
+  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);
+  fitfunATLAS->SetLineColor(kBlue);
   fitfunATLAS->Draw("same");
 
-  TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: pp #rightarrow gg ");
+  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);
+  fitfunATLAS2->SetLineColor(kBlue);
+  //fitfunATLAS2->Draw("same");
+
+  TPaveText *etacut = MakeTPave(0.58,0.51,0.65,0.55,"|#eta|<0.5 ");
+  etacut->Draw();
+
+  TPaveText *events = MakeTPave(0.62,0.43,0.75,0.47,"Events: pp #rightarrow ggX ");
   events->Draw();
 
-  TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
+  TPaveText *algo = MakeTPave(0.64,0.36,0.76,0.40,"k_{T} algorithm #Delta R = 0.6");
+  //TPaveText *algo = MakeTPave(0.64,0.36,0.76,0.40,"Cone algorithm #Delta R = 0.7");
+  algo->Draw();
+
+  TPaveText *Delphes = MakeTPave(0.62,0.29,0.85,0.33,"MG/ME + Pythia + Delphes");
   Delphes->Draw();
 
-  TLegend *legend = new TLegend(0.2,0.6,0.9,0.85,NULL,"NDC");
+  TLegend *legend = new TLegend(0.18,0.72,0.87,0.87,NULL,"NDC");
   legend->AddEntry(fitfunATLAS,"ATLAS resolution","l");
   legend->AddEntry(fitfunDelphes,tempResol,"l");
@@ -157,4 +210,5 @@
   legend->SetBorderSize(0);
   legend->Draw();
+
 
   delete fitfun;
