File indexing completed on 2023-03-17 10:43:45
0001 #define TreeAnalysisHcalScale_cxx
0002 #include "TreeAnalysisHcalScale.h"
0003 #include <TStyle.h>
0004 #include <TCanvas.h>
0005
0006 Bool_t FillChain(TChain *chain, const char *inputFileList);
0007
0008 int main(Int_t argc, Char_t *argv[]) {
0009 if( argc<4 ){
0010 std::cerr << "Please give 4 arguments \n"
0011 << "runList SeedName" << "\n"
0012 << "outputFileName" << "\n"
0013 << "maximum sample size" << "\n"
0014 << "Which sample to do (-1 means all)" << "\n"
0015 << std::endl;
0016 return -1;
0017 }
0018
0019 const char *inputFileList = argv[1];
0020 const char *outFileName = argv[2];
0021 double totalTracks = atof(argv[3]);
0022 int sample = atoi(argv[4]);
0023 if (sample < 0 || sample > 5) sample = -1;
0024
0025 int cut = 0;
0026
0027 std::cout << "runList SeeName " << inputFileList << "\n"
0028 << "outputFileName " << outFileName << "\n"
0029 << "total number of tracks " << totalTracks
0030 << " Sample " << sample << std::endl;
0031
0032
0033 std::cout << "---------------------" << std::endl;
0034 std::cout << "Reading List of input trees from " << inputFileList << std::endl;
0035 std::string particles[6] = {"pi+","pi-","K+", "K-", "p+", "p-"};
0036 double evFrac[6] = {0.695, 0.695, 0.163, 0.163, 0.142, 0.142};
0037 std::vector<std::string> particleNames;
0038 std::vector<double> fraction;
0039 if (sample < 0) {
0040 for (unsigned int i=0; i<6; i++) {
0041 particleNames.push_back(particles[i]);
0042 fraction.push_back(evFrac[i]);
0043 }
0044 } else {
0045 particleNames.push_back(particles[sample]);
0046 fraction.push_back(1.0);
0047 }
0048
0049 TreeAnalysisHcalScale tree(outFileName, particleNames);
0050
0051 for (unsigned int i=0; i<particleNames.size(); i++) {
0052 char fileList[200], treeName[200];
0053 sprintf (fileList, "%s%s.txt", inputFileList, particles[i].c_str());
0054
0055 sprintf (treeName, "/isolatedTracksHcal/tree");
0056 TChain *chain = new TChain(treeName);
0057 std::cout << "try to create a chain for " << fileList << std::endl;
0058 if (! FillChain(chain, fileList) ) {
0059 std::cerr << "Cannot get the tree " << std::endl;
0060 } else {
0061 unsigned int nmax = (unsigned int)(fraction[i]*totalTracks);
0062 tree.Init(chain);
0063 tree.setParticle(i, nmax);
0064 tree.Loop(cut);
0065 tree.weights[i]= (double)nmax/(double)tree.nIsoTrkTotal;
0066 std::cout << "particle " << i << " nmax " << nmax << " nisotrktotal " << tree.nIsoTrkTotal << " weight " << tree.weights[i] << std::endl;
0067 tree.clear();
0068
0069 }
0070 }
0071 tree.AddWeight(particleNames);
0072
0073 return 0;
0074 }
0075
0076 Bool_t FillChain(TChain *chain, const char *inputFileList) {
0077
0078 ifstream infile(inputFileList);
0079 std::string buffer;
0080
0081 if(!infile.is_open()) {
0082 std::cerr << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
0083 return kFALSE;
0084 }
0085
0086
0087 while(1) {
0088 infile >> buffer;
0089 if(!infile.good()) break;
0090
0091
0092 chain->Add(buffer.c_str());
0093 }
0094 std::cout << "No. of Entries in this tree : " << chain->GetEntries() << std::endl;
0095 return kTRUE;
0096 }
0097
0098 TreeAnalysisHcalScale::TreeAnalysisHcalScale(const char *outFileName, std::vector<std::string>& particles) {
0099
0100 ipBin = nmaxBin = 0;
0101 double tempgen_Eta[NPBins+1] = {20.0, 30.0, 40.0, 60.0, 100.0, 1000.0};
0102
0103 for(int i=0; i<NPBins+1; i++) genPartPBins[i] = tempgen_Eta[i];
0104 BookHistograms(outFileName, particles);
0105 }
0106
0107 TreeAnalysisHcalScale::~TreeAnalysisHcalScale() {
0108 fout->cd();
0109 fout->Write();
0110 fout->Close();
0111 }
0112
0113 Int_t TreeAnalysisHcalScale::Cut(Long64_t entry) {
0114
0115
0116
0117 return 1;
0118 }
0119
0120 Int_t TreeAnalysisHcalScale::GetEntry(Long64_t entry) {
0121
0122
0123 if (!fChain) return 0;
0124 return fChain->GetEntry(entry);
0125 }
0126
0127 Long64_t TreeAnalysisHcalScale::LoadTree(Long64_t entry) {
0128
0129
0130 if (!fChain) return -5;
0131 Long64_t centry = fChain->LoadTree(entry);
0132 if (centry < 0) return centry;
0133 if (!fChain->InheritsFrom(TChain::Class())) return centry;
0134 TChain *chain = (TChain*)fChain;
0135 if (chain->GetTreeNumber() != fCurrent) {
0136 fCurrent = chain->GetTreeNumber();
0137 Notify();
0138 }
0139 return centry;
0140 }
0141
0142 void TreeAnalysisHcalScale::Init(TChain *tree) {
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152 t_trackP = 0;
0153 t_trackPt = 0;
0154 t_trackEta = 0;
0155 t_trackPhi = 0;
0156 t_trackHcalEta = 0;
0157 t_trackHcalPhi = 0;
0158 t_hCone = 0;
0159 t_conehmaxNearP = 0;
0160 t_eMipDR = 0;
0161 t_eMipDR_2 = 0;
0162 t_eECALDR = 0;
0163 t_eECALDR_2 = 0;
0164 t_eHCALDR = 0;
0165 t_e11x11_20Sig = 0;
0166 t_e15x15_20Sig = 0;
0167
0168
0169 if (!tree) return;
0170 fChain = tree;
0171 fCurrent = -1;
0172 fChain->SetMakeClass(1);
0173
0174
0175 fChain->SetBranchAddress("t_RunNo", &t_RunNo, &b_t_RunNo);
0176 fChain->SetBranchAddress("t_Lumi", &t_Lumi, &b_t_Lumi);
0177 fChain->SetBranchAddress("t_Bunch", &t_Bunch, &b_t_Bunch);
0178 fChain->SetBranchAddress("t_trackP", &t_trackP, &b_t_trackP);
0179 fChain->SetBranchAddress("t_trackPt", &t_trackPt, &b_t_trackPt);
0180 fChain->SetBranchAddress("t_trackEta", &t_trackEta, &b_t_trackEta);
0181 fChain->SetBranchAddress("t_trackPhi", &t_trackPhi, &b_t_trackPhi);
0182 fChain->SetBranchAddress("t_trackHcalEta", &t_trackHcalEta, &b_t_trackHcalEta);
0183 fChain->SetBranchAddress("t_trackHcalPhi", &t_trackHcalPhi, &b_t_trackHcalPhi);
0184 fChain->SetBranchAddress("t_hCone", &t_hCone, &b_t_hCone);
0185 fChain->SetBranchAddress("t_conehmaxNearP", &t_conehmaxNearP, &b_t_conehmaxNearP);
0186 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0187 fChain->SetBranchAddress("t_eMipDR_2", &t_eMipDR_2, &b_t_eMipDR_2);
0188 fChain->SetBranchAddress("t_eECALDR", &t_eECALDR, &b_t_eECALDR);
0189 fChain->SetBranchAddress("t_eECALDR_2", &t_eECALDR_2, &b_t_eECALDR_2);
0190 fChain->SetBranchAddress("t_eHCALDR", &t_eHCALDR, &b_t_eHCALDR);
0191 fChain->SetBranchAddress("t_e11x11_20Sig", &t_e11x11_20Sig, &b_t_e11x11_20Sig);
0192 fChain->SetBranchAddress("t_e15x15_20Sig", &t_e15x15_20Sig, &b_t_e15x15_20Sig);
0193 Notify();
0194 }
0195
0196 void TreeAnalysisHcalScale::Loop(int cut) {
0197
0198 if (fChain == 0) return;
0199
0200 Long64_t nentries = fChain->GetEntriesFast();
0201 std::cout << "No. of Entries in tree " << nentries << std::endl;
0202
0203 Long64_t nEventsGoodRuns=0, nEventsValidPV=0, nEventsPVTracks=0;
0204
0205 Long64_t nbytes = 0, nb = 0;
0206 std::map<unsigned int, unsigned int> runEvtList, runNTrkList, runNTrkEtaPList, runNTrkMipList, runNTrkMipCharIsoList, runNIsoTrkList;
0207 nIsoTrkTotal=0;
0208
0209 int nIsoTrkEtaBin[NEtaBins], nGoodTrkEtaBin[NEtaBins];
0210 for(int i=0; i<NEtaBins; i++){
0211 nIsoTrkEtaBin[i]=0;
0212 nGoodTrkEtaBin[i]=0;
0213 }
0214 for (Long64_t jentry=0; jentry<nentries;jentry++) {
0215
0216
0217 Long64_t ientry = LoadTree(jentry);
0218 if (ientry < 0) break;
0219 nb = fChain->GetEntry(jentry); nbytes += nb;
0220
0221 if( !(jentry%100000) ) {
0222 std::cout << "procesing event " << jentry+1 << std::endl;
0223 }
0224
0225 bool goodRun=false;
0226 goodRun = true;
0227 bool evtSel = true;
0228 if( runEvtList.find(t_RunNo) != runEvtList.end() ) {
0229 runEvtList[t_RunNo] += 1;
0230 } else {
0231 runEvtList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,1) );
0232 std::cout << "runNo " << t_RunNo <<" "<<runEvtList[t_RunNo]<<std::endl;
0233 }
0234
0235 nEventsGoodRuns++;
0236
0237 if( runNTrkList.find(t_RunNo) != runNTrkList.end() ) {
0238 runNTrkList[t_RunNo] += t_trackP->size();
0239 } else {
0240 runNTrkList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,t_trackP->size()) );
0241 std::cout << "runNo " << t_RunNo <<" "<<runNTrkList[t_RunNo]<<std::endl;
0242 }
0243
0244 unsigned int NIsoTrk=0, NEtaPTrk=0, NMipTrk=0, NMipCharIsoTrk=0;
0245 for(int itrk=0; itrk<t_trackP->size(); itrk++ ){
0246 int iTrkEtaBin=-1, iTrkMomBin=-1;
0247 if(std::abs((*t_trackHcalEta)[itrk]) < 26)
0248 iTrkEtaBin = ((*t_trackHcalEta)[itrk] + 25);
0249 else continue;
0250
0251 for(int ip=0; ip<NPBins; ip++) {
0252 if( (*t_trackP)[itrk]>genPartPBins[ip] && (*t_trackP)[itrk]<genPartPBins[ip+1] ) iTrkMomBin = ip;
0253 }
0254 if(iTrkMomBin==2) nGoodTrkEtaBin[iTrkEtaBin]++;
0255
0256 if(iTrkEtaBin>=0 && iTrkMomBin>=0) {
0257 NEtaPTrk++;
0258 h_trackP[ipBin] ->Fill((*t_trackP)[itrk] );
0259 h_trackPt[ipBin] ->Fill((*t_trackPt)[itrk] );
0260 h_trackEta[ipBin] ->Fill((*t_trackEta)[itrk] );
0261 h_trackPhi[ipBin] ->Fill((*t_trackPhi)[itrk] );
0262 h_trackHcalEta[ipBin] ->Fill((*t_trackHcalEta)[itrk] );
0263 h_trackHcalPhi[ipBin] ->Fill((*t_trackHcalPhi)[itrk] );
0264
0265 h_hCone[ipBin] ->Fill( (*t_hCone)[itrk] );
0266 h_conehmaxNearP[ipBin] ->Fill( (*t_conehmaxNearP)[itrk] );
0267 h_eMipDR[ipBin] ->Fill( (*t_eMipDR)[itrk] );
0268 h_eECALDR[ipBin] ->Fill( (*t_eECALDR)[itrk] );
0269 h_eHCALDR[ipBin] ->Fill( (*t_eHCALDR)[itrk] );
0270 h_e11x11_20Sig[ipBin] ->Fill( (*t_e11x11_20Sig)[itrk] );
0271 h_e15x15_20Sig[ipBin] ->Fill( (*t_e15x15_20Sig)[itrk] );
0272
0273 bool ChargedIso = true, IfMip = true;
0274 bool eNeutIso = true, hNeutIso = true;
0275
0276 if ((*t_conehmaxNearP)[itrk] >= 0.0) ChargedIso = false;
0277 if ((*t_eMipDR)[itrk] >= 1.0) IfMip = false;
0278
0279 if (std::abs((*t_trackEta)[itrk])<1.47 && ((*t_e15x15_20Sig)[itrk]-(*t_e11x11_20Sig)[itrk]) > 0.5) eNeutIso = false;
0280 else if(((*t_e15x15_20Sig)[itrk]-(*t_e11x11_20Sig)[itrk]) > 2.0) eNeutIso = false;
0281
0282 if (((*t_eHCALDR)[itrk]-(*t_hCone)[itrk]) > 3) hNeutIso = false;
0283
0284 if(IfMip){
0285 NMipTrk++;
0286 if(ChargedIso) NMipCharIsoTrk++;
0287 }
0288
0289 if(eNeutIso && hNeutIso && ChargedIso && IfMip){
0290 h_IsotrackPhi[ipBin] ->Fill((*t_trackPhi)[itrk]);
0291 if(iTrkMomBin==2) nIsoTrkEtaBin[iTrkEtaBin]++;
0292 NIsoTrk++;
0293 h_IsotrackHcalIEta[ipBin]->Fill((*t_trackHcalEta)[itrk]);
0294
0295 h_Response[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill(((*t_hCone)[itrk]+(*t_eMipDR_2)[itrk])/(*t_trackP)[itrk]);
0296 h_Response_trunc[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill(((*t_hCone)[itrk]+(*t_eMipDR_2)[itrk])/(*t_trackP)[itrk]);
0297 h_Response_E11x11[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill(((*t_hCone)[itrk]+(*t_e11x11_20Sig)[itrk])/(*t_trackP)[itrk]);
0298 h_Response_E11x11_trunc[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill(((*t_hCone)[itrk]+(*t_e11x11_20Sig)[itrk])/(*t_trackP)[itrk]);
0299
0300 h_eHcalFrac[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0301 h_eHcalFrac[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0302 h_eHcalFrac_trunc[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0303 h_hneutIso[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill((*t_eHCALDR)[itrk] - (*t_hCone)[itrk]);
0304 h_eneutIso[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill((*t_eECALDR)[itrk] - (*t_eMipDR)[itrk]);
0305 h_eneutIsoNxN[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill((*t_e15x15_20Sig)[itrk]-(*t_e11x11_20Sig)[itrk]);
0306 nIsoTrkTotal++;
0307 if (nIsoTrkTotal <= nmaxBin) {
0308 h_Response[0][iTrkMomBin][iTrkEtaBin] ->Fill(((*t_hCone)[itrk]+(*t_eMipDR_2)[itrk])/(*t_trackP)[itrk]);
0309 h_Response_trunc[0][iTrkMomBin][iTrkEtaBin] ->Fill(((*t_hCone)[itrk]+(*t_eMipDR_2)[itrk])/(*t_trackP)[itrk]);
0310 h_Response_E11x11[0][iTrkMomBin][iTrkEtaBin] ->Fill(((*t_hCone)[itrk]+(*t_e11x11_20Sig)[itrk])/(*t_trackP)[itrk]);
0311 h_Response_E11x11_trunc[0][iTrkMomBin][iTrkEtaBin]->Fill(((*t_hCone)[itrk]+(*t_e11x11_20Sig)[itrk])/(*t_trackP)[itrk]);
0312
0313 h_eHcalFrac[0][iTrkMomBin][iTrkEtaBin] ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0314 h_eHcalFrac_trunc[0][iTrkMomBin][iTrkEtaBin] ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0315 h_hneutIso[0][iTrkMomBin][iTrkEtaBin] ->Fill((*t_eHCALDR)[itrk] - (*t_hCone)[itrk]);
0316 h_eneutIso[0][iTrkMomBin][iTrkEtaBin] ->Fill((*t_eECALDR)[itrk] - (*t_eMipDR)[itrk]);
0317 h_eneutIsoNxN[0][iTrkMomBin][iTrkEtaBin] ->Fill((*t_e15x15_20Sig)[itrk]-(*t_e11x11_20Sig)[itrk]);
0318 }
0319 }
0320 }
0321 }
0322 if( runNIsoTrkList.find(t_RunNo) != runNIsoTrkList.end() ) {
0323 runNIsoTrkList[t_RunNo] += NIsoTrk;
0324 } else {
0325 runNIsoTrkList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,NIsoTrk) );
0326 std::cout << "runNo " << t_RunNo <<" "<<runNIsoTrkList[t_RunNo]<<std::endl;
0327 }
0328 if( runNTrkEtaPList.find(t_RunNo) != runNTrkEtaPList.end() ) {
0329 runNTrkEtaPList[t_RunNo] += NEtaPTrk;
0330 } else {
0331 runNTrkEtaPList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,NEtaPTrk) );
0332 std::cout << "runNo " << t_RunNo <<" "<<runNTrkEtaPList[t_RunNo]<<std::endl;
0333 }
0334
0335 if( runNTrkMipList.find(t_RunNo) != runNTrkMipList.end() ) {
0336 runNTrkMipList[t_RunNo] += NMipTrk;
0337 } else {
0338 runNTrkMipList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,NMipTrk) );
0339 std::cout << "runNo " << t_RunNo <<" "<<runNTrkMipList[t_RunNo]<<std::endl;
0340 }
0341
0342 if( runNTrkMipCharIsoList.find(t_RunNo) != runNTrkMipCharIsoList.end() ) {
0343 runNTrkMipCharIsoList[t_RunNo] += NMipCharIsoTrk;
0344 } else {
0345 runNTrkMipCharIsoList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,NMipCharIsoTrk) );
0346 std::cout << "runNo " << t_RunNo <<" "<<runNTrkMipCharIsoList[t_RunNo]<<std::endl;
0347 }
0348
0349
0350 }
0351
0352 std::cout <<"Number of entries in tree "<< nentries <<" # of isolated tracks "<< nIsoTrkTotal
0353 <<" should be at least " << nmaxBin << std::endl;
0354
0355
0356 std::cout << "saved runEvtList " << std::endl;
0357 std::map<unsigned int, unsigned int>::iterator runEvtListItr = runEvtList.begin();
0358 for(runEvtListItr=runEvtList.begin(); runEvtListItr != runEvtList.end(); runEvtListItr++) {
0359 std::cout<<runEvtListItr->first << " "<< runEvtListItr->second << std::endl;
0360 }
0361
0362 std::cout << "Number of tracks in runs " << std::endl;
0363 std::map<unsigned int, unsigned int>::iterator runNTrkListItr = runNTrkList.begin();
0364 for(runNTrkListItr=runNTrkList.begin(); runNTrkListItr != runNTrkList.end(); runNTrkListItr++) {
0365 std::cout<<runNTrkListItr->first << " "<< runNTrkListItr->second << std::endl;
0366 }
0367
0368 std::cout << "Number of tracks in eta range in runs " << std::endl;
0369 std::map<unsigned int, unsigned int>::iterator runNTrkEtaPListItr = runNTrkEtaPList.begin();
0370 for(runNTrkEtaPListItr=runNTrkEtaPList.begin(); runNTrkEtaPListItr != runNTrkEtaPList.end(); runNTrkEtaPListItr++) {
0371 std::cout<<runNTrkEtaPListItr->first << " "<< runNTrkEtaPListItr->second << std::endl;
0372 }
0373
0374 std::cout << "Number of Mip tracks in runs " << std::endl;
0375 std::map<unsigned int, unsigned int>::iterator runNTrkMipListItr = runNTrkMipList.begin();
0376 for(runNTrkMipListItr=runNTrkMipList.begin(); runNTrkMipListItr != runNTrkMipList.end(); runNTrkMipListItr++) {
0377 std::cout<<runNTrkMipListItr->first << " "<< runNTrkMipListItr->second << std::endl;
0378 }
0379
0380
0381 std::cout << "Number of Charged isolated Mip tracks in runs " << std::endl;
0382 std::map<unsigned int, unsigned int>::iterator runNTrkMipCharIsoListItr = runNTrkMipCharIsoList.begin();
0383 for(runNTrkMipCharIsoListItr=runNTrkMipCharIsoList.begin(); runNTrkMipCharIsoListItr != runNTrkMipCharIsoList.end(); runNTrkMipCharIsoListItr++) {
0384 std::cout<<runNTrkMipCharIsoListItr->first << " "<< runNTrkMipCharIsoListItr->second << std::endl;
0385 }
0386
0387 std::cout << "Number of isolated tracks in runs " << std::endl;
0388 std::map<unsigned int, unsigned int>::iterator runNIsoTrkListItr = runNIsoTrkList.begin();
0389 for(runNIsoTrkListItr=runNIsoTrkList.begin(); runNIsoTrkListItr != runNIsoTrkList.end(); runNIsoTrkListItr++) {
0390 std::cout<<runNIsoTrkListItr->first << " "<< runNIsoTrkListItr->second << std::endl;
0391 }
0392 for(int i=0; i<NEtaBins;i++){
0393 std::cout<< "Number of tracks in ieta " << i-25 << ": " << nIsoTrkEtaBin[i] << std::endl;
0394 if(nGoodTrkEtaBin[i]!=0){
0395 double frac = (double)nIsoTrkEtaBin[i]/(double)nGoodTrkEtaBin[i];
0396 h_FracIsotrackHcalIEta[ipBin]->Fill(i-25, frac);
0397 }
0398 }
0399 }
0400
0401 Bool_t TreeAnalysisHcalScale::Notify() {
0402
0403
0404
0405
0406
0407
0408 return kTRUE;
0409 }
0410
0411 void TreeAnalysisHcalScale::Show(Long64_t entry) {
0412
0413
0414 if (!fChain) return;
0415 fChain->Show(entry);
0416 }
0417
0418 void TreeAnalysisHcalScale::BookHistograms(const char *outFileName, std::vector<std::string>& particles) {
0419
0420 fout = new TFile(outFileName, "RECREATE");
0421
0422 fout->cd();
0423
0424 char hname[1000], htit[1000];
0425
0426 for (unsigned int i=0; i<particles.size(); i++) {
0427 sprintf (hname, "trackEtaAll%s", particles[i].c_str());
0428 sprintf (htit, "Eta: All Tracks for %s", particles[i].c_str());
0429 h_trackEtaAll[i] = new TH1F(hname, htit, 100, -3.0, 3.0 );
0430 h_trackEtaAll[i] ->Sumw2();
0431 sprintf (hname, "trackP%s", particles[i].c_str());
0432 sprintf (htit, "P: Tracks in eta region for %s", particles[i].c_str());
0433 h_trackP[i] = new TH1F(hname, htit, 5000, 0.0, 1000 );
0434 h_trackP[i] ->Sumw2();
0435 sprintf (hname, "trackPt%s", particles[i].c_str());
0436 sprintf (htit, "Pt: Tracks in eta region for %s", particles[i].c_str());
0437 h_trackPt[i] = new TH1F(hname, htit, 5000, 0.0, 1000 );
0438 h_trackPt[i] ->Sumw2();
0439 sprintf (hname, "trackEta%s", particles[i].c_str());
0440 sprintf (htit, "Eta: Tracks in eta region for %s", particles[i].c_str());
0441 h_trackEta[i] = new TH1F(hname, htit, 100, -3.0, 3.0 );
0442 h_trackEta[i] ->Sumw2();
0443 sprintf (hname, "trackPhi%s", particles[i].c_str());
0444 sprintf (htit, "Phi: Tracks in eta region for %s", particles[i].c_str());
0445 h_trackPhi[i] = new TH1F(hname, htit, 100, -4.0, 4.0 );
0446 h_trackPhi[i] ->Sumw2();
0447 sprintf (hname, "IsotrackPhi%s", particles[i].c_str());
0448 sprintf (htit, "Phi: Isolated tracks for %s", particles[i].c_str());
0449 h_IsotrackPhi[i] = new TH1F(hname, htit, 100, -4.0, 4.0 );
0450 h_IsotrackPhi[i] ->Sumw2();
0451 sprintf (hname, "trackHcalEta%s", particles[i].c_str());
0452 sprintf (htit, "iEta: Track Hit at Hcal for %s", particles[i].c_str());
0453 h_trackHcalEta[i] = new TH1F(hname, htit, 100, -50.0, 50.0 );
0454 h_trackHcalEta[i] ->Sumw2();
0455 sprintf (hname, "IsotrackHcalIEta%s", particles[i].c_str());
0456 sprintf (htit, "iEta: Isolated Track Hit at Hcal for %s", particles[i].c_str());
0457 h_IsotrackHcalIEta[i] = new TH1F(hname, htit, 60, -30.0, 30.0 );
0458 h_IsotrackHcalIEta[i] ->Sumw2();
0459 sprintf (hname, "FracIsotrackHcalIEta%s", particles[i].c_str());
0460 sprintf (htit, "iEta: Fraction of Isolated Track Hit at Hcal for %s", particles[i].c_str());
0461 h_FracIsotrackHcalIEta[i] = new TH1F(hname, htit, 60, -30.0, 30.0 );
0462 h_FracIsotrackHcalIEta[i] ->Sumw2();
0463
0464 sprintf (hname, "trackHcalPhi%s", particles[i].c_str());
0465 sprintf (htit, "iPhi: Track hit at Hcal for %s", particles[i].c_str());
0466 h_trackHcalPhi[i] = new TH1F(hname, htit, 100, 0.0, 100. );
0467 h_trackHcalPhi[i] ->Sumw2();
0468 sprintf (hname, "h_hcone%s", particles[i].c_str());
0469 sprintf (htit, "Energy in Hcal in the cone for %s", particles[i].c_str());
0470 h_hCone[i] = new TH1F(hname, htit, 5000, -2.0, 1000.0 );
0471 h_hCone[i] ->Sumw2();
0472 sprintf (hname, "h_conehmaxNearP%s", particles[i].c_str());
0473 sprintf (htit, "Max energy in charge isolation for %s", particles[i].c_str());
0474 h_conehmaxNearP[i] = new TH1F(hname, htit, 5000, -2.0, 1000.0 );
0475 h_conehmaxNearP[i] ->Sumw2();
0476 sprintf (hname, "h_eMipDR%s", particles[i].c_str());
0477 sprintf (htit, "Energy in the MIP region for %s", particles[i].c_str());
0478 h_eMipDR[i] = new TH1F(hname, htit, 5000, -2.0, 1000.0 );
0479 h_eMipDR[i] ->Sumw2();
0480 sprintf (hname, "h_eECALDR%s", particles[i].c_str());
0481 sprintf (htit, "Energy in ECAL iso region for %s", particles[i].c_str());
0482 h_eECALDR[i] = new TH1F(hname, htit, 5000, -2.0, 1000.0 );
0483 h_eECALDR[i] ->Sumw2();
0484 sprintf (hname, "h_eHCALDR%s", particles[i].c_str());
0485 sprintf (htit, "Energy in HCAL iso region for %s", particles[i].c_str());
0486 h_eHCALDR[i] = new TH1F(hname, htit, 5000, -2.0, 1000.0 );
0487 h_eHCALDR[i] ->Sumw2();
0488 sprintf (hname, "h_e11x11_20Sig%s", particles[i].c_str());
0489 sprintf (htit, "Energy in 11x11 for %s", particles[i].c_str());
0490 h_e11x11_20Sig[i] = new TH1F(hname ,htit, 5000, -2.0, 1000.0 );
0491 h_e11x11_20Sig[i] ->Sumw2();
0492 sprintf (hname, "h_e15x15_20Sig%s", particles[i].c_str());
0493 sprintf (htit, "Energy in 15x15 for %s", particles[i].c_str());
0494 h_e15x15_20Sig[i] = new TH1F(hname, htit, 5000, -2.0, 1000.0 );
0495 h_e15x15_20Sig[i] ->Sumw2();
0496 }
0497
0498 TDirectory *d_HcalFrac = fout->mkdir( "HcalFrac" );
0499 TDirectory *d_eNeutIso = fout->mkdir( "eNeutIso" );
0500 TDirectory *d_hNeutIso = fout->mkdir( "hNeutIso" );
0501 TDirectory *d_eNeutIsoNxN = fout->mkdir( "eNeutIsoNxN" );
0502 TDirectory *d_Response = fout->mkdir( "Response" );
0503
0504 for (unsigned int i=0; i<particles.size()+1; i++) {
0505 char part[10];
0506 if (i == 0) sprintf (part, "all");
0507 else sprintf (part, "%s", particles[i-1].c_str());
0508 for (int ieta=0; ieta<NEtaBins; ieta++) {
0509 int iEta = ieta;
0510 for (int ip=0; ip<NPBins; ip++) {
0511 double lowMom=-5.0, highMom= 5.0;
0512 lowMom = genPartPBins[ip];
0513 highMom = genPartPBins[ip+1];
0514
0515 d_Response->cd();
0516 sprintf(hname, "h_Response_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0517 sprintf(htit, "Response (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0518 h_Response[i][ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0519 h_Response[i][ip][ieta] ->Sumw2();
0520
0521 sprintf(hname, "h_Response_trunc_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0522 sprintf(htit, "Response_trunc (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0523 h_Response_trunc[i][ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0524 h_Response_trunc[i][ip][ieta] ->Sumw2();
0525
0526 sprintf(hname, "h_Response_E11x11_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0527 sprintf(htit, "Response_E11x11 (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0528 h_Response_E11x11[i][ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0529 h_Response_E11x11[i][ip][ieta] ->Sumw2();
0530
0531 sprintf(hname, "h_Response_E11x11_trunc_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0532 sprintf(htit, "Response_E11x11_trunc (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0533 h_Response_E11x11_trunc[i][ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0534 h_Response_E11x11_trunc[i][ip][ieta] ->Sumw2();
0535
0536 d_HcalFrac->cd();
0537 sprintf(hname, "h_eHcalFrac_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0538 sprintf(htit, "eHcalFrac (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0539 h_eHcalFrac[i][ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0540 h_eHcalFrac[i][ip][ieta] ->Sumw2();
0541
0542 sprintf(hname, "h_eHcalFrac_trunc_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0543 sprintf(htit, "eHcalFrac_trunc (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0544 h_eHcalFrac_trunc[i][ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0545 h_eHcalFrac_trunc[i][ip][ieta] ->Sumw2();
0546
0547 d_hNeutIso->cd();
0548 sprintf(hname, "h_hneutIso_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0549 sprintf(htit, "hneutIso (|i#eta|=%i),(%3.2f<|#p|<%3.2f for %s)", iEta, lowMom, highMom, part);
0550 h_hneutIso[i][ip][ieta] = new TH1F(hname, htit, 1000, -5.0, 20.0);
0551 h_hneutIso[i][ip][ieta] ->Sumw2();
0552
0553 d_eNeutIso->cd();
0554 sprintf(hname, "h_eneutIso_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0555 sprintf(htit, "eneutIso (|i#eta|=%i),(%3.2f<|#p|<%3.2f for %s)", iEta, lowMom, highMom, part);
0556 h_eneutIso[i][ip][ieta] = new TH1F(hname, htit, 1000, -1.0, 25.0);
0557 h_eneutIso[i][ip][ieta] ->Sumw2();
0558
0559 d_eNeutIsoNxN->cd();
0560 sprintf(hname, "h_eneutIsoNxN_etaBin%i_pBin%i_particle%s",ieta,ip,part);
0561 sprintf(htit, "eneutIsoNxN (|i#eta|=%i),(%3.2f<|#p|<%3.2f for %s)", iEta, lowMom, highMom, part);
0562 h_eneutIsoNxN[i][ip][ieta] = new TH1F(hname, htit, 1000, -1.0, 10.0);
0563 h_eneutIsoNxN[i][ip][ieta] ->Sumw2();
0564 }
0565 }
0566 }
0567 for (int ieta=0; ieta<NEtaBins; ieta++) {
0568 int iEta = ieta;
0569 for (int ip=0; ip<NPBins; ip++) {
0570 double lowMom=-5.0, highMom= 5.0;
0571 lowMom = genPartPBins[ip];
0572 highMom = genPartPBins[ip+1];
0573
0574 d_HcalFrac->cd();
0575 sprintf(hname, "h_eHcalFrac_all_etaBin%i_pBin%i", ieta,ip);
0576 sprintf(htit, "eHcalFrac allWeighted (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0577 h_eHcalFrac_all[ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0578 h_eHcalFrac_all[ip][ieta] ->Sumw2();
0579
0580 sprintf(hname, "h_eHcalFrac_trunc_all_etaBin%i_pBin%i", ieta,ip);
0581 sprintf(htit, "eHcalFrac_trunc allWeigthed (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0582 h_eHcalFrac_trunc_all[ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0583 h_eHcalFrac_trunc_all[ip][ieta] ->Sumw2();
0584
0585 d_Response->cd();
0586 sprintf(hname, "h_Response_all_etaBin%i_pBin%i", ieta,ip);
0587 sprintf(htit, "Response allWeighted (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0588 h_Response_all[ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0589 h_Response_all[ip][ieta] ->Sumw2();
0590
0591 sprintf(hname, "h_Response_trunc_all_etaBin%i_pBin%i", ieta,ip);
0592 sprintf(htit, "Response_trunc allWeigthed (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0593 h_Response_trunc_all[ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0594 h_Response_trunc_all[ip][ieta] ->Sumw2();
0595
0596 sprintf(hname, "h_Response_E11x11_all_etaBin%i_pBin%i", ieta,ip);
0597 sprintf(htit, "Response_E11x11 allWeighted (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0598 h_Response_E11x11_all[ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0599 h_Response_E11x11_all[ip][ieta] ->Sumw2();
0600
0601 sprintf(hname, "h_Response_E11x11_trunc_all_etaBin%i_pBin%i", ieta,ip);
0602 sprintf(htit, "Response_E11x11_trunc allWeigthed (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0603 h_Response_E11x11_trunc_all[ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0604 h_Response_E11x11_trunc_all[ip][ieta] ->Sumw2();
0605
0606 }
0607 }
0608 fout->cd();
0609
0610 }
0611
0612 void TreeAnalysisHcalScale::clear() {
0613 if (!fChain) return;
0614 delete fChain->GetCurrentFile();
0615 }
0616
0617 void TreeAnalysisHcalScale::setParticle(unsigned int ip, unsigned int nmax) {
0618 ipBin = ip;
0619 nmaxBin = nmax;
0620 }
0621
0622 void TreeAnalysisHcalScale::AddWeight( std::vector<std::string> particleNames){
0623
0624 for (unsigned int i=0; i<particleNames.size(); i++) {
0625 char fileList[200], treeName[200];
0626 std::cout << "particle " << i << " weight " << weights[i] << std::endl;
0627 for (int ieta=0; ieta<NEtaBins; ieta++) {
0628 for (int ip=0; ip<NPBins; ip++) {
0629 h_eHcalFrac_all[ip][ieta]->Add(h_eHcalFrac[i+1][ip][ieta], weights[i]);
0630 h_eHcalFrac_trunc_all[ip][ieta]->Add(h_eHcalFrac_trunc[i+1][ip][ieta], weights[i]);
0631
0632 h_Response_all[ip][ieta]->Add(h_Response[i+1][ip][ieta], weights[i]);
0633 h_Response_trunc_all[ip][ieta]->Add(h_Response_trunc[i+1][ip][ieta], weights[i]);
0634
0635 h_Response_E11x11_all[ip][ieta]->Add(h_Response_E11x11[i+1][ip][ieta], weights[i]);
0636 h_Response_E11x11_trunc_all[ip][ieta]->Add(h_Response_E11x11_trunc[i+1][ip][ieta], weights[i]);
0637
0638 }
0639 }
0640 }
0641 }
0642