Analysis2010/macros/analyzeLoop.C

#include "TNtuple.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TCanvas.h"
#include <iostream>

using namespace std;

#define MAXHITS 100000


class MyRecHit{
public:
   int n;

   float e[MAXHITS];
   float et[MAXHITS];
   float eta[MAXHITS];
   float phi[MAXHITS];
   bool isjet[MAXHITS];

   float jtpt;
   float jteta;
   float jtphi;

   int depth;

};

class MyBkg{
public:
   int n;
   float rho[50];
   float sigma[50];
};


void analyzeLoop(const char* infile = "/d101/yetkin/work/zdcScint/r150476.root", const char* outfile = "plotsRun150476.root"){

   cout<<"a"<<endl;

   const char* bit0 = "HLT_HIMinBiasHF";
   const char* bit1 = "L1_BscMinBiasOR_BptxAND";

   double towerCut = 3;

   bool MC = false;
   bool OpenHLTonly = true;

   TFile * inf = new TFile(infile);
   TFile* outf = new TFile(outfile,"recreate");

   TNtuple* nt = new TNtuple("evNT","","ETmyTowers:NmyTowers:EThfRecHit:NhfRecHit:NgenParts:ETgenParts");

   TH1D* ha[100];
   TH1D* hb[100];
   TH1D* hc[100];
   TH1D* he[100];
   TH1D* hn[100];

   TH1D* hneb[100];
   TH1D* hnee[100];

   TH2D* h2o[100];
   TH2D* h2a[100];
   TH2D* h2b[100];
   TH2D* h2c[100];
   TH2D* h2d[100];
   TH2D* h2pm[100];
   TH2D* h2ebhf[100];
   TH2D* h2z[100];
   TH2D* h2zpm[100];
   TH2D* h2zz[100];
   TH2D* h2rNpart[100];
   TH2D* h2sNpart[100];
   TH2D* h2rsNpart[100];

   TH2D* h2pmN[100];

   TH2D* h2rHF[100];
   TH2D* h2sHF[100];
   TH2D* h2rsHF[100];

   for(int i = 0; i < 20; ++i){
      ha[i] = new TH1D(Form("ha%d",i),";N HF towers above throshold;Events",100,0,1000);
      hb[i] = new TH1D(Form("hb%d",i),";Sum Tower E_{T} [GeV];Events",100,0,4000);
      hc[i] = new TH1D(Form("hc%d",i),";Sum HF E_{T} [GeV];Events",100,0,5000);


      he[i] = new TH1D(Form("he%d",i),"",1000,0,4000);

      h2z[i] = new TH2D(Form("h2z%d",i),";Sum HF E_{T} [GeV]; Sum ZDC Energy [GeV]",100,0,5000,100,0,350000);
      h2zpm[i] = new TH2D(Form("h2zpm%d",i),";Sum ZDC- Energy [GeV]; Sum ZDC+ Energy [GeV]",100,-100,200000,100,-100,200000);
      h2zz[i] = new TH2D(Form("h2zz%d",i),";Sum HF E_{T} [GeV]; Sum ZDC Energy [GeV]",100,0,500,100,0,350000);

      hn[i] = new TH1D(Form("hn%d",i),"",1000,0,10000);

      hneb[i] = new TH1D(Form("hneb%d",i),";N_{hits} EE > 0.3 GeV; Events",100,0,10000);
      hnee[i] = new TH1D(Form("hnee%d",i),";N_{hits} EE > 0.3 GeV; Events",100,0,10000);

      h2o[i] = new TH2D(Form("h2o%d",i),"",420,0,420,600,0,6000);
      h2a[i] = new TH2D(Form("h2a%d",i),"",420,0,420,900,0,9000);
      h2b[i] = new TH2D(Form("h2b%d",i),"",900,0,9000,600,0,6000);
      h2c[i] = new TH2D(Form("h2c%d",i),"",420,0,420,370,0,3700);
      h2d[i] = new TH2D(Form("h2d%d",i),"",370,0,3700,600,0,6000);

      h2pm[i] = new TH2D(Form("h2pm%d",i),"; E_{T} Sum HF- [GeV]; E_{T} Sum HF+ [GeV]",100,0,2500,100,0,2500);
      h2ebhf[i] = new TH2D(Form("h2ebhf%d",i),"; E_{T} Sum HF [GeV]; E_{T} Sum EB [GeV]",100,0,5000,100,0,500);

      h2rHF[i] = new TH2D(Form("h2rHF%d",i),Form("; E_{T} Sum HF [GeV]; #rho[%d] [GeV]",i),100,0,5000,100,0,500);
      h2sHF[i] = new TH2D(Form("h2sHF%d",i),Form("; E_{T} Sum HF [GeV]; #sigma[%d] [GeV]",i),100,0,5000,100,0,500);
      h2rsHF[i] = new TH2D(Form("h2rsHF%d",i),Form("; E_{T} Sum HF [GeV]; #rho[%d] + #sigma[%d] [GeV]",i,i),100,0,5000,100,0,500);

      h2rNpart[i] = new TH2D(Form("h2rNpart%d",i),Form("; N_{part}; #rho[%d] [GeV]",i,i),100,0,500,100,0,500);
      h2sNpart[i] = new TH2D(Form("h2sNpart%d",i),Form("; N_{part}; #sigma[%d] [GeV]",i,i),100,0,500,100,0,500);
      h2rsNpart[i] = new TH2D(Form("h2rsNpart%d",i),Form("; N_{part}; #rho[%d] + #sigma[%d] [GeV]",i,i),100,0,500,100,0,500);

      h2pmN[i] = new TH2D(Form("h2pmN%d",i),"; N towers HF-; N towers HF+",1000,0,1000,1000,0,1000);

   }

   ha[1]->SetLineColor(2);
   ha[1]->SetMarkerColor(2);

   TTree* t1 = (TTree*)inf->Get("hltanalysis/HltTree");
   TTree* t2 = (TTree*)inf->Get("rechits/hf");
   TTree* t3 = (TTree*)inf->Get("rechits/hbhe");
   TTree* t4 = (TTree*)inf->Get("rechits/ee");
   TTree* t5 = (TTree*)inf->Get("rechits/eb");
   TTree* t6 = (TTree*)inf->Get("rechits/bkg");
   TTree* t7 = (TTree*)inf->Get("rechits/tower");
   TTree* t8;
   if(MC){
      t8 = (TTree*)inf->Get("mc/hi");
      t1->AddFriend(t8);
   }

   cout<<"a"<<endl;

   t1->AddFriend(t2);
   t1->AddFriend(t3);
   t1->AddFriend(t4);
   t1->AddFriend(t5);
   t1->AddFriend(t6);
   t1->AddFriend(t7);
   cout<<"a"<<endl;

   MyRecHit hbheRecHit;
   MyRecHit hfRecHit;
   MyRecHit ebRecHit;
   MyRecHit eeRecHit;
   MyRecHit myBC;
   MyRecHit myTowers;
   MyRecHit genParts;
   MyBkg bkg;
   cout<<"a"<<endl;


   int hiBin;
   int hiNpix, hiNtracks, hiNtracksPtCut, hiNtracksEtaCut, hiNtracksEtaPtCut;
   float hiHF, hiHFplus, hiHFminus, hiHFhit, hiHFhitPlus, hiHFhitMinus, hiEB, hiET, hiEE, hiEEplus, hiEEminus, hiZDC, hiZDCplus, hiZDCminus;
   cout<<"a"<<endl;
  
   float Npart;
   float Ncoll;
   float Nhard;
   float Phi0;
   float b;
   cout<<"a"<<endl;

   int fNcharged;
   int fNchargedMR;
   float fMeanPt;
   float fMeanPtMR;
   float fEtMR;
   int fNchargedPtCut;
   int fNchargedPtCutMR;
   cout<<"a"<<endl;

   int trig[20];
   cout<<"a"<<endl;

   t1->SetBranchAddress(bit0,&(trig[0]));
   t1->SetBranchAddress(bit1,&(trig[1]));
   cout<<"a"<<endl;

   t1->SetBranchAddress("Npart",&Npart);
   t1->SetBranchAddress("Ncoll",&Ncoll);
   t1->SetBranchAddress("hiHF",&hiHF);
   t1->SetBranchAddress("hiHFplus",&hiHFplus);
   t1->SetBranchAddress("hiHFminus",&hiHFminus);
   t1->SetBranchAddress("hiZDC",&hiZDC);
   t1->SetBranchAddress("hiZDCplus",&hiZDCplus);
   t1->SetBranchAddress("hiZDCminus",&hiZDCminus);
   cout<<"a"<<endl;

   if(!OpenHLTonly){
   t6->SetBranchAddress("n",&bkg.n);
   t6->SetBranchAddress("rho",bkg.rho);
   t6->SetBranchAddress("sigma",bkg.sigma);

   t3->SetBranchAddress("e",hbheRecHit.e);
   t3->SetBranchAddress("et",hbheRecHit.et);
   t3->SetBranchAddress("eta",hbheRecHit.eta);
   t3->SetBranchAddress("phi",hbheRecHit.phi);
   t3->SetBranchAddress("n",&hbheRecHit.n);

   t2->SetBranchAddress("e",hfRecHit.e);
   t2->SetBranchAddress("et",hfRecHit.et);
   t2->SetBranchAddress("eta",hfRecHit.eta);
   t2->SetBranchAddress("phi",hfRecHit.phi);
   t2->SetBranchAddress("n",&hfRecHit.n);

   t5->SetBranchAddress("e",ebRecHit.e);
   t5->SetBranchAddress("et",ebRecHit.et);
   t5->SetBranchAddress("eta",ebRecHit.eta);
   t5->SetBranchAddress("phi",ebRecHit.phi);
   t5->SetBranchAddress("n",&ebRecHit.n);

   t4->SetBranchAddress("e",eeRecHit.e);
   t4->SetBranchAddress("et",eeRecHit.et);
   t4->SetBranchAddress("eta",eeRecHit.eta);
   t4->SetBranchAddress("phi",eeRecHit.phi);
   t4->SetBranchAddress("n",&eeRecHit.n);

   t7->SetBranchAddress("e",myTowers.e);
   t7->SetBranchAddress("et",myTowers.et);
   t7->SetBranchAddress("eta",myTowers.eta);
   t7->SetBranchAddress("phi",myTowers.phi);
   t7->SetBranchAddress("n",&myTowers.n);

   if(MC){
   t8->SetBranchAddress("pt",genParts.et);
   t8->SetBranchAddress("eta",genParts.eta);
   t8->SetBranchAddress("phi",genParts.phi);
   t8->SetBranchAddress("mult",&genParts.n);
   }
   }

   cout<<"a"<<endl;

   for(int iev = 0; iev < t1->GetEntries(); ++iev){
      t1->GetEntry(iev);
      if(!OpenHLTonly){

        t2->GetEntry(iev);
        t3->GetEntry(iev);
        t4->GetEntry(iev);
        t5->GetEntry(iev);
        t6->GetEntry(iev);
        t7->GetEntry(iev);
        if(MC)t8->GetEntry(iev);
      }


      cout<<"a"<<endl;

      if(!trig[0]) continue;

      cout<<"a"<<endl;

      float EThbheRecHit = 0;
      for(int i = 0; i < hbheRecHit.n; ++i){
         float e = hbheRecHit.e[i];
         float et = hbheRecHit.et[i];
         float eta = hbheRecHit.eta[i];
         float phi = hbheRecHit.phi[i];
         EThbheRecHit += et;
      }

      float ETmyTowers = 0;
      int NmyTowers = 0;
      int NmyTowersP = 0;
      int NmyTowersM = 0;
      int NmyTowersT = 0;
      for(int i = 0; i < myTowers.n; ++i){
         float e = myTowers.e[i];
         float et = myTowers.et[i];
         float eta = myTowers.eta[i];
         float phi = myTowers.phi[i];

         ETmyTowers += et;
         if(fabs(eta) < 2.87) continue;
         if(e < towerCut) continue;
         NmyTowers++;
         if(eta > 0 ) NmyTowersP++;
         if(eta < 0 ) NmyTowersM++;

      }

      float EThfRecHit = 0;
      int NhfRecHit = 0;

      for(int i = 0; i < hfRecHit.n; ++i){
         float e = hfRecHit.e[i];
         float et = hfRecHit.et[i];
         float eta = hfRecHit.eta[i];
         float phi = hfRecHit.phi[i];
         EThfRecHit += et;
         if(e < towerCut) continue;
         NhfRecHit++;
      }

      float ETgenParts = 0;
      int NgenParts = 0;

      for(int i = 0; i < genParts.n; ++i){
         float e = genParts.e[i];
         float et = genParts.et[i];
         float eta = genParts.eta[i];
         float phi = genParts.phi[i];

         if(fabs(eta) < 2.87) continue;
         NgenParts++;
         ETgenParts += et;
      }

      ha[0]->Fill(NmyTowers);
      h2a[0]->Fill(Npart,NgenParts);
      h2b[0]->Fill(NgenParts,hiHF);
      h2c[0]->Fill(Npart,ETgenParts);
      h2d[0]->Fill(ETgenParts,hiHF);
      
      if(!trig[0] && trig[1]) h2pmN[1]->Fill(NmyTowersM,NmyTowersP);
      h2pmN[0]->Fill(NmyTowersM,NmyTowersP);


      h2pm[0]->Fill(hiHFminus,hiHFplus);
      h2ebhf[0]->Fill(hiHF,hiEB);
      h2z[0]->Fill(hiHF,hiZDC);
      h2zz[0]->Fill(hiHF,hiZDC);
      h2o[0]->Fill(Npart,hiHF);

      hneb[0]->Fill(0);
      hnee[0]->Fill(0);
      hb[0]->Fill(ETmyTowers);

      for(int i = 0; i < bkg.n; ++i){
         h2rHF[i]->Fill(hiHF,bkg.rho[i]);
         h2sHF[i]->Fill(hiHF,bkg.sigma[i]);
         h2rsHF[i]->Fill(hiHF,bkg.rho[i]+bkg.sigma[i]);
         h2rNpart[i]->Fill(Npart,bkg.rho[i]);
         h2sNpart[i]->Fill(Npart,bkg.sigma[i]);
         h2rsNpart[i]->Fill(Npart,bkg.rho[i]+bkg.sigma[i]);
      }

      nt->Fill(ETmyTowers,NmyTowers,EThfRecHit,NhfRecHit,NgenParts,ETgenParts);
   }


   cout<<"Main loop complete"<<endl;

   t1->AddFriend(nt);


   TCanvas* c1 = new TCanvas("c1","",1200,800);
   c1->Divide(5,3);

   c1->cd(1);

   c1->cd(2);
  c1->cd(3);
  c1->cd(4);
  c1->cd(5);

  c1->cd(6);
  h2rsNpart[1]->Draw("colz");

  c1->cd(7);
  h2rsNpart[6]->Draw("colz");

  c1->cd(8);
  //  t1->Draw("hiEB:Npart",trig,"colz");

  c1->cd(9);
  //  t1->Draw("hiNpix:Npart",trig,"colz");

  c1->cd(10);
  //  t1->Draw(Form("NmyTowers>>ha0"));
  //  t1->Draw(Form("NmyTowers>>ha1"),trig,"same");
  ha[0]->Draw();

  c1->cd(11);
  h2a[0]->Draw("colz");
  
  c1->cd(12);
  h2b[0]->Draw("colz");

  c1->cd(2);

  c1->cd(3);

  c1->cd(4);
  h2c[0]->Draw("colz");

  c1->cd(5);
  h2d[0]->Draw("colz");

  outf->Write();

}

