HcalCalibAlgos/macros/jet01.C

0001 {
0002 gROOT->Reset();
0003 gROOT->SetStyle("Plain");
0004
0005 gStyle->SetOptStat(1111);
0006 gStyle->SetOptFit(111);
0007
0008 Int_t n = 10;
0009 Double_t entries[n];
0010 TH1F *histResp[n];
0011
0012 histResp[0]  = new TH1F("h_resp_jf1_eta01_etgamma18_22"," ",50,0., 2.);
0013 histResp[1]  = new TH1F("h_resp_jf1_eta01_etgamma22_26"," ",50,0., 2.);
0014 histResp[2]  = new TH1F("h_resp_jf1_eta01_etgamma27_33"," ",50,0., 2.);
0015 histResp[3]  = new TH1F("h_resp_jf1_eta01_etgamma36_44"," ",50,0., 2.);
0016 histResp[4]  = new TH1F("h_resp_jf1_eta01_etgamma54_66"," ",50,0., 2.);
0017 histResp[5]  = new TH1F("h_resp_jf1_eta01_etgamma90_110"," ",50,0., 2.);
0018 histResp[6]  = new TH1F("h_resp_jf1_eta01_etgamma135_165"," ",50,0., 2.);
0019 histResp[7]  = new TH1F("h_resp_jf1_eta01_etgamma180_220"," ",50,0., 2.);
0020 histResp[8]  = new TH1F("h_resp_jf1_eta01_etgamma270_330"," ",50,0., 2.);
0021 histResp[9]  = new TH1F("h_resp_jf1_eta01_etgamma450_550"," ",50,0., 2.);
0022
0023 TH1F *histEtG[n];
0024 histEtG[0]  = new TH1F("h_etg_jf1_eta01_etgamma18_22"," ",100,0., 40.);
0025 histEtG[1]  = new TH1F("h_etg_jf1_eta01_etgamma22_26"," ",100,0., 50.);
0026 histEtG[2]  = new TH1F("h_etg_jf1_eta01_etgamma27_33"," ",100,0., 150.);
0027 histEtG[3]  = new TH1F("h_etg_jf1_eta01_etgamma36_44"," ",100,0., 150.);
0028 histEtG[4]  = new TH1F("h_etg_jf1_eta01_etgamma54_66"," ",100,20., 200.);
0029 histEtG[5]  = new TH1F("h_etg_jf1_eta01_etgamma90_110"," ",100,30.,300.);
0030 histEtG[6]  = new TH1F("h_etg_jf1_eta01_etgamma135_165"," ",100,50., 400.);
0031 histEtG[7]  = new TH1F("h_etg_jf1_eta01_etgamma180_220"," ",100,50., 600.);
0032 histEtG[8]  = new TH1F("h_etg_jf1_eta01_etgamma270_330"," ",100,100., 1000.);
0033 histEtG[9]  = new TH1F("h_etg_jf1_eta01_etgamma450_550"," ",100,100., 1500.);
0034
0035 TH1F *histEtJ[n];
0036 histEtJ[0]  = new TH1F("h_etj_jf1_eta01_etgamma18_22"," ",100,0., 100.);
0037 histEtJ[1]  = new TH1F("h_etj_jf1_eta01_etgamma22_26"," ",100,0., 100.);
0038 histEtJ[2]  = new TH1F("h_etj_jf1_eta01_etgamma27_33"," ",100,0., 150.);
0039 histEtJ[3]  = new TH1F("h_etj_jf1_eta01_etgamma36_44"," ",100,0., 150.);
0040 histEtJ[4]  = new TH1F("h_etj_jf1_eta01_etgamma54_66"," ",100,20., 200.);
0041 histEtJ[5]  = new TH1F("h_etj_jf1_eta01_etgamma90_110"," ",100,30.,300.);
0042 histEtJ[6]  = new TH1F("h_etj_jf1_eta01_etgamma135_165"," ",100,50., 400.);
0043 histEtJ[7]  = new TH1F("h_etj_jf1_eta01_etgamma180_220"," ",100,50., 600.);
0044 histEtJ[8]  = new TH1F("h_etj_jf1_eta01_etgamma270_330"," ",100,100., 1000.);
0045 histEtJ[9]  = new TH1F("h_etj_jf1_eta01_etgamma450_550"," ",100,100., 1500.);
0046
0047 // Histogram names needed for output filenemas
0048 char *name[n];
0049 name[0] = "JetResponseEt20Eta01";
0050 name[1] = "JetResponseEt24Eta01";
0051 name[2] = "JetResponseEt30Eta01";
0052 name[3] = "JetResponseEt40Eta01";
0053 name[4] = "JetResponseEt60Eta01";
0054 name[5] = "JetResponseEt100Eta01";
0055 name[6] = "JetResponseEt150Eta01";
0056 name[7] = "JetResponseEt200Eta01";
0057 name[8] = "JetResponseEt300Eta01";
0058 name[9] = "JetResponseEt500Eta01";
0059
0060 // Fit quantities
0061 TF1 *fit[n];
0062 Double_t resp[n], genEtG[n], respErr[n],genEtJ[n], genEtGErr[n], genEtJErr[n],respSigma[n], respSigmaErr[n] ;
0063 Double_t mean, rms, xmin, xmax;
0064 TH1F *histClone;
0065 TString *gif_file;
0066 TString *eps_file;
0067
0068 Int_t respPoints=0;
0069 Double_t pi=3.1415927;
0070
0071 // 20 GeV, Histo 0
0072 std::ifstream in20( "d20/1-01" );
0073 string line;
0074 entries[0] = 0;
0075 while( std::getline( in20, line)){
0076   istringstream linestream(line);
0077   double Etgamma,Etjet,Etajet;
0078   linestream>>Etgamma>>Etjet>>Etajet;
0079
0080   histResp[0]->Fill(Etjet/Etgamma);
0081   histEtG[0]->Fill(Etgamma);
0082   histEtJ[0]->Fill(Etjet);
0083   entries[0]++;
0084 }
0085
0086 // 24 GeV, Histo 1
0087 std::ifstream in24( "d24/1-01" );
0088 string line;
0089 entries[1] = 0;
0090 while( std::getline( in24, line)){
0091   istringstream linestream(line);
0092   double Etgamma,Etjet,Etajet;
0093   linestream>>Etgamma>>Etjet>>Etajet;
0094   histEtG[1]->Fill(Etgamma);
0095   histEtJ[1]->Fill(Etjet);
0096
0097   histResp[1]->Fill(Etjet/Etgamma);
0098   entries[1]++;
0099 }
0100
0101 // 30 GeV, Histo 2
0102 std::ifstream in30( "d30/1-01" );
0103 string line;
0104 entries[2] = 0;
0105 while( std::getline( in30, line)){
0106   istringstream linestream(line);
0107   double Etgamma,Etjet,Etajet;
0108   linestream>>Etgamma>>Etjet>>Etajet;
0109   histEtG[2]->Fill(Etgamma);
0110   histEtJ[2]->Fill(Etjet);
0111
0112   histResp[2]->Fill(Etjet/Etgamma);
0113   entries[2]++;
0114 }
0115 // 40 GeV, Histo 3
0116 std::ifstream in40( "d40/1-01" );
0117 string line;
0118 entries[3] = 0;
0119 while( std::getline( in40, line)){
0120   istringstream linestream(line);
0121   double Etgamma,Etjet,Etajet;
0122   linestream>>Etgamma>>Etjet>>Etajet;
0123   histEtG[3]->Fill(Etgamma);
0124   histEtJ[3]->Fill(Etjet);
0125
0126   histResp[3]->Fill(Etjet/Etgamma);
0127   entries[3]++;
0128 }
0129 // 60 GeV, Histo 4
0130 std::ifstream in60( "d60/1-01" );
0131 string line;
0132 entries[4] = 0;
0133 while( std::getline( in60, line)){
0134   istringstream linestream(line);
0135   double Etgamma,Etjet,Etajet;
0136   linestream>>Etgamma>>Etjet>>Etajet;
0137   histEtG[4]->Fill(Etgamma);
0138   histEtJ[4]->Fill(Etjet);
0139
0140   histResp[4]->Fill(Etjet/Etgamma);
0141   entries[4]++;
0142 }
0143
0144 // 100 GeV, Histo 5
0145 std::ifstream in100( "d100/1-01" );
0146 string line;
0147 entries[5] = 0;
0148 while( std::getline( in100, line)){
0149   istringstream linestream(line);
0150   double Etgamma,Etjet,Etajet;
0151   linestream>>Etgamma>>Etjet>>Etajet;
0152   histEtG[5]->Fill(Etgamma);
0153   histEtJ[5]->Fill(Etjet);
0154
0155   histResp[5]->Fill(Etjet/Etgamma);
0156   entries[5]++;
0157 }
0158 // 150 GeV, Histo 6
0159 std::ifstream in150( "d150/1-01" );
0160 string line;
0161 entries[6] = 0;
0162 while( std::getline( in150, line)){
0163   istringstream linestream(line);
0164   double Etgamma,Etjet,Etajet;
0165   linestream>>Etgamma>>Etjet>>Etajet;
0166   histEtG[6]->Fill(Etgamma);
0167   histEtJ[6]->Fill(Etjet);
0168
0169   histResp[6]->Fill(Etjet/Etgamma);
0170   entries[6]++;
0171 }
0172
0173
0174 // 200 GeV, Histo 7
0175 std::ifstream in200( "d200/1-01" );
0176 string line;
0177 entries[7] = 0;
0178 while( std::getline( in200, line)){
0179   istringstream linestream(line);
0180   double Etgamma,Etjet,Etajet;
0181   linestream>>Etgamma>>Etjet>>Etajet;
0182   histEtG[7]->Fill(Etgamma);
0183   histEtJ[7]->Fill(Etjet);
0184
0185   histResp[7]->Fill(Etjet/Etgamma);
0186   entries[7]++;
0187 }
0188
0189 // 300 GeV, Histo 8
0190 std::ifstream in300( "d300/1-01" );
0191 string line;
0192 entries[8] = 0;
0193 while( std::getline( in300, line)){
0194   istringstream linestream(line);
0195   double Etgamma,Etjet,Etajet;
0196   linestream>>Etgamma>>Etjet>>Etajet;
0197   histEtG[8]->Fill(Etgamma);
0198   histEtJ[8]->Fill(Etjet);
0199
0200   histResp[8]->Fill(Etjet/Etgamma);
0201   entries[8]++;
0202 }
0203 // 500 GeV, Histo 9
0204 std::ifstream in500( "d500/1-01" );
0205 string line;
0206 entries[9] = 0;
0207 while( std::getline( in500, line)){
0208   istringstream linestream(line);
0209   double Etgamma,Etjet,Etajet;
0210   linestream>>Etgamma>>Etjet>>Etajet;
0211   histEtG[9]->Fill(Etgamma);
0212   histEtJ[9]->Fill(Etjet);
0213
0214   histResp[9]->Fill(Etjet/Etgamma);
0215   entries[9]++;
0216 }
0217
0218 for (Int_t i=0; i<10;i++)
0219 {
0220   cout<<" Number of measurements "<<entries[i]<< " at point "<<i<<endl;
0221 }
0222 // Draw the second graphic with the three continuous functions
0223 TCanvas c1("c1"," ",10,10,800,600);
0224
0225 Int_t entryCut = 30;
0226 Int_t respPoints = 0;
0227
0228 for (Int_t i=0; i<n; i++)
0229 {
0230   cout<<" Here "<<i<<endl;
0231   entries[respPoints] = histResp[respPoints]->GetEntries();
0232   if(entries[respPoints] > entryCut)
0233   {
0234
0235     fit[i] = new TF1("mygaus","[0]*exp(-0.5*((x-[1])/[2])^2)");
0236     //Make the points red with error bars and save a copy to overlay
0237     histResp[i]->SetLineWidth(2);
0238     histResp[i]->SetLineColor(2);
0239     histClone = (TH1F*)histResp[i]->Clone();
0240     //1st fit using 1 sigma range from full dist mean and rms
0241     mean=histResp[i]->GetMean();
0242     rms=histResp[i]->GetRMS();
0243     xmin=mean-rms;
0244     xmax=mean+rms;
0245     cout << "Full dist: mean=" << mean <<", rms="<<rms<<", xmin="<<xmin<<", xmax="<<xmax<<endl;
0246     Double_t binWidth=histResp[i]->GetBinWidth(1);
0247     fit[i]->SetParameter(0,binWidth*entries[i]/(rms*sqrt(2*pi)));
0248     fit[i]->SetParLimits(0,0,entries[i]);
0249     fit[i]->SetParameter(1,mean);
0250     fit[i]->SetParLimits(1,xmin,xmax);
0251     fit[i]->SetParameter(2,rms);
0252     fit[i]->SetParLimits(2,0,1);
0253     histResp[i]->Fit("mygaus","","",xmin,xmax);
0254     //fit[i] = histResp[i]->GetFunction("gaus");
0255     mean = fit[i]->GetParameter(1);
0256     rms = fit[i]->GetParameter(2);
0257
0258     //2nd fit using 1 sigma range from last Gaussian Fit
0259     xmin=mean-rms;
0260     xmax=mean+rms;
0261     cout << "1st Gaus Fit: mean=" << mean <<", rms="<<rms<<", xmin="<<xmin<<", xmax="<<xmax<<endl;
0262     fit[i]->SetParLimits(1,xmin,xmax);
0263     histResp[i]->Fit("mygaus","","",xmin,xmax);
0264     //fit[i] = histResp[i]->GetFunction("gaus");
0265     mean = fit[i]->GetParameter(1);
0266     rms = fit[i]->GetParameter(2);
0267     //Final fit using 1 sigma range from last Gaussian Fit
0268     xmin=mean-rms;
0269     xmax=mean+rms;
0270     cout << "2nd Gaus Fit: mean=" << mean <<", rms="<<rms<<", xmin="<<xmin<<", xmax="<<xmax<<endl;
0271     fit[i]->SetParLimits(1,xmin,xmax);
0272     histResp[i]->Fit("mygaus","","e",xmin,xmax);
0273     histClone->Draw("eSAME");
0274     //fit[i] = histResp[i]->GetFunction("gaus");
0275     mean = fit[i]->GetParameter(1);
0276     rms = fit[i]->GetParameter(2);
0277     cout << "Final Gaus Fit: mean=" << mean <<", rms="<<rms<<endl;
0278     if( i == 0 )
0279     {
0280       histResp[i]->Fit("mygaus","","",0.15,0.6);
0281     }
0282     if( i == 1 )
0283     {
0284       histResp[i]->Fit("mygaus","","",0.15,0.6);
0285     }
0286
0287     if( i == 2 )
0288     {
0289       histResp[i]->Fit("mygaus","","",0.1,0.55);
0290     }
0291
0292     if( i == 4 )
0293     {
0294       histResp[i]->Fit("mygaus","","",0.15,0.9);
0295     }
0296
0297
0298     histResp[i]->Draw();
0299     gif_file = new TString("tmp/");
0300     *gif_file+= name[i];
0301     *gif_file+= ".gif";
0302     c1.Print(*gif_file);
0303 // Try to find curve
0304
0305
0306     //Get the genJetEt, response, and their errors
0307     genEtG[respPoints]=histEtG[i]->GetMean();
0308     cout<<" Mean "<<histEtG[i]->GetMean()<<endl;
0309     genEtGErr[respPoints]=histEtG[i]->GetRMS()/sqrt(entries[i]);
0310
0311     genEtJ[respPoints]=histEtJ[i]->GetMean();
0312     cout<<" Mean "<<histEtJ[i]->GetMean()<<endl;
0313     genEtJErr[respPoints]=histEtJ[i]->GetRMS()/sqrt(entries[i]);
0314
0315     resp[respPoints]=fit[i]->GetParameter(1);
0316     respErr[respPoints]=fit[i]->GetParError(1);
0317     respSigma[respPoints]=fit[i]->GetParameter(2);
0318     respSigmaErr[respPoints]=fit[i]->GetParError(2);
0319     respPoints++;
0320    }
0321 }
0322 // Save usefule fit quantities for later
0323
0324 Float_t k[10];
0325 ofstream outFile("Response_01.dat");
0326 for (Int_t i=0; i<respPoints; i++){
0327   k[i] =  genEtJ[i]/genEtG[i];
0328   outFile << " EtG  " << genEtG[i] <<" EtJ "<<genEtJ[i] <<" K "<<k[i]<< "  " << resp[i] << "  " << respErr[i] << " " << respSigma[i] << "  " <<  respSigmaErr[i] <<endl;
0329 }
0330 // Two graphics, the first is for the fit, the second we will save
0331 TGraphErrors* gr1 = new TGraphErrors (respPoints,genEtJ,resp,genEtJErr,respErr);
0332 TGraphErrors* gr2 = new TGraphErrors (respPoints,genEtG,resp,genEtGErr,respErr);
0333 TGraphErrors* gr3 = new TGraphErrors (respPoints,genEtJ,resp,genEtJErr,respErr);
0334 /*
0335 Double_t par[5];
0336 TF1 *func1 = new TF1("func1","[0]+[1]*log(x+[2])-[3]/(x+20.)",10.,600.);
0337 func1->SetParameter(0,1.);
0338 func1->SetParameter(1,1.);
0339 func1->SetParameter(2,1.);
0340 func1->SetParameter(3,1.);
0341 //func1->SetParameter(4,10.);
0342 gr1->Fit("func1", "R",  "r", 10., 600.);
0343 TF1 *fitM1 = gr1->GetFunction("func1");
0344 func1->GetParameters(&par[0]);
0345
0346 // Write out the parameters of the two functions which define the response.
0347 Float_t etamax = 0.226;
0348
0349 FILE *Out1 = fopen("h01.txt", "w+");
0350
0351 fprintf(Out1," %.5f %.5f %.5f %.5f %.5f\n",
0352 etamax, par[0], par[1], par[2], par[3]);
0353
0354 fclose(Out1);
0355 */
0356 Double_t par[7];
0357
0358 Double_t por1, por2;
0359
0360 por1 = 18.;
0361 por2 = 23.;
0362
0363 // The three functions that define the response.
0364 // The first two functions will be fit to the data at low and high pt.
0365 // The third function (line) will connect the first two functions together continuously in Et and Response.
0366 TF1 *func1 = new TF1("func1","[0]*sqrt(x +[1]) + [2]",10.,por2);
0367 TF1 *func2 = new TF1("func2","[0]/(sqrt([1]*x + [2])) + [3] ",por1,5000.);
0368 TF1 *func3 = new TF1("func3","pol1 ",por1,por2);
0369
0370 func2->SetParameter(0,-3.);
0371 func2->SetParameter(3,1.);
0372 func2->SetParLimits(1,0.1,10.);
0373 func2->SetParLimits(2,-1.5,10000.);
0374
0375 // Fit the response graphic with the first two functions. Fits overlap in region por1 < Et < por2 !!
0376 gr1->Fit("func2", "R",  "r", por1, 5000);
0377 TF1 *fitM1 = gr1->GetFunction("func2");
0378
0379 gr1->Fit("func1", "R+", "r", 10, por2);
0380 TF1 *fitM2 = gr1->GetFunction("func1");
0381
0382 // Load the fit parameters into a single array
0383 Double_t par[7];
0384 func1->GetParameters(&par[0]);
0385 func2->GetParameters(&par[3]);
0386
0387 cout<<"  "<< par[0]<<"  "<< par[1]<<"  "<< par[2]
0388 <<"  "<< par[3]<<"  "<< par[4]<<"  "<< par[5]
0389 <<"  "<< par[6] << endl;
0390
0391 // Calculate the line that connects the two functions in the region por1 < Et < por2
0392 Double_t a1,a2;
0393
0394 a2 = par[3]/(sqrt(fabs(par[4]*por2 + par[5]))) + par[6];
0395 a1 = par[0]*sqrt(por1 +par[1]) + par[2];
0396
0397 func3->SetParameter(0,(a1*por2 - a2*por1)/(por2-por1));
0398 func3->SetParameter(1,(a2-a1)/(por2-por1));
0399
0400 func3->Draw("same");
0401
0402 cout<<"  "<< a1 <<"  "<< a2 << endl;
0403
0404 // Redefine the first two functions so they cover Et<Por1 and Et>por2 respectively.
0405 TF1 *funcA1 = new TF1("funcA1","[0]*sqrt(x +[1]) + [2]",10.,por1);
0406 TF1 *funcA2 = new TF1("funcA2","[0]/(sqrt([1]*x + [2])) + [3] ",por2,5000.);
0407
0408 funcA1->SetParameter(0,par[0]);
0409 funcA1->SetParameter(1,par[1]);
0410 funcA1->SetParameter(2,par[2]);
0411
0412 funcA2->SetParameter(0,par[3]);
0413 funcA2->SetParameter(1,par[4]);
0414 funcA2->SetParameter(2,par[5]);
0415 funcA2->SetParameter(3,par[6]);
0416
0417
0418
0419 // Write out the parameters of the two functions which define the response.
0420 Float_t etamax = 0.226;
0421
0422 FILE *Out1 = fopen("h01.txt", "w+");
0423
0424 fprintf(Out1," %.3f %.1f %.1f %.5f %.5f %.5f %.5f %.5f %.5f %.5f\n",
0425 etamax, por1, por2, par[0], par[1], par[2], par[3], par[4], par[5], par[6]);
0426
0427 fclose(Out1);
0428
0429
0430
0431 // Draw the second graphic with the three continuous functions
0432
0433 TCanvas c2("c2"," ",10,10,800,600);
0434 c2->SetLogx();
0435 gPad->SetTicks(1,1);
0436 gr2->SetMaximum(1.0);
0437 gr2->SetMinimum(0.0);
0438 gr2->GetYaxis()->SetTitle("Jet Response");
0439 gr2->GetYaxis()->SetTitleSize(0.05);
0440 gr2->GetXaxis()->SetTitle("GammaJet E_{T} (GeV)");
0441 gr2->GetXaxis()->SetTitleOffset(1.2);
0442 gr2->GetXaxis()->SetRangeUser(0.95*genEtG[0],1.05*genEtG[respPoints]);
0443 gr2->Draw("AP");
0444 c2.Print("GammaJetEta01_1.gif");
0445
0446 TCanvas c3("c3"," ",10,10,800,600);
0447 c3->SetLogx();
0448 gPad->SetTicks(1,1);
0449 gr3->SetMaximum(1.0);
0450 gr3->SetMinimum(0.0);
0451 gr3->GetYaxis()->SetTitle("Jet Response");
0452 gr3->GetYaxis()->SetTitleSize(0.05);
0453 gr3->GetXaxis()->SetTitle("RecJet E_{T} (GeV)");
0454 gr3->GetXaxis()->SetTitleOffset(1.2);
0455 gr3->GetXaxis()->SetRangeUser(0.95*genEtJ[0],1.05*genEtJ[respPoints]);
0456 gr3->Draw("AP");
0457 //func1->Draw("same");
0458 funcA1->Draw("same");
0459 funcA2->Draw("same");
0460 func3->Draw("same");
0461
0462 c3.Print("GammaJetEta01_2.gif");
0463
0464 }