HcalCalibAlgos/macros/extract_pfcorrs.C

0001
0002 {
0003   #include <map>
0004
0005   gStyle -> SetOptStat(0);
0006 //  gStyle->SetPalette(1);
0007   gStyle->SetPadGridX(1);
0008   gStyle->SetPadGridY(1);
0009
0010 TString imgpath("~/afs/public_html/pfcorrs/v02/");
0011
0012  TFile* fcorrs = new TFile("~/nobackup/arch/hcalCorrsFile_resp01.root","OPEN");
0013  TProfile* respcorrs = (TProfile*)fcorrs->Get("corrs1");
0014
0015  TFile* tf = new TFile("~/nobackup/arch/hcalCorrPFv13_35cm.root","OPEN");
0016  //TFile* tf = new TFile("~/afs/arch/hcalCorrPFv2_26.2cm.root","OPEN");
0017 //TFile* tf = new TFile("./HcalCorrPF.root","OPEN");
0018  ofstream new_pfcorrs("newPFcorrs.txt");
0019  ofstream new_pfcorrs_subnoise("newPFcorrsNoiseSubtracted.txt");
0020
0021 // ofstream new_pfcorrs("newPFcorrs26.2cm.txt");
0022 // ofstream new_pfcorrs_subnoise("newPFcorrsNoiseSubtracted26.2cm.txt");
0023
0024
0025 /*
0026  TString imgpath("~/afs/public_html/validation/pfcorrs/30cm/");
0027  TFile* tf = new TFile("./HcalCorrPFcone30cm.root","OPEN");
0028  ofstream new_pfcorrs("newPFcorrs30cm.txt");
0029  ofstream new_pfcorrs_subnoise("newPFcorrsNoiseSubtracted30cm.txt");
0030 */
0031 /*
0032  TString imgpath("~/afs/public_html/validation/pfcorrs/40cm/");
0033  TFile* tf = new TFile("./HcalCorrPFcone40cm.root","OPEN");
0034  ofstream new_pfcorrs("newPFcorrs40cm.txt");
0035  ofstream new_pfcorrs_subnoise("newPFcorrsNoiseSubtracted.txt");
0036   */
0037
0038 /*
0039  TString imgpath("~/afs/public_html/validation/pfcorrs/50cm/");
0040  TFile* tf = new TFile("./HcalCorrPFcone50cm.root","OPEN");
0041  ofstream new_pfcorrs("newPFcorrs50cm.txt");
0042  ofstream new_pfcorrs_subnoise("newPFcorrsNoiseSubtracted50cm.txt");
0043 */
0044
0045 /*
0046  TProfile* p1 = (TProfile*)tf->Get("hcalRecoAnalyzer/enHcal");
0047  TProfile* p2 = (TProfile*)tf->Get("hcalRecoAnalyzer/enHcalNoise");
0048  TProfile* p3 = (TProfile*)tf->Get("hcalRecoAnalyzer/nCells");
0049  TProfile* p4 = (TProfile*)tf->Get("hcalRecoAnalyzer/nCellsNoise");
0050  TProfile* p5 = (TProfile*)tf->Get("hcalRecoAnalyzer/enEcal");
0051 */
0052
0053  TCut tr_cut = "eParticle>45 && eParticle<55";  label = new TText(0.03,0.2, "50 GeV");
0054
0055  label -> SetNDC(); label->SetTextAlign(11); label->SetTextSize(0.05); label->SetTextAngle(90); label->SetTextColor(kRed);
0056
0057   TCut trkQual = "";
0058   //  TCut trkQual = "(trkQual[0]==1 && abs(etaParticle)<2.4) || abs(etaParticle)>2.4";
0059   TCut mip_cut = "eECAL09cm<1.";
0060     TCut hit_dist = "delR<16";
0061   TCut mc_dist = "";
0062   // TCut mc_dist = "(abs(etaParticle)<2.4 && delRmc[0]<2) || abs(etaParticle)>2.4";
0063   //TCut low_resp_cut = "";
0064   TCut low_resp_cut = "eHcalCone/eParticle>0.2";
0065   TCut maxPNearBy = "";
0066   //  TCut neutral_iso_cut ="";
0067   TCut neutral_iso_cut = "(abs(iEta)<=13&&(eECAL40cm-eECAL09cm)<6.8)||(abs(iEta)==14&&(eECAL40cm-eECAL09cm)<6.6)||(abs(iEta)==15&&(eECAL40cm-eECAL09cm)<6.6)||(abs(iEta)==16&&(eECAL40cm-eECAL09cm)<9.8)||(abs(iEta)==17&&(eECAL40cm-eECAL09cm)<10.4)||(abs(iEta)==18&&(eECAL40cm-eECAL09cm)<9.8)||(abs(iEta)==19&&(eECAL40cm-eECAL09cm)<11.0)||(abs(iEta)==20&&(eECAL40cm-eECAL09cm)<12.3)||(abs(iEta)==21&&(eECAL40cm-eECAL09cm)<13.6)||(abs(iEta)==22&&(eECAL40cm-eECAL09cm)<15.2)||(abs(iEta)==23&&(eECAL40cm-eECAL09cm)<15.4)||(abs(iEta)==24&&(eECAL40cm-eECAL09cm)<16.3)||(abs(iEta)==25&&(eECAL40cm-eECAL09cm)<16.1)||(abs(iEta)==26&&(eECAL40cm-eECAL09cm)<15.4)||(abs(iEta)==27&&(eECAL40cm-eECAL09cm)<15.4)||abs(iEta)>27";
0068
0069   TCut selection = trkQual && neutral_iso_cut && tr_cut && mip_cut && hit_dist && mc_dist && maxPNearBy && low_resp_cut;
0070
0071
0072   TTree* pftree = (TTree*)tf->Get("hcalPFcorrs/pfTree");
0073 //  pftree -> Draw("eHcalCone/eParticle:iEta>>p1(82, -40.5, 41.5)", "", "prof goff");
0074   pftree -> Draw("eHcalCone/eParticle:iEta>>p1(82, -40.5, 41.5)", selection, "prof goff");
0075   pftree -> Draw("(eHcalCone-eHcalConeNoise)/eParticle:iEta>>p2(82, -40.5, 41.5)", selection, "prof goff");
0076   pftree -> Draw("UsedCells:iEta>>p3(82, -40.5, 41.5)", selection, "prof goff");
0077   pftree -> Draw("UsedCellsNoise:iEta>>p4(82, -40.5, 41.5)", selection, "prof goff");
0078   pftree -> Draw("eECAL:iEta>>p5(82, -40.5, 41.5)", selection, "prof goff");
0079
0080 TProfile* diff1 = new TProfile("diff1", "(oldvalue - new)/oldvalue", 84,-42,42);
0081 TProfile* diff2 = new TProfile("diff2", "(pfcorr - pfSubNoise)/pfcorr", 84,-42,42);
0082
0083
0084 TProfile* corrs2 = new TProfile("corrs2", "PF corrs", 84,-42,42);
0085 TProfile* corrs3 = new TProfile("corrs3", "PF corrs", 84,-42,42);
0086 TProfile* corrs4 = new TProfile("corrs4", "PF corrs", 84,-42,42);
0087
0088 //ifstream old_pfcorrs("../data/HcalPFCorrs_v1.03_mc.txt");
0089 //ifstream old_pfcorrs("../data/HcalPFCorrs_v2.00_mc.txt");
0090 ifstream old_pfcorrs("../data/response_corrections.txt");
0091
0092
0093   //FILE *new_pfcorrs = fopen("new_pfcorrs.txt", "w+");
0094   Int_t   iEta;
0095   UInt_t  iPhi;
0096   Int_t  depth;
0097   //TString sdName;
0098   string sdName;
0099   UInt_t  detId;
0100   Float_t value;
0101 //  HcalSubdetector sd;
0102   map<Int_t, Float_t> CorrValues;
0103   map<Int_t, Float_t> CorrValuesSubNoise;
0104
0105 for (Int_t i=0; i<=82; i++)
0106 {
0107
0108   if (p1->GetBinContent(i)!=0) CorrValues[p1->GetBinCenter(i)]= 1./p1->GetBinContent(i);
0109   if (p2->GetBinContent(i)!=0) CorrValuesSubNoise[p2->GetBinCenter(i)]= 1./p2->GetBinContent(i);
0110
0111   cout<<"bin: "<< p1->GetBinCenter(i)<<"   response: "<< p1->GetBinContent(i)<<"   corrs: "<< CorrValues[p1->GetBinCenter(i)]<<"   subnoise corrs: "<< CorrValuesSubNoise[p1->GetBinCenter(i)]<<endl;
0112
0113
0114 //if(depth>0) {
0115
0116 //diff1-> Fill(iEta, (value-CorrValues[iEta])/value);
0117 //   diff2-> Fill(i, (CorrValues[i] - CorrValuesSubNoise[i])/CorrValues[i]);
0118
0119    //corrs1-> Fill(iEta, value);
0120  }
0121  // }
0122
0123 for (Int_t i=-41; i<=41; i++)
0124 {
0125   corrs2-> Fill(i, CorrValues[i]);
0126   corrs3-> Fill(i, CorrValuesSubNoise[i]);
0127
0128  }
0129
0130
0131
0132 TCanvas* c1 = new TCanvas("c1","all",0,0,350,350);
0133   c1-> cd();
0134   p1 -> Draw("");
0135   // p1 -> SetMaximum(60.);
0136   p1 -> SetXTitle("iEta");
0137   c1->SaveAs(imgpath+"a01.png");
0138
0139
0140   p2 -> Draw("");
0141   p2 -> SetXTitle("iEta");
0142    c1->SaveAs(imgpath+"a02.png");
0143
0144
0145   p3 -> Draw("");
0146   p3 -> SetXTitle("iEta");
0147    c1->SaveAs(imgpath+"a03.png");
0148
0149   p4 -> Draw("");
0150   p4 -> SetXTitle("iEta");
0151    c1->SaveAs(imgpath+"a04.png");
0152
0153
0154   p5 -> Draw("");
0155   p5 -> SetXTitle("iEta");
0156   c1->SaveAs(imgpath+"a05.png");
0157
0158
0159   respcorrs -> Draw("");
0160   respcorrs -> SetMaximum(1.6);
0161   respcorrs -> SetMinimum(0.6);
0162   corrs2 -> Draw("same");
0163   corrs2 -> SetLineColor(kBlue);
0164   c1->SaveAs(imgpath+"c01.png");
0165
0166   respcorrs -> Draw("");
0167   corrs3 -> Draw("same");
0168   respcorrs -> SetMaximum(1.6);
0169   respcorrs -> SetMinimum(0.6);
0170   //corrs2 -> SetLineColor(kBlue);
0171   corrs3 -> SetLineColor(kRed+1);
0172   c1->SaveAs(imgpath+"c02.png");
0173
0174   //diff1 -> Draw("");
0175   //c1->SaveAs(imgpath+"c03.png");
0176
0177   //diff2 -> Draw("");
0178   //c1->SaveAs(imgpath+"c04.png");
0179
0180
0181   tf->Close();
0182
0183
0184
0185   std::string line;
0186
0187   while (getline(old_pfcorrs, line))
0188     {
0189       if(!line.size() || line[0]=='#')
0190     {
0191       //fprintf(new_pfcorrs, "%1s%16s%16s%16s%16s%9s%11s\n", "#", "eta", "phi", "depth", "det", "value", "DetId");
0192       new_pfcorrs<<"#             eta             phi           depth        det           value       DetId"<<endl;
0193       continue;
0194     }
0195
0196       std::istringstream linestream(line);
0197       linestream >> iEta >> iPhi >> depth >> sdName >> value >> hex >> detId;
0198       // corrs1-> Fill(abs(iEta), value);
0199
0200       //cout<<" I'm HERE"<<iEta<<iPhi<<endl;
0201
0202
0203       //format:
0204       //               -1               1               1              HB    0.97635     42004081
0205
0206       //  fprintf(new_pfcorrs, "%17i%16i%16i%16s%9.5f%11X\n", iEta, iPhi, depth, sdName, value, detId);
0207
0208
0209       new_pfcorrs.width(17);
0210       new_pfcorrs<<dec<<iEta;
0211       new_pfcorrs.width(16);
0212       new_pfcorrs<<iPhi;
0213       new_pfcorrs.width(16);
0214       new_pfcorrs<<depth;
0215       new_pfcorrs.width(11);
0216       new_pfcorrs<<sdName;
0217       new_pfcorrs.width(16);
0218       new_pfcorrs.setf(ios::fixed, ios::floatfield);
0219       new_pfcorrs.precision(5);
0220       if (sdName=="HO" || depth==-99) {new_pfcorrs<<value;} else new_pfcorrs<<CorrValues[iEta];
0221       new_pfcorrs.width(13);
0222       new_pfcorrs.setf(ios::uppercase);
0223       new_pfcorrs<<hex<<detId<<endl;
0224
0225
0226       new_pfcorrs_subnoise.width(17);
0227       new_pfcorrs_subnoise<<dec<<iEta;
0228       new_pfcorrs_subnoise.width(16);
0229       new_pfcorrs_subnoise<<iPhi;
0230       new_pfcorrs_subnoise.width(16);
0231       new_pfcorrs_subnoise<<depth;
0232       new_pfcorrs_subnoise.width(11);
0233       new_pfcorrs_subnoise<<sdName;
0234       new_pfcorrs_subnoise.width(16);
0235       new_pfcorrs_subnoise.setf(ios::fixed, ios::floatfield);
0236       new_pfcorrs_subnoise.precision(5);
0237       if (sdName=="HO" || depth==-99) {new_pfcorrs_subnoise<<value;} else new_pfcorrs_subnoise<<CorrValues[iEta];
0238       new_pfcorrs_subnoise.width(13);
0239       new_pfcorrs_subnoise.setf(ios::uppercase);
0240       new_pfcorrs_subnoise<<hex<<detId<<endl;
0241
0242     }
0243 }
0244