Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-31 23:17:17

0001 // derive MET pileup mitigation
0002 
0003 // 1) produce nVtx vs nTT4 plot to find pu bins
0004 //    root -l -b -q 'deriveMETPUM.cxx(1,0,0,1)' # doTow, doLUT, doFit, doMC
0005 // 2) fit nVtx vs nTT4 plot
0006 //    root -l -b -q 'deriveMETPUM.cxx(0,0,1,1)' # doTow, doLUT, doFit, doMC
0007 // 3) check fitting of plot and calculate pu binning, copy into "vector<int> puBinBs   = "
0008 //    and run the filling of the pileup binned tower energy histogram arrays
0009 //    root -l -b -q 'deriveMETPUM.cxx(0,0,0,1)' # doTow, doLUT, doFit, doMC 
0010 // 4) produce LUT containing eta and pileup binned MET tower energy thresholds
0011 //    root -l -b -q 'deriveMETPUM.cxx(0,1,0,1)' # doTow, doLUT, doFit, doMC 
0012 
0013 // nb for MC only RAWEMU is needed. for data, AOD reco vertex information is needed
0014 
0015 
0016 #include "TMath.h"
0017 #include "TFile.h"
0018 #include "TTree.h"
0019 #include "TH1F.h"
0020 #include "TChain.h"
0021 #include <TCanvas.h>
0022 #include <iostream>
0023 #include <iomanip>
0024 #include <fstream>
0025 #include <string>
0026 #include "L1Trigger/L1TNtuples/interface/L1AnalysisEventDataFormat.h"
0027 #include "L1Trigger/L1TNtuples/interface/L1AnalysisL1UpgradeDataFormat.h"
0028 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoVertexDataFormat.h"
0029 #include "L1Trigger/L1TNtuples/interface/L1AnalysisCaloTPDataFormat.h"
0030 #include "L1Trigger/L1TNtuples/interface/L1AnalysisL1CaloTowerDataFormat.h"
0031 
0032 
0033 void fitProfile(TH1D* prof){
0034  
0035   TCanvas* profCanvas = new TCanvas("profCanvas","",600,600);
0036  
0037   TH1D *prof_fit2 = (TH1D*)prof->Clone("prof_fit2");
0038   TH1D *prof_fit3 = (TH1D*)prof->Clone("prof_fit3");
0039   TH1D *prof_fit4 = (TH1D*)prof->Clone("prof_fit4");
0040   prof->SetMarkerStyle(7);
0041   prof_fit2->SetMarkerStyle(7);
0042   prof_fit3->SetMarkerStyle(7);
0043   prof_fit4->SetMarkerStyle(7);
0044   
0045   TF1 *fit1 = new TF1("fit1","[0]*x+[1]");
0046   fit1->SetParameter(0,1.0); 
0047   fit1->SetParameter(1,20);   
0048   fit1->SetRange(20,120);     
0049   prof->Fit("fit1","R");
0050   TF1 *fit2 = new TF1("fit2","[0]*x+[1]");
0051   fit2->SetParameter(0,0.8); 
0052   fit2->SetParameter(1,60);   
0053   fit2->SetRange(90,155);
0054   prof_fit2->Fit("fit2","R");
0055   TF1 *fit3 = new TF1("fit3","[0]*x+[1]");
0056   fit3->SetParameter(0,0.5);
0057   fit3->SetParameter(1,4);  
0058   fit3->SetRange(110,140);
0059   prof_fit3->Fit("fit3","R");   
0060   TF1 *fit4 = new TF1("fit4","[0]*x+[1]");
0061   fit4->SetParameter(0,0.5);
0062   fit4->SetParameter(1,4);  
0063   fit4->SetRange(140,160);   
0064   prof_fit4->Fit("fit4","R");
0065   prof->Draw();  
0066   //prof_fit2->Draw("sames");  
0067   //prof_fit3->Draw("sames");  
0068   //prof_fit4->Draw("sames");  
0069   gPad->Update();
0070   TPaveStats *st1 = (TPaveStats*)prof->FindObject("stats");
0071   TPaveStats *st2 = (TPaveStats*)prof_fit2->FindObject("stats");
0072   //TPaveStats *st3 = (TPaveStats*)prof_fit3->FindObject("stats");
0073   //TPaveStats *st4 = (TPaveStats*)prof_fit4->FindObject("stats");
0074   st1->SetY1NDC(0.75); st1->SetY2NDC(0.95); 
0075   //st2->SetY1NDC(0.55); st2->SetY2NDC(0.75);
0076   //st3->SetY1NDC(0.35); st3->SetY2NDC(0.55);  
0077   //st4->SetY1NDC(0.15); st4->SetY2NDC(0.35);  
0078   fit1->Draw("sames"); //fit2->Draw("sames"); //fit3->Draw("sames"); fit4->Draw("sames");
0079   prof->Write();
0080   profCanvas->SaveAs("ProfNVtxNTowemu4.pdf");
0081 
0082 }
0083 
0084 void makeLUT(TFile* file, int puBins){
0085 
0086   vector<vector<TH1D*> > hTowEtPU(puBins, vector<TH1D*> (41));
0087 
0088   for(uint pu=0;pu<puBins;pu++){
0089     for(uint eta=0; eta<41; ++eta){
0090       stringstream tempPUEt;
0091       tempPUEt  << "hTowEt_"     << pu << "_" << eta+1;
0092       hTowEtPU[pu][eta] = (TH1D*)file->Get(tempPUEt.str().c_str());
0093     }
0094   }  
0095 
0096   ofstream one;
0097   ofstream p3;
0098   ofstream p5;
0099   ofstream p05;
0100 
0101   one.open("one.txt");
0102   p3.open("p3.txt");
0103   p5.open("p5.txt");
0104   p05.open("p05.txt");
0105 
0106       
0107   stringstream intro;
0108   intro << "# address to et sum tower Et threshold LUT\n" \
0109     "# maps 11 bits to 9 bits\n" \
0110     "# 11 bits = (compressedPileupEstimate << 6) | abs(ieta)\n" \
0111     "# compressedPileupEstimate is unsigned 5 bits, abs(ieta) is unsigned 6 bits\n" \
0112     "# data: tower energy threshold returned has 9 bits \n" \
0113     "# anything after # is ignored with the exception of the header\n" \
0114     "# the header is first valid line starting with #<header> versionStr nrBitsAddress nrBitsData </header>\n" \
0115     "#<header> v1 11 9 </header>\n"; 
0116 
0117   one.write(intro.str().c_str(), intro.str().length() );
0118   p3.write(intro.str().c_str(), intro.str().length() ); 
0119   p5.write(intro.str().c_str(), intro.str().length() );
0120   p05.write(intro.str().c_str(), intro.str().length() );
0121 
0122   
0123   int addr = 0;
0124 
0125   int lastFilled = 0;
0126 
0127   for(uint pu=0;pu<puBins;pu++){
0128     for(int eta=-1; eta<64; ++eta){ 
0129     
0130       if(eta==28) continue;
0131 
0132       if(eta==-1){
0133           one << addr  << " 0             # nTT4 = " << pu*5 << "-" << pu*5+5 << " ieta = 0\n";
0134           p3 << addr   << " 0             # nTT4 = " << pu*5 << "-" << pu*5+5 << " ieta = 0\n";
0135           p5 << addr   << " 0             # nTT4 = " << pu*5 << "-" << pu*5+5 << " ieta = 0\n";
0136           p05 << addr  << " 0             # nTT4 = " << pu*5 << "-" << pu*5+5 << " ieta = 0\n";
0137 
0138           ++addr;
0139           continue;
0140       }
0141 
0142       if(pu < 4 && eta<41){
0143           one << addr  << " " << 0  << "             # ieta = " << eta+1 << "\n";
0144           p3 << addr   << " " << 0  << "             # ieta = " << eta+1 << "\n";
0145           p5 << addr   << " " << 0  << "             # ieta = " << eta+1 << "\n";
0146           p05 << addr   << " " << 0  << "             # ieta = " << eta+1 << "\n";
0147           ++addr;  
0148           continue;
0149       }
0150     
0151       if(eta>40){
0152           one << addr  << " 0 #dummy\n";
0153           p3 << addr   << " 0 #dummy\n";
0154           p5 << addr   << " 0 #dummy\n";
0155           p05 << addr   << " 0 #dummy\n";
0156           ++addr;  
0157           continue;
0158       }
0159     
0160       double pass = 0;
0161       double d1(999.),dp3(999.),dp5(999.),dp05(999.);
0162       double t1(999.),tp3(999.),tp5(999.),tp05(999.);
0163 
0164       cout << "pu =  " << pu << ", eta = " << eta << endl;
0165 
0166       int nBins = hTowEtPU[pu][eta]->GetNbinsX();
0167       
0168       for(uint t=0; t<nBins; t++){
0169         int puHist = pu;
0170         if(hTowEtPU[pu][eta]->Integral(1,nBins)==0){
0171           if(lastFilled == 0){
0172             t1=tp3=tp5=tp05=0;
0173             break;
0174           } else puHist = lastFilled;
0175         } else {
0176             lastFilled = pu;
0177           }
0178 
0179         pass = hTowEtPU[puHist][eta]->Integral(t+1,nBins)/hTowEtPU[puHist][eta]->Integral(1,nBins);
0180         if( abs(pass-0.01) < d1  ){
0181           t1  = t;
0182           d1 = pass - 0.01;
0183         }
0184         if( abs(pass-0.003) < dp3  ){
0185           tp3  = t;
0186           dp3 = pass - 0.003;
0187         }
0188         if( abs(pass-0.005) < dp5  ){
0189           tp5  = t;
0190           dp5 = pass - 0.005;
0191         }
0192         if( abs(pass-0.0005) < dp05  ){
0193           tp05  = t;
0194           dp05 = pass - 0.0005;
0195         }
0196       }
0197 
0198       //int sat = 32;
0199 
0200       // double rai = 1.2;
0201       // double div = 40;
0202       // double off = 0.2;
0203 
0204       // t1   = round( t1  *((pow(float(pu),rai)/div)+off) );
0205       // tp3   = round( tp3  *((pow(float(pu),rai)/div)+off) );
0206       // tp5  = round( tp5 *((pow(float(pu),rai)/div)+off) );
0207 
0208       //if(eta==27) cout << pu << "   " << ((pow(float(pu),rai)/div)+off) << endl; 
0209 
0210       // if(eta<15){
0211       //    t1 = round(t1*(((double)eta)/14));
0212       //    tp3 = round(tp3*(((double)eta)/14));
0213       //    tp5 = round(tp5*(((double)eta)/14));
0214       // }
0215       
0216       //if(t1>sat)  t1 =sat;
0217       //if(tp3>sat)  tp3 =sat;
0218       //if(tp5>sat)   tp5=sat;
0219 
0220       if(eta<30){
0221         one  << addr  << " "  << t1   << "             # ieta = " << eta+1 << "\n";
0222         p3   << addr  << " "  << tp3  << "             # ieta = " << eta+1 << "\n";
0223         p5   << addr  << " "  << tp5  << "             # ieta = " << eta+1 << "\n"; 
0224         p05  << addr  << " "  << tp05 << "             # ieta = " << eta+1 << "\n";
0225       } else {
0226     int ft1=t1*2, ftp3=tp3*2, ftp5=tp5*2, ftp05=tp05*2;
0227         //int ft1=floor(t1*1.5), ftp3=floor(tp3*1.5), ftp5=floor(tp5*1.5), ftp05=floor(tp05*1.5);
0228         one  << addr  << " "  << ft1   << "             # ieta = " << eta+1 << "\n";
0229         p3   << addr  << " "  << ftp3  << "             # ieta = " << eta+1 << "\n";
0230         p5   << addr  << " "  << ftp5  << "             # ieta = " << eta+1 << "\n"; 
0231         p05  << addr  << " "  << ftp05 << "             # ieta = " << eta+1 << "\n";
0232       }
0233 
0234       ++addr;
0235 
0236     }   
0237   }
0238 
0239   one.close();
0240   p3.close();
0241   p5.close();
0242   p05.close();
0243 
0244 
0245 }
0246 
0247 
0248 // 1d formatter
0249 void formatPlot1D(TH1D* plot1d, int colour){
0250 
0251   plot1d->GetXaxis()->SetTitleOffset(1.2);
0252   plot1d->GetYaxis()->SetTitleOffset(1.2);
0253   plot1d->SetMarkerStyle(7);
0254   plot1d->SetMarkerColor(colour);
0255   plot1d->SetLineColor(colour);
0256   plot1d->SetLineWidth(2);
0257   //plot1d->SetStats(false);
0258 
0259 }
0260 
0261 //2d formatter
0262 void formatPlot2D(TH2D* plot2d){
0263 
0264   plot2d->GetXaxis()->SetTitleOffset(1.4);
0265   plot2d->GetYaxis()->SetTitleOffset(1.4);
0266   //plot2d->SetStats(false);
0267   
0268 }
0269 
0270 
0271 //main plotting function
0272 void deriveMETPUM(bool doTow, bool doLUT, bool doFit, bool doMC){
0273 
0274   gStyle->SetStatW(0.1);
0275   gStyle->SetOptStat("ne");
0276   gStyle->SetOptFit(0001);
0277 
0278   //output filename
0279   string outFilename = "metPUM_out.root";
0280   //vector<int> puBinBs   = {0,0,0,0,0, 0, 0, 0, 1, 5,10,14,18,23,27,32,36,41,45,50, 56, 62, 68, 74, 80, 86, 93, 99,105,111,117,123,999};
0281   vector<int> puBinBs   = {0,0,0,0,4,10,16,22,29,35,41,47,53,59,66,72,78,84,90,96,103,109,115,121,127,133,140,146,152,158,164,170,999};
0282 
0283   int puBins = puBinBs.size()-1;
0284 
0285   TFile* file = TFile::Open( outFilename.c_str() );
0286 
0287   //check file exists
0288   if(doLUT || doFit){
0289     if (file==0){
0290       cout << "TERMINATE: input file does not exist: " << outFilename << endl;
0291       return;
0292     }
0293     if(doFit){
0294       TH1D* prof = (TH1D*)file->Get("hProfNVtxNTowemu4");
0295       fitProfile(prof);
0296       return;
0297     }
0298     if(doLUT){
0299       makeLUT(file, puBins);
0300       return;
0301     }  
0302   }
0303 
0304   //input ntuple
0305   string  inputFile01 = "root://eoscms.cern.ch//eos/cms/store/group/dpg_trigger/comm_trigger/L1Trigger/bundocka/condor/nu_v9NewSF_1637715475/*.root";
0306   
0307   //check file doesn't exist
0308   if (file!=0){
0309     cout << "TERMINATE: not going to overwrite file " << outFilename << endl;
0310     return;
0311   }
0312 
0313   cout << "Loading up the TChain..." << endl;
0314   TChain * eventTree = new TChain("l1EventTree/L1EventTree");
0315   eventTree->Add(inputFile01.c_str());
0316   TChain * treeL1Towemu = new TChain("l1CaloTowerEmuTree/L1CaloTowerTree");
0317   treeL1Towemu->Add(inputFile01.c_str());
0318   TChain * treeTPhw = new TChain("l1CaloTowerTree/L1CaloTowerTree");
0319   treeTPhw->Add(inputFile01.c_str());
0320   TChain * vtxTree;
0321   if(!doMC)
0322   {
0323       vtxTree = new TChain("l1RecoTree/RecoTree");
0324       vtxTree->Add(inputFile01.c_str());
0325   }
0326 
0327   L1Analysis::L1AnalysisEventDataFormat           *event_ = new L1Analysis::L1AnalysisEventDataFormat();
0328   eventTree->SetBranchAddress("Event", &event_);
0329   L1Analysis::L1AnalysisL1CaloTowerDataFormat   *l1Towemu_ = new L1Analysis::L1AnalysisL1CaloTowerDataFormat();
0330   treeL1Towemu->SetBranchAddress("L1CaloTower", &l1Towemu_);
0331   L1Analysis::L1AnalysisCaloTPDataFormat   *l1TPhw_ = new L1Analysis::L1AnalysisCaloTPDataFormat();
0332   treeTPhw->SetBranchAddress("CaloTP", &l1TPhw_);
0333   L1Analysis::L1AnalysisRecoVertexDataFormat        *vtx_ = new L1Analysis::L1AnalysisRecoVertexDataFormat();
0334   if(!doMC)
0335   {
0336       vtxTree->SetBranchAddress("Vertex", &vtx_);
0337   }
0338 
0339 
0340   //initialise histograms
0341   TH1D* hNTow_emu = new TH1D("nTower_emu", ";nTowers; # Events",  150, -0.5, 1499.5);  
0342   TH1D* hNTow4_emu = new TH1D("nTower_emu4", ";nTowers; # Events",  100, -0.5, 499.5);  
0343   TH1D* hNVtx = new TH1D("nVtx_reco", ";nVtx; # Events", 120, -0.5, 119.5);
0344   
0345   TH2D* hNTowemuVsNVtx = new TH2D("nTowemu_vs_nVtx", ";nRecoVtx; nL1Tow", 100, -0.5, 99.5, 125, -0.5, 1499.5);
0346   TH2D* hNTowemuVsNVtx4 = new TH2D("nTowemu_vs_nVtx4", ";nRecoVtx; nL1Tow |i#eta|<5", 70, -0.5, 69.5, 50, -0.5, 199.5);  
0347   
0348   TProfile* hProfNVtxNTowemu = new TProfile("hProfNVtxNTowemu",";nVtx; <nTT>", 60,-0.5,119.5);
0349   TProfile* hProfNVtxNTowemu4 = new TProfile("hProfNVtxNTowemu4",";nVtx; <nTT> abs(i#eta)<5", 128,-0.5,127.5);
0350   TProfile* hProfNTowemuNVtx4 = new TProfile("hProfNTowemuNVtx4",";<nTT> abs(i#eta)<5;nVtx", 40,0,200);
0351 
0352   TH1D* hAllTowEtemu = new TH1D("towerEtEmu", ";Tower E_{T}; # Towers", 40, -0.5, 39.5);
0353   TH1D* hAllTowEtaemu = new TH1D("towerEtaEmu", ";Tower E_{T}; # Towers", 100, -50.5, 49.5);
0354 
0355   TProfile* hProfTowEtEta = new TProfile("hProfTowEtEta",";Tower i#eta;<Tower E_{T}> (GeV)", 100, -50.5, 49.5);
0356   vector<TProfile*>  hProfTowEtEtaPU(puBins);
0357 
0358   TProfile* hProfHCALTPEtEta = new TProfile("hProfHCALTPEtEta",";TP i#eta;<TP E_{T}> (GeV)", 100, -50.5, 49.5);
0359   vector<TProfile*>  hProfHCALTPEtEtaPU(puBins);
0360 
0361   TProfile* hProfECALTPEtEta = new TProfile("hProfECALTPEtEta",";TP i#eta;<TP E_{T}> (GeV)", 100, -50.5, 49.5);
0362   vector<TProfile*>  hProfECALTPEtEtaPU(puBins);
0363 
0364   TH1D* hTowEt = new TH1D("hTowEt",";Tower E_{T}; #Towers", 512, 0, 256);
0365   vector<vector<TH1D*> > hTowEtPU(puBins, vector<TH1D*> (41));
0366 
0367   TH1D* hECALTPEt = new TH1D("hECALTPEt",";ECAL TP E_{T}; #TPs", 256, 0, 128);
0368   vector<vector<TH1D*> > hECALTPEtPU(puBins, vector<TH1D*> (28));
0369 
0370   TH1D* hHCALTPEt = new TH1D("hHCALTPEt",";HCAL TP E_{T}; #TPs", 256, 0, 128);
0371   vector<vector<TH1D*> > hHCALTPEtPU(puBins, vector<TH1D*> (41));
0372 
0373   TH2D* hTowEtaPhi = new TH2D("hTowEtaPhi","; eta; phi", 100, -50, 50, 72, 0, 72);
0374   
0375   for(uint pu=0; pu<puBins;++pu){
0376     stringstream tempPU;
0377     stringstream tempPUE;
0378     stringstream tempPUH;
0379     stringstream tempPUEt;
0380     stringstream tempPUEEt;
0381     stringstream tempPUHEt;
0382     
0383     tempPU << "hProfTowEtEta_" << pu;
0384     tempPUE << "hProfECALTPEtEta_" << pu;
0385     tempPUH << "hProfHCALTPEtEta_" << pu;
0386 
0387     hProfTowEtEtaPU[pu] = new TProfile(tempPU.str().c_str(),";Tower i#eta;<Tower E_{T}> (GeV)", 100, -50.5, 49.5);
0388     hProfECALTPEtEtaPU[pu] = new TProfile(tempPUE.str().c_str(),";ECAL TP i#eta;<TP E_{T}> (GeV)", 100, -50.5, 49.5);
0389     hProfHCALTPEtEtaPU[pu] = new TProfile(tempPUH.str().c_str(),";HCAL TP i#eta;<TP E_{T}> (GeV)", 100, -50.5, 49.5);
0390     
0391     for(uint eta=0; eta<41; eta++){ 
0392       stringstream tempPUEt;
0393       stringstream tempPUEEt;
0394       stringstream tempPUHEt;
0395       tempPUEt  << "hTowEt_"     << pu << "_" << eta+1;
0396       tempPUEEt << "hECALTPEt_" << pu << "_" << eta+1;
0397       tempPUHEt << "hHCALTPEt_" << pu << "_" << eta+1;
0398 
0399       hTowEtPU[pu][eta] = new TH1D(tempPUEt.str().c_str(),";Tower E_{T} (GeV); #Towers", 512, 0, 256);
0400       if(eta<28) hECALTPEtPU[pu][eta] = new TH1D(tempPUEEt.str().c_str(),";ECAL TP E_{T} (GeV); #TPs", 256, 0, 128);
0401       hHCALTPEtPU[pu][eta] = new TH1D(tempPUHEt.str().c_str(),";HCAL TP E_{T} (GeV); #TPs", 256, 0, 128);
0402     }
0403   }
0404       
0405   
0406   //initialise some variables
0407   int nTowemu(-1), nVtx(-1), nTowemu4(-1);
0408   int nVtxBin(-1);
0409   int nHCALTP(-1), nECALTP(-1);
0410   float towEt(-1);
0411   vector<int> towOcc(41);
0412   
0413 
0414   //main loop
0415 
0416   // get number of entries
0417   Long64_t nentries;
0418   nentries = treeL1Towemu->GetEntries();
0419 
0420   // Number of events to run over
0421   int nEvents = nentries; // lol
0422 
0423   for (Long64_t jentry=0; jentry<nEvents; jentry++){
0424   
0425     nTowemu = -1;
0426     nVtx = -1; 
0427     nTowemu4 = -1;
0428     nVtxBin = -1;
0429     nHCALTP = -1;
0430     nECALTP = -1;
0431     towEt = -1;
0432     fill(towOcc.begin(), towOcc.end(), 0);
0433  
0434     //counter
0435     if((jentry%10000)==0) cout << "Done " << jentry  << " events of " << nEvents << endl;
0436     if((jentry%1000)==0) cout << "." << flush;
0437 
0438     eventTree->GetEntry(jentry);
0439     if(!doMC)
0440     {
0441       vtxTree->GetEntry(jentry);
0442       nVtx = vtx_->nVtx;
0443     } else {
0444       eventTree->GetEntry(jentry);
0445       nVtx = event_->nPV;
0446     }
0447 
0448     hNVtx->Fill(nVtx);
0449 
0450     if(nVtx==0) continue;
0451 
0452     treeL1Towemu->GetEntry(jentry);
0453     nTowemu = l1Towemu_->nTower;
0454 
0455     if(!doTow){
0456       //treeTPhw->GetEntry(jentry); 
0457       //nHCALTP = l1TPhw_->nHCALTP;
0458       //nECALTP = l1TPhw_->nECALTP;    
0459       for(uint puIt=0; puIt<puBins;puIt++){
0460           if(nVtx >= puBinBs[puIt] && nVtx < puBinBs[puIt+1]){
0461           nVtxBin = puIt;
0462           break;
0463             }
0464       }
0465     }
0466 
0467     
0468     for(uint towIt=0; towIt<nTowemu; ++towIt){
0469 
0470       hTowEtaPhi->Fill(l1Towemu_->ieta[towIt], l1Towemu_->iphi[towIt]);
0471       hAllTowEtaemu->Fill(l1Towemu_->ieta[towIt]);
0472       hAllTowEtemu->Fill(l1Towemu_->iet[towIt]);
0473       if(abs(l1Towemu_->ieta[towIt])<5) nTowemu4++;
0474       if(!doTow){ 
0475           towOcc[abs(l1Towemu_->ieta[towIt])-1] +=1;
0476           towEt = l1Towemu_->iet[towIt];
0477           hTowEt->Fill(towEt/2);
0478           if(towEt/2 > 120) cout << "Tow Et = " << towEt << endl;
0479               hTowEtPU[nVtxBin][abs(l1Towemu_->ieta[towIt])-1]->Fill(towEt/2);
0480               hProfTowEtEta->Fill(l1Towemu_->ieta[towIt],l1Towemu_->iet[towIt]);
0481               hProfTowEtEtaPU[nVtxBin]->Fill(l1Towemu_->ieta[towIt],l1Towemu_->iet[towIt]);
0482           }
0483       }
0484 
0485     if(!doTow){ 
0486       for(uint eta=0;eta<41;eta++){
0487           hTowEtPU[nVtxBin][eta]->Fill(0.,(144-towOcc[eta]));
0488       }
0489     
0490       
0491       // //fill ECAL TP histos
0492       // for(uint tpIt=0; tpIt<nECALTP; ++tpIt){
0493       //   hECALTPEt->Fill(l1TPhw_->ecalTPet[tpIt]);
0494       //   hECALTPEtPU[nVtxBin][abs(l1TPhw_->ecalTPieta[tpIt])-1]->Fill(l1TPhw_->ecalTPet[tpIt]);
0495       //   hProfECALTPEtEta->Fill(l1TPhw_->ecalTPieta[tpIt],l1TPhw_->ecalTPet[tpIt]);
0496       //   hProfECALTPEtEtaPU[nVtxBin]->Fill(l1TPhw_->ecalTPieta[tpIt],l1TPhw_->ecalTPet[tpIt]);
0497       // }
0498       // //fill HCAL TP histos
0499       // for(uint tpIt=0; tpIt<nHCALTP; ++tpIt){
0500       //   hHCALTPEt->Fill(l1TPhw_->hcalTPet[tpIt]);
0501       //   hHCALTPEtPU[nVtxBin][abs(l1TPhw_->hcalTPieta[tpIt])-1]->Fill(l1TPhw_->hcalTPet[tpIt]);
0502       //   hProfHCALTPEtEta->Fill(l1TPhw_->hcalTPieta[tpIt],l1TPhw_->hcalTPet[tpIt]);
0503       //   hProfHCALTPEtEtaPU[nVtxBin]->Fill(l1TPhw_->hcalTPieta[tpIt],l1TPhw_->hcalTPet[tpIt]);
0504       // }
0505 
0506     }
0507 
0508     //fill emulated tower histos
0509     hNTowemuVsNVtx->Fill(nVtx,nTowemu);
0510     hProfNVtxNTowemu->Fill(nVtx,nTowemu);
0511     hNTow_emu->Fill(nTowemu);
0512     hNTow4_emu->Fill(nTowemu4);
0513     hNTowemuVsNVtx4->Fill(nVtx,nTowemu4);
0514     hProfNVtxNTowemu4->Fill(nVtx,nTowemu4);
0515     hProfNTowemuNVtx4->Fill(nTowemu4,nVtx);
0516 
0517   }
0518      
0519   //end event loop, now plot histos  
0520   TFile outFile( outFilename.c_str() , "new");
0521 
0522   TCanvas* canvas = new TCanvas("canvas","",600,600);
0523   TLegend* legend = new TLegend(0.5,0.6,0.7,0.88);
0524   
0525   legend->SetLineColor(0);
0526 
0527   stringstream saveName("");
0528  
0529   formatPlot1D(hNTow_emu,4); 
0530   formatPlot1D(hNTow4_emu,4); 
0531   formatPlot1D(hNVtx,4);
0532   formatPlot2D(hTowEtaPhi);
0533 
0534   hTowEtaPhi->Draw("colz");
0535   hTowEtaPhi->Write();
0536   canvas->SaveAs("towEtaPhi.pdf");
0537 
0538   formatPlot2D(hNTowemuVsNVtx);
0539  
0540   hNTow_emu->Draw();
0541   hNTow_emu->Write();
0542   canvas->SaveAs("nTowemu.pdf");
0543 
0544   hNTow4_emu->Draw();
0545   hNTow4_emu->Write();
0546   canvas->SaveAs("nTowemu4.pdf");
0547 
0548   formatPlot1D(hAllTowEtemu,4);
0549   hAllTowEtemu->Draw();
0550   hAllTowEtemu->Write();
0551   canvas->SaveAs("towEtEmu.pdf");
0552 
0553   formatPlot1D(hAllTowEtaemu,4);
0554   hAllTowEtaemu->Draw();
0555   hAllTowEtaemu->Write();
0556   canvas->SaveAs("towEtaEmu.pdf");
0557   
0558   hNVtx->Draw();
0559   hNVtx->Write();
0560   canvas->SaveAs("nVtx.pdf");
0561 
0562   hNTowemuVsNVtx->Draw("colz");
0563   hNTowemuVsNVtx->Write();
0564   canvas->SaveAs("nTowemuVsNVtx.pdf");
0565 
0566   hNTowemuVsNVtx4->Draw("colz");
0567   hNTowemuVsNVtx4->Write();
0568   canvas->SaveAs("nTowemuVsNVtx4.pdf");
0569 
0570   hProfNVtxNTowemu->SetMarkerStyle(7);
0571   hProfNVtxNTowemu->Draw("");
0572   hProfNVtxNTowemu->Write();
0573   canvas->SaveAs("ProfNVtxNTowemu.pdf");
0574 
0575   hProfNVtxNTowemu4->SetMarkerStyle(7);  
0576   hProfNVtxNTowemu4->Draw("");
0577   hProfNVtxNTowemu4->Write();
0578   canvas->SaveAs("ProfNVtxNTowemu4.pdf");
0579 
0580   fitProfile(hProfNTowemuNVtx4);
0581 
0582   if(!doTow){ 
0583 
0584     formatPlot1D(hTowEt,4);
0585     //formatPlot1D(hECALTPEt,4);
0586     //formatPlot1D(hHCALTPEt,4);
0587     
0588     hProfTowEtEta->SetMarkerStyle(7);
0589     hProfTowEtEta->SetMaximum(4.0);
0590     hProfTowEtEta->Scale(0.5);
0591     hProfTowEtEta->Draw("");
0592     hProfTowEtEta->Write();
0593     canvas->SaveAs("ProfTowEtEta.pdf");
0594       
0595 
0596     for(uint i=0;i<hProfTowEtEtaPU.size();i++){
0597       hProfTowEtEtaPU[i]->SetMarkerStyle(7);
0598       hProfTowEtEtaPU[i]->SetMaximum(4.0);
0599       hProfTowEtEtaPU[i]->Scale(0.5);
0600       hProfTowEtEtaPU[i]->Draw("");
0601       hProfTowEtEtaPU[i]->Write();
0602       stringstream fn("");
0603       fn << "ProfTowEtEtaPU_" << i << ".pdf";
0604       //canvas->SaveAs(fn.str().c_str());
0605     }
0606     hProfECALTPEtEta->SetMarkerStyle(7);
0607     hProfECALTPEtEta->SetMaximum(4.0);
0608     hProfECALTPEtEta->Draw("");
0609     hProfECALTPEtEta->Write();
0610     //canvas->SaveAs("ProfECALTPEtEta.pdf");
0611 
0612     for(uint i=0;i<hProfECALTPEtEtaPU.size();i++){
0613       hProfECALTPEtEtaPU[i]->SetMarkerStyle(7);
0614       hProfECALTPEtEtaPU[i]->SetMaximum(4.0);
0615       hProfECALTPEtEtaPU[i]->Draw("");
0616       hProfECALTPEtEtaPU[i]->Write();
0617       stringstream fn("");
0618       fn << "ProfECALTPEtEtaPU_" << i << ".pdf";
0619       //canvas->SaveAs(fn.str().c_str());
0620     }
0621 
0622     hProfHCALTPEtEta->SetMarkerStyle(7);
0623     hProfHCALTPEtEta->SetMaximum(4.0);
0624     hProfHCALTPEtEta->Draw("");
0625     hProfHCALTPEtEta->Write();
0626     //canvas->SaveAs("ProfHCALTPEtEta.pdf");
0627 
0628     for(uint i=0;i<hProfHCALTPEtEtaPU.size();i++){
0629       hProfHCALTPEtEtaPU[i]->SetMarkerStyle(7);
0630       hProfHCALTPEtEtaPU[i]->SetMaximum(4.0);
0631       hProfHCALTPEtEtaPU[i]->Draw("");
0632       hProfHCALTPEtEtaPU[i]->Write();
0633       stringstream fn("");
0634       fn << "ProfHCALTPEtEtaPU_" << i << ".pdf";
0635       //canvas->SaveAs(fn.str().c_str());
0636     }
0637   
0638     canvas->SetLogy(1);
0639 
0640     hTowEt->Draw("");
0641     hTowEt->Write("");
0642     canvas->SaveAs("TowEt.pdf");
0643 
0644     hECALTPEt->Draw("");
0645     hECALTPEt->Write("");
0646     //canvas->SaveAs("ECALTPEt.pdf");
0647 
0648     hHCALTPEt->Draw("");
0649     hHCALTPEt->Write("");
0650     //canvas->SaveAs("HCALTPEt.pdf");
0651 
0652   
0653     for(uint pu=0;pu<puBins;pu++){
0654       for(uint eta=0; eta<41; ++eta){
0655       
0656         stringstream etaS("");
0657         etaS << "_Eta" << eta+1;
0658       
0659     formatPlot1D(hTowEtPU[pu][eta],4);
0660     hTowEtPU[pu][eta]->GetXaxis()->SetLimits(0.,50.);
0661     hTowEtPU[pu][eta]->Draw("");
0662     hTowEtPU[pu][eta]->Write();
0663     stringstream fn("");
0664     fn << "TowEtPU_" << pu << etaS.str() << ".pdf";
0665     //canvas->SaveAs(fn.str().c_str());
0666 
0667     // for(uint pu=0;pu<hECALTPEtPU.size();pu++){
0668     //   formatPlot1D(hECALTPEtPU[pu][eta],4);
0669     //   hECALTPEtPU[pu][eta]->Draw("");
0670     //   hECALTPEtPU[pu][eta]->Write();
0671     //   stringstream fn("");
0672     //   fn << "ECALTPEtPU_" << pu << etaS.str() << ".pdf";
0673     //   //canvas->SaveAs(fn.str().c_str());
0674     // }
0675      
0676     // for(uint pu=0;pu<hHCALTPEtPU.size();pu++){
0677     //   formatPlot1D(hHCALTPEtPU[pu][eta],4);
0678     //   hHCALTPEtPU[pu][eta]->Draw("");
0679     // hHCALTPEtPU[pu][eta]->Write();
0680     // stringstream fn("");
0681     // fn << "HCALTPEtPU_" << pu << etaS.str() << ".pdf";
0682     // //canvas->SaveAs(fn.str().c_str());
0683     // }
0684       }
0685     }
0686     
0687     
0688     
0689     
0690     canvas->Close();
0691     outFile.Close();
0692     
0693   }
0694 }