Revision:	1.1
Committed:	Fri Nov 12 17:03:15 2010 UTC (14 years, 5 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Log Message:	macros
#	User	Rev	Content
1	yilmaz	1.1	#include "TNtuple.h"
2			#include "TFile.h"
3			#include "TH1D.h"
4			#include "TH2D.h"
5			#include "TCanvas.h"
6			#include <iostream>
7
8			using namespace std;
9
10			#define MAXHITS 100000
11
12
13			class MyRecHit{
14			public:
15			int n;
16
17			float e[MAXHITS];
18			float et[MAXHITS];
19			float eta[MAXHITS];
20			float phi[MAXHITS];
21			bool isjet[MAXHITS];
22
23			float jtpt;
24			float jteta;
25			float jtphi;
26
27			int depth;
28
29			};
30
31			class MyBkg{
32			public:
33			int n;
34			float rho[50];
35			float sigma[50];
36			};
37
38
39			void analyzeLoop(const char* infile = "/d101/yetkin/work/zdcScint/r150476.root", const char* outfile = "plotsRun150476.root"){
40
41			cout<<"a"<<endl;
42
43			const char* bit0 = "HLT_HIMinBiasHF";
44			const char* bit1 = "L1_BscMinBiasOR_BptxAND";
45
46			double towerCut = 3;
47
48			bool MC = false;
49			bool OpenHLTonly = true;
50
51			TFile * inf = new TFile(infile);
52			TFile* outf = new TFile(outfile,"recreate");
53
54			TNtuple* nt = new TNtuple("evNT","","ETmyTowers:NmyTowers:EThfRecHit:NhfRecHit:NgenParts:ETgenParts");
55
56			TH1D* ha[100];
57			TH1D* hb[100];
58			TH1D* hc[100];
59			TH1D* he[100];
60			TH1D* hn[100];
61
62			TH1D* hneb[100];
63			TH1D* hnee[100];
64
65			TH2D* h2o[100];
66			TH2D* h2a[100];
67			TH2D* h2b[100];
68			TH2D* h2c[100];
69			TH2D* h2d[100];
70			TH2D* h2pm[100];
71			TH2D* h2ebhf[100];
72			TH2D* h2z[100];
73			TH2D* h2zpm[100];
74			TH2D* h2zz[100];
75			TH2D* h2rNpart[100];
76			TH2D* h2sNpart[100];
77			TH2D* h2rsNpart[100];
78
79			TH2D* h2pmN[100];
80
81			TH2D* h2rHF[100];
82			TH2D* h2sHF[100];
83			TH2D* h2rsHF[100];
84
85			for(int i = 0; i < 20; ++i){
86			ha[i] = new TH1D(Form("ha%d",i),";N HF towers above throshold;Events",100,0,1000);
87			hb[i] = new TH1D(Form("hb%d",i),";Sum Tower E_{T} [GeV];Events",100,0,4000);
88			hc[i] = new TH1D(Form("hc%d",i),";Sum HF E_{T} [GeV];Events",100,0,5000);
89
90
91			he[i] = new TH1D(Form("he%d",i),"",1000,0,4000);
92
93			h2z[i] = new TH2D(Form("h2z%d",i),";Sum HF E_{T} [GeV]; Sum ZDC Energy [GeV]",100,0,5000,100,0,350000);
94			h2zpm[i] = new TH2D(Form("h2zpm%d",i),";Sum ZDC- Energy [GeV]; Sum ZDC+ Energy [GeV]",100,-100,200000,100,-100,200000);
95			h2zz[i] = new TH2D(Form("h2zz%d",i),";Sum HF E_{T} [GeV]; Sum ZDC Energy [GeV]",100,0,500,100,0,350000);
96
97			hn[i] = new TH1D(Form("hn%d",i),"",1000,0,10000);
98
99			hneb[i] = new TH1D(Form("hneb%d",i),";N_{hits} EE > 0.3 GeV; Events",100,0,10000);
100			hnee[i] = new TH1D(Form("hnee%d",i),";N_{hits} EE > 0.3 GeV; Events",100,0,10000);
101
102			h2o[i] = new TH2D(Form("h2o%d",i),"",420,0,420,600,0,6000);
103			h2a[i] = new TH2D(Form("h2a%d",i),"",420,0,420,900,0,9000);
104			h2b[i] = new TH2D(Form("h2b%d",i),"",900,0,9000,600,0,6000);
105			h2c[i] = new TH2D(Form("h2c%d",i),"",420,0,420,370,0,3700);
106			h2d[i] = new TH2D(Form("h2d%d",i),"",370,0,3700,600,0,6000);
107
108			h2pm[i] = new TH2D(Form("h2pm%d",i),"; E_{T} Sum HF- [GeV]; E_{T} Sum HF+ [GeV]",100,0,2500,100,0,2500);
109			h2ebhf[i] = new TH2D(Form("h2ebhf%d",i),"; E_{T} Sum HF [GeV]; E_{T} Sum EB [GeV]",100,0,5000,100,0,500);
110
111			h2rHF[i] = new TH2D(Form("h2rHF%d",i),Form("; E_{T} Sum HF [GeV]; #rho[%d] [GeV]",i),100,0,5000,100,0,500);
112			h2sHF[i] = new TH2D(Form("h2sHF%d",i),Form("; E_{T} Sum HF [GeV]; #sigma[%d] [GeV]",i),100,0,5000,100,0,500);
113			h2rsHF[i] = new TH2D(Form("h2rsHF%d",i),Form("; E_{T} Sum HF [GeV]; #rho[%d] + #sigma[%d] [GeV]",i,i),100,0,5000,100,0,500);
114
115			h2rNpart[i] = new TH2D(Form("h2rNpart%d",i),Form("; N_{part}; #rho[%d] [GeV]",i,i),100,0,500,100,0,500);
116			h2sNpart[i] = new TH2D(Form("h2sNpart%d",i),Form("; N_{part}; #sigma[%d] [GeV]",i,i),100,0,500,100,0,500);
117			h2rsNpart[i] = new TH2D(Form("h2rsNpart%d",i),Form("; N_{part}; #rho[%d] + #sigma[%d] [GeV]",i,i),100,0,500,100,0,500);
118
119			h2pmN[i] = new TH2D(Form("h2pmN%d",i),"; N towers HF-; N towers HF+",1000,0,1000,1000,0,1000);
120
121			}
122
123			ha[1]->SetLineColor(2);
124			ha[1]->SetMarkerColor(2);
125
126			TTree* t1 = (TTree*)inf->Get("hltanalysis/HltTree");
127			TTree* t2 = (TTree*)inf->Get("rechits/hf");
128			TTree* t3 = (TTree*)inf->Get("rechits/hbhe");
129			TTree* t4 = (TTree*)inf->Get("rechits/ee");
130			TTree* t5 = (TTree*)inf->Get("rechits/eb");
131			TTree* t6 = (TTree*)inf->Get("rechits/bkg");
132			TTree* t7 = (TTree*)inf->Get("rechits/tower");
133			TTree* t8;
134			if(MC){
135			t8 = (TTree*)inf->Get("mc/hi");
136			t1->AddFriend(t8);
137			}
138
139			cout<<"a"<<endl;
140
141			t1->AddFriend(t2);
142			t1->AddFriend(t3);
143			t1->AddFriend(t4);
144			t1->AddFriend(t5);
145			t1->AddFriend(t6);
146			t1->AddFriend(t7);
147			cout<<"a"<<endl;
148
149			MyRecHit hbheRecHit;
150			MyRecHit hfRecHit;
151			MyRecHit ebRecHit;
152			MyRecHit eeRecHit;
153			MyRecHit myBC;
154			MyRecHit myTowers;
155			MyRecHit genParts;
156			MyBkg bkg;
157			cout<<"a"<<endl;
158
159
160			int hiBin;
161			int hiNpix, hiNtracks, hiNtracksPtCut, hiNtracksEtaCut, hiNtracksEtaPtCut;
162			float hiHF, hiHFplus, hiHFminus, hiHFhit, hiHFhitPlus, hiHFhitMinus, hiEB, hiET, hiEE, hiEEplus, hiEEminus, hiZDC, hiZDCplus, hiZDCminus;
163			cout<<"a"<<endl;
164
165			float Npart;
166			float Ncoll;
167			float Nhard;
168			float Phi0;
169			float b;
170			cout<<"a"<<endl;
171
172			int fNcharged;
173			int fNchargedMR;
174			float fMeanPt;
175			float fMeanPtMR;
176			float fEtMR;
177			int fNchargedPtCut;
178			int fNchargedPtCutMR;
179			cout<<"a"<<endl;
180
181			int trig[20];
182			cout<<"a"<<endl;
183
184			t1->SetBranchAddress(bit0,&(trig[0]));
185			t1->SetBranchAddress(bit1,&(trig[1]));
186			cout<<"a"<<endl;
187
188			t1->SetBranchAddress("Npart",&Npart);
189			t1->SetBranchAddress("Ncoll",&Ncoll);
190			t1->SetBranchAddress("hiHF",&hiHF);
191			t1->SetBranchAddress("hiHFplus",&hiHFplus);
192			t1->SetBranchAddress("hiHFminus",&hiHFminus);
193			t1->SetBranchAddress("hiZDC",&hiZDC);
194			t1->SetBranchAddress("hiZDCplus",&hiZDCplus);
195			t1->SetBranchAddress("hiZDCminus",&hiZDCminus);
196			cout<<"a"<<endl;
197
198			if(!OpenHLTonly){
199			t6->SetBranchAddress("n",&bkg.n);
200			t6->SetBranchAddress("rho",bkg.rho);
201			t6->SetBranchAddress("sigma",bkg.sigma);
202
203			t3->SetBranchAddress("e",hbheRecHit.e);
204			t3->SetBranchAddress("et",hbheRecHit.et);
205			t3->SetBranchAddress("eta",hbheRecHit.eta);
206			t3->SetBranchAddress("phi",hbheRecHit.phi);
207			t3->SetBranchAddress("n",&hbheRecHit.n);
208
209			t2->SetBranchAddress("e",hfRecHit.e);
210			t2->SetBranchAddress("et",hfRecHit.et);
211			t2->SetBranchAddress("eta",hfRecHit.eta);
212			t2->SetBranchAddress("phi",hfRecHit.phi);
213			t2->SetBranchAddress("n",&hfRecHit.n);
214
215			t5->SetBranchAddress("e",ebRecHit.e);
216			t5->SetBranchAddress("et",ebRecHit.et);
217			t5->SetBranchAddress("eta",ebRecHit.eta);
218			t5->SetBranchAddress("phi",ebRecHit.phi);
219			t5->SetBranchAddress("n",&ebRecHit.n);
220
221			t4->SetBranchAddress("e",eeRecHit.e);
222			t4->SetBranchAddress("et",eeRecHit.et);
223			t4->SetBranchAddress("eta",eeRecHit.eta);
224			t4->SetBranchAddress("phi",eeRecHit.phi);
225			t4->SetBranchAddress("n",&eeRecHit.n);
226
227			t7->SetBranchAddress("e",myTowers.e);
228			t7->SetBranchAddress("et",myTowers.et);
229			t7->SetBranchAddress("eta",myTowers.eta);
230			t7->SetBranchAddress("phi",myTowers.phi);
231			t7->SetBranchAddress("n",&myTowers.n);
232
233			if(MC){
234			t8->SetBranchAddress("pt",genParts.et);
235			t8->SetBranchAddress("eta",genParts.eta);
236			t8->SetBranchAddress("phi",genParts.phi);
237			t8->SetBranchAddress("mult",&genParts.n);
238			}
239			}
240
241			cout<<"a"<<endl;
242
243			for(int iev = 0; iev < t1->GetEntries(); ++iev){
244			t1->GetEntry(iev);
245			if(!OpenHLTonly){
246
247			t2->GetEntry(iev);
248			t3->GetEntry(iev);
249			t4->GetEntry(iev);
250			t5->GetEntry(iev);
251			t6->GetEntry(iev);
252			t7->GetEntry(iev);
253			if(MC)t8->GetEntry(iev);
254			}
255
256
257
258			cout<<"a"<<endl;
259
260			if(!trig[0]) continue;
261
262			cout<<"a"<<endl;
263
264			float EThbheRecHit = 0;
265			for(int i = 0; i < hbheRecHit.n; ++i){
266			float e = hbheRecHit.e[i];
267			float et = hbheRecHit.et[i];
268			float eta = hbheRecHit.eta[i];
269			float phi = hbheRecHit.phi[i];
270			EThbheRecHit += et;
271			}
272
273			float ETmyTowers = 0;
274			int NmyTowers = 0;
275			int NmyTowersP = 0;
276			int NmyTowersM = 0;
277			int NmyTowersT = 0;
278			for(int i = 0; i < myTowers.n; ++i){
279			float e = myTowers.e[i];
280			float et = myTowers.et[i];
281			float eta = myTowers.eta[i];
282			float phi = myTowers.phi[i];
283
284			ETmyTowers += et;
285			if(fabs(eta) < 2.87) continue;
286			if(e < towerCut) continue;
287			NmyTowers++;
288			if(eta > 0 ) NmyTowersP++;
289			if(eta < 0 ) NmyTowersM++;
290
291			}
292
293			float EThfRecHit = 0;
294			int NhfRecHit = 0;
295
296			for(int i = 0; i < hfRecHit.n; ++i){
297			float e = hfRecHit.e[i];
298			float et = hfRecHit.et[i];
299			float eta = hfRecHit.eta[i];
300			float phi = hfRecHit.phi[i];
301			EThfRecHit += et;
302			if(e < towerCut) continue;
303			NhfRecHit++;
304			}
305
306			float ETgenParts = 0;
307			int NgenParts = 0;
308
309			for(int i = 0; i < genParts.n; ++i){
310			float e = genParts.e[i];
311			float et = genParts.et[i];
312			float eta = genParts.eta[i];
313			float phi = genParts.phi[i];
314
315			if(fabs(eta) < 2.87) continue;
316			NgenParts++;
317			ETgenParts += et;
318			}
319
320			ha[0]->Fill(NmyTowers);
321			h2a[0]->Fill(Npart,NgenParts);
322			h2b[0]->Fill(NgenParts,hiHF);
323			h2c[0]->Fill(Npart,ETgenParts);
324			h2d[0]->Fill(ETgenParts,hiHF);
325
326			if(!trig[0] && trig[1]) h2pmN[1]->Fill(NmyTowersM,NmyTowersP);
327			h2pmN[0]->Fill(NmyTowersM,NmyTowersP);
328
329
330
331			h2pm[0]->Fill(hiHFminus,hiHFplus);
332			h2ebhf[0]->Fill(hiHF,hiEB);
333			h2z[0]->Fill(hiHF,hiZDC);
334			h2zz[0]->Fill(hiHF,hiZDC);
335			h2o[0]->Fill(Npart,hiHF);
336
337			hneb[0]->Fill(0);
338			hnee[0]->Fill(0);
339			hb[0]->Fill(ETmyTowers);
340
341			for(int i = 0; i < bkg.n; ++i){
342			h2rHF[i]->Fill(hiHF,bkg.rho[i]);
343			h2sHF[i]->Fill(hiHF,bkg.sigma[i]);
344			h2rsHF[i]->Fill(hiHF,bkg.rho[i]+bkg.sigma[i]);
345			h2rNpart[i]->Fill(Npart,bkg.rho[i]);
346			h2sNpart[i]->Fill(Npart,bkg.sigma[i]);
347			h2rsNpart[i]->Fill(Npart,bkg.rho[i]+bkg.sigma[i]);
348			}
349
350			nt->Fill(ETmyTowers,NmyTowers,EThfRecHit,NhfRecHit,NgenParts,ETgenParts);
351			}
352
353
354			cout<<"Main loop complete"<<endl;
355
356			t1->AddFriend(nt);
357
358
359			TCanvas* c1 = new TCanvas("c1","",1200,800);
360			c1->Divide(5,3);
361
362			c1->cd(1);
363
364			c1->cd(2);
365			c1->cd(3);
366			c1->cd(4);
367			c1->cd(5);
368
369			c1->cd(6);
370			h2rsNpart[1]->Draw("colz");
371
372			c1->cd(7);
373			h2rsNpart[6]->Draw("colz");
374
375			c1->cd(8);
376			// t1->Draw("hiEB:Npart",trig,"colz");
377
378			c1->cd(9);
379			// t1->Draw("hiNpix:Npart",trig,"colz");
380
381			c1->cd(10);
382			// t1->Draw(Form("NmyTowers>>ha0"));
383			// t1->Draw(Form("NmyTowers>>ha1"),trig,"same");
384			ha[0]->Draw();
385
386			c1->cd(11);
387			h2a[0]->Draw("colz");
388
389			c1->cd(12);
390			h2b[0]->Draw("colz");
391
392			c1->cd(2);
393
394			c1->cd(3);
395
396			c1->cd(4);
397			h2c[0]->Draw("colz");
398
399			c1->cd(5);
400			h2d[0]->Draw("colz");
401
402			outf->Write();
403
404			}
405