Corrections implemented into AliAnalysisHadEtReconstructed, simulated hadronic and...
[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisHadEtMonteCarlo.cxx
1 //_________________________________________________________________________
2 //  Utility Class for transverse energy studies, charged hadrons
3 //  Base class for MC analysis
4 //  - MC output
5 // implementation file
6 //
7 //Created by Christine Nattrass, Rebecca Scott, Irakli Martashvili
8 //University of Tennessee at Knoxville
9 //_________________________________________________________________________
10 #include "AliAnalysisHadEtMonteCarlo.h"
11 #include "AliAnalysisEtCuts.h"
12
13 #include "AliStack.h"
14 #include "AliMCEvent.h"
15 #include "AliESDEvent.h"
16 #include "AliESDtrackCuts.h"
17 #include "AliESDpid.h"
18 #include <iostream>
19
20 using namespace std;
21
22 ClassImp(AliAnalysisHadEtMonteCarlo);
23
24 AliAnalysisHadEtMonteCarlo::AliAnalysisHadEtMonteCarlo():AliAnalysisHadEt()
25                             ,fSimHadEt(0)
26                             ,fSimTotEt(0) 
27 {
28 }
29 void AliAnalysisHadEtMonteCarlo::ResetEventValues(){
30   AliAnalysisHadEt::ResetEventValues();
31     fSimHadEt=0.0;
32     fSimTotEt=0.0;
33 }
34 Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
35 { // analyse MC and real event info
36   FillHisto1D("NEvents",0.5,1);
37
38   AnalyseEvent(ev);
39   AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
40   AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
41   AliStack *stack = mcEvent->Stack();
42
43   //for PID
44   AliESDpid *pID = new AliESDpid();
45   pID->MakePID(realEvent);
46
47   //This code taken from https://twiki.cern.ch/twiki/bin/view/ALICE/SelectionOfPrimaryTracksForPp2009DataAnalysis
48   //Gets good tracks
49   //=============================================
50   // Primary vertex
51   const AliESDVertex *vertex = realEvent->GetPrimaryVertexTracks();
52   if(vertex->GetNContributors()<1) {
53     // SPD vertex
54     vertex = realEvent->GetPrimaryVertexSPD();
55     if(vertex->GetNContributors()<1) {
56       // NO GOOD VERTEX, SKIP EVENT 
57     }
58   }
59   // apply a cut |zVertex| < CUT, if needed
60
61   
62
63   //fEsdtrackCutsITSTPC->SetEtaRange(-0.8,0.8); // normally, |eta|<0.8
64   //=============================================
65
66   //Roughly following $ALICE_ROOT/PWG0/dNdEta/AlidNdEtaCorrectionTask
67
68   //=============================================TPC&&ITS=============================================
69   TString *strTPC = new TString("TPC");
70   TString *strITS = new TString("ITS");
71   TString *strTPCITS = new TString("TPCITS");
72   for(Int_t cutset=0;cutset<3;cutset++){
73     TString *cutName;
74     TObjArray* list;
75     switch(cutset){
76     case 0:
77       cutName = strTPC;
78       list = fEsdtrackCutsTPC->GetAcceptedTracks(realEvent);
79       break;
80     case 1:
81       cutName = strITS;
82       list = fEsdtrackCutsITS->GetAcceptedTracks(realEvent);
83       break;
84     case 2:
85       cutName = strTPCITS;
86       list = fEsdtrackCutsITSTPC->GetAcceptedTracks(realEvent);
87       break;
88     default:
89       cerr<<"Error:  cannot fill histograms!"<<endl;
90       return -1;
91     }
92     Int_t nGoodTracks = list->GetEntries();
93     for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
94       {
95         AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
96         if (!track)
97           {
98             Printf("ERROR: Could not get track %d", iTrack);
99             continue;
100           }
101         else{
102           Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
103           if(cutset!=1){
104             nSigmaPion = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kPion));
105             nSigmaProton = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kProton));
106             nSigmaKaon = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kKaon));
107             nSigmaElectron = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kElectron));
108           }
109           else{
110             nSigmaPion = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kPion));
111             nSigmaProton = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kProton));
112             nSigmaKaon = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kKaon));
113             nSigmaElectron = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kElectron));
114           }
115           bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
116           bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
117           bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
118           bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
119
120           //bool IsElectron = false;
121           bool unidentified = (!isProton && !isKaon && !isElectron);
122           Float_t dEdx = track->GetTPCsignal();
123           if(cutset==1) dEdx = track->GetITSsignal();
124
125           FillHisto2D(Form("dEdxAll%s",cutName->Data()),track->P(),dEdx,1.0);
126           //if(cutset==1) cout<<"filling "<<track->P()<<" "<<dEdx<<endl;
127
128           UInt_t label = (UInt_t)TMath::Abs(track->GetLabel());
129           TParticle  *simPart  = stack->Particle(label);
130           if(!simPart) {
131             Printf("no MC particle\n");                 
132             continue;           
133           }
134           else{//analysis
135             if(stack->IsPhysicalPrimary(label)){
136               if (TMath::Abs(simPart->Eta()) < fCuts->GetCommonEtaCut())            {
137
138                 Int_t pdgCode =  simPart->GetPDG(0)->PdgCode();
139                 Int_t mypid = 0;
140                 if(pdgCode==fPiPlusCode) mypid = 1;
141                 if(pdgCode==fProtonCode) mypid = 2;
142                 if(pdgCode==fKPlusCode) mypid = 3;
143                 if(pdgCode==fEPlusCode) mypid = 4;
144                 if(pdgCode==fPiMinusCode) mypid = 1;
145                 if(pdgCode==fAntiProtonCode) mypid = 2;
146                 if(pdgCode==fKMinusCode) mypid = 3;
147                 if(pdgCode==fEMinusCode) mypid = 4;
148                 //cout<<pdgCode->PdgCode()<<" ";
149                 //fPdgDB->GetSimParticle("pi+")->PdgCode();
150                 bool filled = false;      
151                 //============Charged hadrons===================================
152                 //identified...
153                 if(isPion){
154                   if(pdgCode!=fPiPlusCode && pdgCode!=fPiMinusCode){
155                     FillHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),1,mypid,1);
156                     //if(mypid==0)cerr<<"I was misidentified! I'm not a pion! I am a "<<simPart->GetName()<<endl;
157                   }
158                   float myEt = Et(simPart);
159                   if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedPiPlus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
160                   else{ FillHisto2D(Form("EtReconstructed%sIdentifiedPiMinus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
161                   FillHisto2D(Form("dEdxPion%s",cutName->Data()),track->P(),dEdx,1.0);
162                 }
163                 if(isProton){
164                   if(pdgCode!=fProtonCode && pdgCode!=fAntiProtonCode){
165                     FillHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),2,mypid,1);
166                     // if(mypid==0)cerr<<"I was misidentified!  I'm not a proton! I am a "<<simPart->GetName()<<endl;
167                   }
168                   float myEt = Et(simPart);
169                   if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedProton",cutName->Data()),track->Pt(),track->Eta(),myEt);}
170                   else{ FillHisto2D(Form("EtReconstructed%sIdentifiedAntiProton",cutName->Data()),track->Pt(),track->Eta(),myEt);}
171                   FillHisto2D(Form("dEdxProton%s",cutName->Data()),track->P(),dEdx,1.0);
172                 }
173                 if(isKaon){
174                   if(pdgCode!=fKMinusCode && pdgCode!=fKPlusCode){
175                     FillHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),3,mypid,1);
176                     //if(mypid==0)cerr<<"I was misidentified!  I'm not a kaon! I am a "<<simPart->GetName()<<" p "<<track->P()<<" nSigmaProton "<<nSigmaProton<<" nSigmaPion "<<nSigmaPion<<" nSigmaKaon "<<nSigmaKaon<<endl;
177                   }
178                   float myEt = Et(simPart);
179                   if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedKPlus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
180                   else{ FillHisto2D(Form("EtReconstructed%sIdentifiedKMinus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
181                   FillHisto2D(Form("dEdxKaon%s",cutName->Data()),track->P(),dEdx,1.0);
182                 }
183                 if(isElectron){
184                   if(pdgCode!=fEMinusCode && pdgCode!=fEPlusCode){
185                     FillHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),4,mypid,1);
186                     //cerr<<"I was misidentified!  I'm not an electron! I am a "<<simPart->GetName()<<endl;
187                   }
188                   float myEt = Et(simPart);
189                   if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedEPlus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
190                   else{ FillHisto2D(Form("EtReconstructed%sIdentifiedEMinus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
191                   FillHisto2D(Form("dEdxElectron%s",cutName->Data()),track->P(),dEdx,1.0);
192                 }
193                 if(unidentified){
194                   if(pdgCode!=fEMinusCode && pdgCode!=fEPlusCode){
195                     float myEtPi = Et(simPart,fPionMass);
196                     float myEt = Et(simPart);
197                     FillHisto2D(Form("EtReconstructed%sUnidentifiedAssumingPion",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
198                     FillHisto2D(Form("EtReconstructed%sUnidentified",cutName->Data()),track->Pt(),track->Eta(),myEt);
199                     FillHisto2D(Form("EtNReconstructed%sUnidentified",cutName->Data()),track->Pt(),track->Eta(),1.0);
200                   }
201                   FillHisto2D(Form("dEdxUnidentified%s",cutName->Data()),track->P(),dEdx,1.0);
202                   //cout<<"I was not identified.  I am a "<<simPart->GetName()<<" PID "<<pdgCode<<endl;
203                   //track what was not identified successfully
204                   FillHisto1D(Form("UnidentifiedPIDs%s",cutName->Data()),mypid,1);
205                 }
206                 //...simulated
207                 if(pdgCode == fPiPlusCode){
208                   //cout<<"I'm a real primary "<<simPart->GetName()<<"! "<<"my label is "<<simPart->GetFirstMother()<<" track no "<<iTrack<<"/"<<realEvent->GetNumberOfTracks()<<endl;//<<" "<<label<<" "<<pdgCode<<endl;
209                 
210                   float myEt = Et(simPart);
211                   FillHisto2D(Form("EtReconstructed%sPiPlus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
212                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
213                   FillHisto2D(Form("EtNReconstructed%sPiPlus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
214                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
215                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
216                   filled = true;
217                 }
218                 if(pdgCode == fPiMinusCode){
219                   float myEt = Et(simPart);
220                   FillHisto2D(Form("EtReconstructed%sPiMinus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
221                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
222                   FillHisto2D(Form("EtNReconstructed%sPiMinus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
223                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
224                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
225                   filled = true;
226                 }
227                 if(pdgCode == fKPlusCode){
228                   float myEt = Et(simPart);
229                   float myEtPi = Et(simPart,fPionMass);
230                   FillHisto2D(Form("EtReconstructed%sKPlus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
231                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
232                   FillHisto2D(Form("EtNReconstructed%sKPlus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
233                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
234                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEtPi);
235                   FillHisto2D(Form("EtReconstructed%sKPlusAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEtPi);
236                   filled = true;
237                 }
238                 if(pdgCode == fKMinusCode){
239                   float myEt = Et(simPart);
240                   float myEtPi = Et(simPart,fPionMass);
241                   FillHisto2D(Form("EtReconstructed%sKMinus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
242                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
243                   FillHisto2D(Form("EtNReconstructed%sKMinus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
244                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
245                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEtPi);
246                   FillHisto2D(Form("EtReconstructed%sKMinusAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEtPi);
247                   filled = true;
248                 }
249                 if(pdgCode == fProtonCode){
250                   float myEt = Et(simPart);
251                   float myEtPi = Et(simPart,fPionMass);
252                   FillHisto2D(Form("EtReconstructed%sProton",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
253                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
254                   FillHisto2D(Form("EtNReconstructed%sProton",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
255                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
256                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEtPi);
257                   FillHisto2D(Form("EtReconstructed%sProtonAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEtPi);
258                   filled = true;
259                 }
260                 if(pdgCode == fAntiProtonCode){
261                   float myEt = Et(simPart);
262                   float myEtPi = Et(simPart,fPionMass);
263                   FillHisto2D(Form("EtReconstructed%sAntiProton",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
264                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
265                   FillHisto2D(Form("EtNReconstructed%sAntiProton",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
266                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
267                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEtPi);
268                   FillHisto2D(Form("EtReconstructed%sAntiProtonAssumingPion",cutName->Data()),simPart->Pt(),simPart->Eta(),myEtPi);
269                   filled = true;
270                 }
271                 if(pdgCode == fEPlusCode){
272                   float myEt = Et(simPart);
273                   FillHisto2D(Form("EtReconstructed%sEPlus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
274                   if(!isElectron || unidentified){
275                     float myEtPi = Et(simPart,fPionMass);
276                     FillHisto2D(Form("EtReconstructed%sMisidentifiedElectrons",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
277                   }
278                   filled = true;
279                 }
280                 if(pdgCode == fEMinusCode){
281                   if(!isElectron || unidentified){
282                     float myEtPi = Et(simPart,fPionMass);
283                     FillHisto2D(Form("EtReconstructed%sMisidentifiedElectrons",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
284                   }
285                   float myEt = Et(simPart);
286                     FillHisto2D(Form("EtReconstructed%sEMinus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
287                     filled = true;
288                 }
289                 //if(!filled){
290                   //TParticlePDG *pc = simPart->GetPDG(0);
291                   //if( strcmp(pc->ParticleClass(),"Baryon")==0 || strcmp(pc->ParticleClass(),"Meson")==0 ){
292                   //cout<<"Did not find a place for "<<simPart->GetName()<<" "<<pdgCode<<" which is a "<<pc->ParticleClass()<<endl;
293                   //}
294                   //}
295               }
296               
297             }
298             else{//not a primary - we're after V0 daughters!
299               //cout<<"I'm a secondary "<<simPart->GetName()<<"!";//<<endl;
300               TParticle *mom = stack->Particle(simPart->GetFirstMother());
301               if(mom){
302                 TParticlePDG *pc = mom->GetPDG(0);
303                 if(pc){
304                   Int_t pdgCode =  mom->GetPDG(0)->PdgCode();
305                   if(pdgCode == fLambdaCode){
306                     float myEt = Et(simPart);
307                     FillHisto2D(Form("EtReconstructed%sLambdaDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
308                   }
309                   if(pdgCode == fAntiLambdaCode){
310                     float myEt = Et(simPart);
311                     FillHisto2D(Form("EtReconstructed%sAntiLambdaDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
312                   }
313                   if(pdgCode == fK0SCode){
314                     float myEt = Et(simPart);
315                     FillHisto2D(Form("EtReconstructed%sK0SDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
316                   }
317                   if(pdgCode == fXiCode){
318                     float myEt = Et(simPart);
319                     FillHisto2D(Form("EtReconstructed%sXiDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
320                   }
321                   if(pdgCode == fAntiXiCode){
322                     float myEt = Et(simPart);
323                     FillHisto2D(Form("EtReconstructed%sAntiXiDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
324                   }
325                   if(pdgCode == fOmegaCode){
326                     float myEt = Et(simPart);
327                     FillHisto2D(Form("EtReconstructed%sOmegaDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
328                   }
329                   if(pdgCode == fXiCode){
330                     float myEt = Et(simPart);
331                     FillHisto2D(Form("EtReconstructed%sAntiOmegaDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
332                   }
333
334                   if(mom->GetFirstMother()>0){
335                     TParticle *grandma = stack->Particle(mom->GetFirstMother());
336                     if(grandma){
337                       Int_t pdgCodeMom =  mom->GetPDG(0)->PdgCode();
338                       if(pdgCodeMom==fPiPlusCode || pdgCodeMom==fPiMinusCode || pdgCodeMom==fProtonCode ||pdgCodeMom==fAntiProtonCode || pdgCodeMom==fKPlusCode || pdgCode==fKMinusCode){
339                         //cout<<" my grandmother is "<<grandma->GetName()<<" "<<endl;
340                         Int_t pdgCodeGrandma =  grandma->GetPDG(0)->PdgCode();
341                       
342                         if(pdgCodeGrandma == fXiCode){
343                           float myEt = Et(simPart);
344                           FillHisto2D(Form("EtReconstructed%sXiDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
345                         }
346                         if(pdgCodeGrandma == fAntiXiCode){
347                           float myEt = Et(simPart);
348                           FillHisto2D(Form("EtReconstructed%sAntiXiDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
349                         }
350                         if(pdgCodeGrandma == fOmegaCode){
351                           float myEt = Et(simPart);
352                           FillHisto2D(Form("EtReconstructed%sOmegaDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
353                         }
354                         if(pdgCodeGrandma == fXiCode){
355                           float myEt = Et(simPart);
356                           FillHisto2D(Form("EtReconstructed%sAntiOmegaDaughters",cutName->Data()),track->Pt(),track->Eta(),myEt);
357                         }
358
359                       }
360                     }
361                   }
362                 }
363               }
364             }
365           }
366
367         }
368       }
369     delete list;
370   }
371   //delete AliESDpid;
372   return 1;
373 }
374 Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
375 { // analyse MC event
376      ResetEventValues();
377      
378     // Get us an mc event
379     AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
380
381     // Let's play with the stack!
382     AliStack *stack = mcEvent->Stack();
383
384     Int_t nPrim = stack->GetNtrack();
385
386     //=================Tracks which may or may not have been reconstructed=================
387
388     for (Int_t iPart = 0; iPart < nPrim; iPart++)
389     {
390
391         TParticle *part = stack->Particle(iPart);
392
393         if (!part)
394           {
395             Printf("ERROR: Could not get particle %d", iPart);
396             continue;
397           }
398
399         TParticlePDG *pc = part->GetPDG(0);
400
401         // Check if it is a primary particle
402         if (stack->IsPhysicalPrimary(iPart)){//primaries
403
404           if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())           {
405
406             Int_t pdgCode =  part->GetPDG(0)->PdgCode();
407             //cout<<pdgCode->PdgCode()<<" ";
408             //fPdgDB->GetParticle("pi+")->PdgCode();
409             bool filled = false;
410             //============Charged hadrons===================================
411             if(pdgCode == fPiPlusCode){
412               //cout<<"I'm a simulated primary "<<part->GetName()<<"! "<<"my label is "<<part->GetFirstMother()<<" pt "<<part->Pt()<<endl;
413               float myEt = Et(part);
414               fSimHadEt += myEt;
415               fSimTotEt += myEt;
416               FillHisto2D("EtSimulatedPiPlus",part->Pt(),part->Eta(),myEt);
417               FillHisto2D("EtNSimulatedPiPlus",part->Pt(),part->Eta(),1.0);
418               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
419               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
420               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEt);
421               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
422               filled = true;
423             }
424             if(pdgCode == fPiMinusCode){
425               float myEt = Et(part);
426               fSimHadEt += myEt;
427               fSimTotEt += myEt;
428               FillHisto2D("EtSimulatedPiMinus",part->Pt(),part->Eta(),myEt);
429               FillHisto2D("EtNSimulatedPiMinus",part->Pt(),part->Eta(),1.0);
430               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
431               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
432               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEt);
433               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
434               filled = true;
435             }
436             if(pdgCode == fKPlusCode){
437               float myEt = Et(part);
438               float myEtPi = Et(part,fPionMass);
439               fSimHadEt += myEt;
440               fSimTotEt += myEt;
441               FillHisto2D("EtSimulatedKPlus",part->Pt(),part->Eta(),myEt);
442               FillHisto2D("EtNSimulatedKPlus",part->Pt(),part->Eta(),1.0);
443               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
444               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
445               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEtPi);
446               FillHisto2D("EtSimulatedKPlusAssumingPion",part->Pt(),part->Eta(),myEtPi);
447               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
448               filled = true;
449             }
450             if(pdgCode == fKMinusCode){
451               float myEt = Et(part);
452               float myEtPi = Et(part,fPionMass);
453               fSimHadEt += myEt;
454               fSimTotEt += myEt;
455               FillHisto2D("EtSimulatedKMinus",part->Pt(),part->Eta(),myEt);
456               FillHisto2D("EtNSimulatedKMinus",part->Pt(),part->Eta(),1.0);
457               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
458               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
459               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEtPi);
460               FillHisto2D("EtSimulatedKMinusAssumingPion",part->Pt(),part->Eta(),myEtPi);
461               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
462               filled = true;
463             }
464             if(pdgCode == fProtonCode){
465               float myEt = Et(part);
466               float myEtPi = Et(part,fPionMass);
467               fSimHadEt += myEt;
468               fSimTotEt += myEt;
469               FillHisto2D("EtSimulatedProton",part->Pt(),part->Eta(),myEt);
470               FillHisto2D("EtNSimulatedProton",part->Pt(),part->Eta(),1.0);
471               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
472               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
473               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEtPi);
474               FillHisto2D("EtSimulatedProtonAssumingPion",part->Pt(),part->Eta(),myEtPi);
475               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
476               filled = true;
477             }
478             if(pdgCode == fAntiProtonCode){
479               float myEt = Et(part);
480               float myEtPi = Et(part,fPionMass);
481               fSimHadEt += myEt;
482               fSimTotEt += myEt;
483               FillHisto2D("EtSimulatedAntiProton",part->Pt(),part->Eta(),myEt);
484               FillHisto2D("EtNSimulatedAntiProton",part->Pt(),part->Eta(),1.0);
485               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
486               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
487               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEtPi);
488               FillHisto2D("EtSimulatedAntiProtonAssumingPion",part->Pt(),part->Eta(),myEtPi);
489               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
490               filled = true;
491             }
492             //============Other hadrons===================================
493
494             if(pdgCode == fNeutronCode){
495               float myEt = Et(part);
496               fSimHadEt += myEt;
497               fSimTotEt += myEt;
498               FillHisto2D("EtSimulatedNeutron",part->Pt(),part->Eta(),myEt);
499               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
500               filled = true;
501             }
502             if(pdgCode == fAntiNeutronCode){
503               float myEt = Et(part);
504               fSimHadEt += myEt;
505               fSimTotEt += myEt;
506               FillHisto2D("EtSimulatedAntiNeutron",part->Pt(),part->Eta(),myEt);
507               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
508               filled = true;
509             }
510             if(pdgCode == fLambdaCode){
511               float myEt = Et(part);
512               fSimHadEt += myEt;
513               fSimTotEt += myEt;
514               //cout<<"I am a simulated lambda! pt "<<part->Pt()<<" eta "<<part->Eta()<<endl;
515               FillHisto2D("EtSimulatedLambda",part->Pt(),part->Eta(),myEt);
516               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
517               Int_t ndaughters = part->GetNDaughters();
518               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
519                 Int_t daughterindex = part->GetDaughter(idaughter);
520                 if(daughterindex<0 || daughterindex>1e5) continue;
521                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
522                 if(daughter){
523                   if(daughter->GetPDG(0)){
524                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
525                     if(daughtercode==fPiMinusCode || daughtercode==fProtonCode){
526                       myEt = Et(daughter);
527                       FillHisto2D("EtSimulatedLambdaDaughters",daughter->Pt(),daughter->Eta(),myEt);
528                       //cout<<"Lambda daughter is a "<<daughter->GetName()<<endl;
529                     }
530                   }
531                   else{
532                     //cout<<"Lambda daughter is a "<<daughter->GetName()<<endl;
533                   }
534                 }
535               }
536               filled = true;
537             }
538             if(pdgCode == fAntiLambdaCode){
539               float myEt = Et(part);
540               fSimHadEt += myEt;
541               fSimTotEt += myEt;
542               FillHisto2D("EtSimulatedAntiLambda",part->Pt(),part->Eta(),myEt);
543               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
544               Int_t ndaughters = part->GetNDaughters();
545               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
546                 Int_t daughterindex = part->GetDaughter(idaughter);
547                 if(daughterindex<0 || daughterindex>1e5) continue;
548                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
549                 if(daughter){
550                   if(daughter->GetPDG(0)){
551                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
552                     if(daughtercode==fPiPlusCode || daughtercode==fAntiProtonCode){
553                       myEt = Et(daughter);
554                       FillHisto2D("EtSimulatedAntiLambdaDaughters",daughter->Pt(),daughter->Eta(),myEt);
555                       //cout<<"AntiLambda daughter is a "<<daughter->GetName()<<endl;
556                     }
557                   }
558                   else{
559                     //cout<<"AntiLambda daughter is a "<<daughter->GetName()<<endl;
560                   }
561                 }
562               }
563               filled = true;
564             }
565             if(pdgCode == fK0SCode){
566               float myEt = Et(part);
567               fSimHadEt += myEt;
568               fSimTotEt += myEt;
569               FillHisto2D("EtSimulatedK0S",part->Pt(),part->Eta(),myEt);
570               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
571               Int_t ndaughters = part->GetNDaughters();
572               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
573                 Int_t daughterindex = part->GetDaughter(idaughter);
574                 if(daughterindex<0 || daughterindex>1e5) continue;
575                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
576                 if(daughter){
577                   if(daughter->GetPDG(0)){
578
579                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
580                     if(daughtercode==fPiMinusCode || daughtercode==fPiPlusCode){
581                       myEt = Et(daughter);
582                       FillHisto2D("EtSimulatedK0SDaughters",daughter->Pt(),daughter->Eta(),myEt);
583                       //cout<<"K0S daughter is a "<<daughter->GetName()<<endl;
584                     }
585                   }
586                   else{
587                     //cout<<"K0S daughter is a "<<daughter->GetName()<<endl;
588                   }
589                 }
590               }
591               filled = true;
592             }
593             if(pdgCode == fK0LCode){
594               float myEt = Et(part);
595               fSimHadEt += myEt;
596               fSimTotEt += myEt;
597               FillHisto2D("EtSimulatedK0L",part->Pt(),part->Eta(),myEt);
598               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
599               filled = true;
600             }
601             if(pdgCode == fOmegaCode){
602               float myEt = Et(part);
603               fSimHadEt += myEt;
604               fSimTotEt += myEt;
605               FillHisto2D("EtSimulatedOmega",part->Pt(),part->Eta(),myEt);
606               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
607               Int_t ndaughters = part->GetNDaughters();
608               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
609                 Int_t daughterindex = part->GetDaughter(idaughter);
610                 if(daughterindex<0 || daughterindex>1e5) continue;
611                 TParticle *daughter = stack->Particle(daughterindex);
612                 if(daughter){
613                   if(daughter->GetPDG(0)){
614
615                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
616                     if(daughtercode==fPiPlusCode || daughtercode==fProtonCode || daughtercode==fKMinusCode){
617                       myEt = Et(daughter);
618                       FillHisto2D("EtSimulatedOmegaDaughters",daughter->Pt(),daughter->Eta(),myEt);
619                     //cout<<"Omega daughter is a "<<daughter->GetName()<<endl;
620                     }
621                   }
622                   else{
623                     //cout<<"Omega daughter is a "<<daughter->GetName()<<endl;
624                   }
625                 }
626               }
627               filled = true;
628             }
629             if(pdgCode == fAntiOmegaCode){
630               float myEt = Et(part);
631               fSimHadEt += myEt;
632               fSimTotEt += myEt;
633               FillHisto2D("EtSimulatedOmega",part->Pt(),part->Eta(),myEt);
634               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
635               Int_t ndaughters = part->GetNDaughters();
636               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
637                 Int_t daughterindex = part->GetDaughter(idaughter);
638                 if(daughterindex<0 || daughterindex>1e5) continue;
639                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
640                 if(daughter){
641                   if(daughter->GetPDG(0)){
642                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
643                     if(daughtercode==fPiMinusCode || daughtercode==fAntiProtonCode || daughtercode==fKPlusCode){
644                       myEt = Et(daughter);
645                       FillHisto2D("EtSimulatedAntiOmegaDaughters",daughter->Pt(),daughter->Eta(),myEt);
646                       //cout<<"AntiOmega daughter is a "<<daughter->GetName()<<endl;
647                     }
648                   }
649                   else{
650                     //cout<<"AntiOmega daughter is a "<<daughter->GetName()<<endl;
651                   }
652                 }
653               }
654               filled = true;
655             }
656             //There are two codes for Sigmas
657             if(pdgCode == fSigmaCode || pdgCode == -3222){
658               float myEt = Et(part);
659               fSimHadEt += myEt;
660               fSimTotEt += myEt;
661               FillHisto2D("EtSimulatedSigma",part->Pt(),part->Eta(),myEt);
662               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
663               filled = true;
664             }
665             if(pdgCode == fAntiSigmaCode || pdgCode == 3222){
666               float myEt = Et(part);
667               fSimHadEt += myEt;
668               fSimTotEt += myEt;
669               FillHisto2D("EtSimulatedAntiSigma",part->Pt(),part->Eta(),myEt);
670               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
671               filled = true;
672             }
673             if(pdgCode == fXiCode){
674               float myEt = Et(part);
675               fSimHadEt += myEt;
676               fSimTotEt += myEt;
677               FillHisto2D("EtSimulatedXi",part->Pt(),part->Eta(),myEt);
678               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
679               Int_t ndaughters = part->GetNDaughters();
680               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
681                 Int_t daughterindex = part->GetDaughter(idaughter);
682                 if(daughterindex<0 || daughterindex>1e5 || daughterindex>1e5) continue;
683                 //cerr<<"Daughter index "<<daughterindex<<" npart "<<nPrim<<endl;
684                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
685                 if(daughter){
686                   if(daughter->GetPDG(0)){
687
688                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
689                     if(daughtercode==fPiPlusCode || daughtercode==fProtonCode || daughtercode==fPiMinusCode){
690                       myEt = Et(daughter);
691                       FillHisto2D("EtSimulatedXiDaughters",daughter->Pt(),daughter->Eta(),myEt);
692                     //cout<<"Xi daughter is a "<<daughter->GetName()<<endl;
693                     }
694                   }
695                   else{
696                     //cout<<"Xi daughter is a "<<daughter->GetName()<<endl;
697                   }
698                 }
699               }
700               filled = true;
701             }
702             if(pdgCode == fAntiXiCode){
703               float myEt = Et(part);
704               fSimHadEt += myEt;
705               fSimTotEt += myEt;
706               FillHisto2D("EtSimulatedAntiXi",part->Pt(),part->Eta(),myEt);
707               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
708               Int_t ndaughters = part->GetNDaughters();
709               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
710                 Int_t daughterindex = part->GetDaughter(idaughter);
711                 if(daughterindex<0 || daughterindex>1e5) continue;
712                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
713                 if(daughter){
714                   if(daughter->GetPDG(0)){
715                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
716                     if(daughtercode==fPiPlusCode || daughtercode==fAntiProtonCode || daughtercode==fPiMinusCode){
717                       myEt = Et(daughter);
718                       FillHisto2D("EtSimulatedAntiXiDaughters",daughter->Pt(),daughter->Eta(),myEt);
719                       //cout<<"AntiXi daughter is a "<<daughter->GetName()<<endl;
720                     }
721                   }
722                   else{
723                     //cout<<"AntiXi daughter is a "<<daughter->GetName()<<endl;
724                   }
725                 }
726               }
727               filled = true;
728             }
729             if(pdgCode == fXi0Code){
730               float myEt = Et(part);
731               fSimHadEt += myEt;
732               fSimTotEt += myEt;
733               FillHisto2D("EtSimulatedXi0",part->Pt(),part->Eta(),myEt);
734               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
735               filled = true;
736             }
737             if(pdgCode == fAntiXi0Code){
738               float myEt = Et(part);
739               fSimHadEt += myEt;
740               fSimTotEt += myEt;
741               FillHisto2D("EtSimulatedAntiXi0",part->Pt(),part->Eta(),myEt);
742               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
743               filled = true;
744             }
745             //============electrons===================================
746
747             if(pdgCode == fEPlusCode){
748               float myEt = Et(part);
749               fSimTotEt += myEt;
750               FillHisto2D("EtSimulatedEPlus",part->Pt(),part->Eta(),myEt);
751               filled = true;
752             }
753             if(pdgCode == fEMinusCode){
754               float myEt = Et(part);
755               fSimTotEt += myEt;
756               FillHisto2D("EtSimulatedEMinus",part->Pt(),part->Eta(),myEt);
757               filled = true;
758             }
759             //============neutrals===================================
760             if(pdgCode == fGammaCode){
761               TParticle *mom = stack->Particle(part->GetFirstMother());
762               Int_t pdgCodeMom =  mom->GetPDG(0)->PdgCode();
763               //cout<<"I am a gamma and my mom is "<<mom->GetName()<<endl;
764               if(pdgCodeMom == fEtaCode){
765                 float myEt = Et(mom);
766               fSimTotEt += myEt;
767                 FillHisto2D("EtSimulatedEta",mom->Pt(),mom->Eta(),myEt);
768                 filled = true;
769               }
770               if(pdgCodeMom == fPi0Code){
771                 float myEt = Et(mom);
772               fSimTotEt += myEt;
773                 FillHisto2D("EtSimulatedPi0",mom->Pt(),mom->Eta(),myEt);
774                 filled = true;
775               }
776               if(pdgCodeMom == fOmega0Code){
777                 float myEt = Et(mom);
778               fSimTotEt += myEt;
779                 FillHisto2D("EtSimulatedOmega0",mom->Pt(),mom->Eta(),myEt);
780                 filled = true;
781               }
782               if(!filled){
783                 float myEt = Et(part);
784               fSimTotEt += myEt;
785                 FillHisto2D("EtSimulatedGamma",part->Pt(),part->Eta(),myEt);
786                 filled = true;
787               }
788             }
789             if(pdgCode == fEtaCode){
790               float myEt = Et(part);
791               fSimTotEt += myEt;
792               FillHisto2D("EtSimulatedEta",part->Pt(),part->Eta(),myEt);
793               filled = true;
794             }
795             if(pdgCode == fPi0Code){
796               float myEt = Et(part);
797               fSimTotEt += myEt;
798               FillHisto2D("EtSimulatedPi0",part->Pt(),part->Eta(),myEt);
799               filled = true;
800             }
801             if(pdgCode == fOmega0Code){
802               float myEt = Et(part);
803               fSimTotEt += myEt;
804               FillHisto2D("EtSimulatedOmega0",part->Pt(),part->Eta(),myEt);
805               filled = true;
806             }
807             if(!filled){
808               //if( strcmp(pc->ParticleClass(),"Baryon")==0 || strcmp(pc->ParticleClass(),"Meson")==0 ){
809                 cout<<"Did not find a place for "<<part->GetName()<<" "<<pdgCode<<" which is a "<<pc->ParticleClass()<<endl;
810                 //}
811             }
812           }
813         }
814     }
815
816
817
818
819 //     fTotNeutralEtAcc = fTotNeutralEt;
820 //     fTotEt = fTotChargedEt + fTotNeutralEt;
821 //     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
822     
823 //     FillHistograms();
824
825     return 1;
826     
827 }
828
829 void AliAnalysisHadEtMonteCarlo::Init()
830 { // Init
831     AliAnalysisHadEt::Init();
832 }
833
834 void AliAnalysisHadEtMonteCarlo::CreateHistograms(){
835   //for simulated Et only (no reconstruction)
836   CreateEtaPtHisto2D(TString("EtSimulatedPiPlus"),TString("Simulated E_{T} from #pi^{+}"));
837   CreateEtaPtHisto2D("EtSimulatedPiMinus","Simulated E_{T} from #pi^{-}");
838   CreateEtaPtHisto2D("EtSimulatedKPlus","Simulated E_{T} from K^{+}");
839   CreateEtaPtHisto2D("EtSimulatedKMinus","Simulated E_{T} from K^{-}");
840   CreateEtaPtHisto2D("EtSimulatedProton","Simulated E_{T} from p");
841   CreateEtaPtHisto2D("EtSimulatedAntiProton","Simulated E_{T} from #bar{p}");
842   CreateEtaPtHisto2D("EtSimulatedChargedHadron","Simulated E_{T} from charged hadrons");
843   CreateEtaPtHisto2D(TString("EtNSimulatedPiPlus"),TString("Number of Simulated #pi^{+}"));
844   CreateEtaPtHisto2D("EtNSimulatedPiMinus","Number of simulated #pi^{-}");
845   CreateEtaPtHisto2D("EtNSimulatedKPlus","Number of simulated K^{+}");
846   CreateEtaPtHisto2D("EtNSimulatedKMinus","Number of simulated K^{-}");
847   CreateEtaPtHisto2D("EtNSimulatedProton","Number of simulated p");
848   CreateEtaPtHisto2D("EtNSimulatedAntiProton","Number of simulated #bar{p}");
849   CreateEtaPtHisto2D("EtNSimulatedChargedHadron","Number of simulated charged hadrons");
850
851   CreateEtaPtHisto2D("EtSimulatedChargedHadronAssumingPion","Simulated E_{T} from charged hadrons assuming they are all pions");
852   CreateEtaPtHisto2D("EtSimulatedKPlusAssumingPion","Simulated E_{T} from K^{+} assuming #pi mass");
853   CreateEtaPtHisto2D("EtSimulatedKMinusAssumingPion","Simulated E_{T} from K^{-} assuming #pi mass");
854   CreateEtaPtHisto2D("EtSimulatedProtonAssumingPion","Simulated E_{T} from p assuming #pi mass");
855   CreateEtaPtHisto2D("EtSimulatedAntiProtonAssumingPion","Simulated E_{T} from #bar{p} assuming #pi mass");
856
857   CreateEtaPtHisto2D("EtSimulatedLambda","Simulated E_{T} from #Lambda");
858   CreateEtaPtHisto2D("EtSimulatedAntiLambda","Simulated E_{T} from #bar{#Lambda}");
859   CreateEtaPtHisto2D("EtSimulatedK0S","Simulated E_{T} from K^{0}_{S}");
860   CreateEtaPtHisto2D("EtSimulatedK0L","Simulated E_{T} from K^{0}_{L}");
861   CreateEtaPtHisto2D("EtSimulatedNeutron","Simulated E_{T} from neutrons");
862   CreateEtaPtHisto2D("EtSimulatedAntiNeutron","Simulated E_{T} from #bar{n}");
863   CreateEtaPtHisto2D("EtSimulatedEPlus","Simulated E_{T} from e^{+}");
864   CreateEtaPtHisto2D("EtSimulatedEMinus","Simulated E_{T} from e^{-}");
865   CreateEtaPtHisto2D("EtSimulatedOmega","Simulated E_{T} from #Omega^{-}");
866   CreateEtaPtHisto2D("EtSimulatedAntiOmega","Simulated E_{T} from #Omega^{+}");
867   CreateEtaPtHisto2D("EtSimulatedXi","Simulated E_{T} from #Xi^{-}");
868   CreateEtaPtHisto2D("EtSimulatedAntiXi","Simulated E_{T} from #Xi^{+}");
869   CreateEtaPtHisto2D("EtSimulatedSigma","Simulated E_{T} from #Xi^{-}");
870   CreateEtaPtHisto2D("EtSimulatedAntiSigma","Simulated E_{T} from #Xi^{+}");
871   CreateEtaPtHisto2D("EtSimulatedXi0","Simulated E_{T} from #Xi^{0}");
872   CreateEtaPtHisto2D("EtSimulatedAntiXi0","Simulated E_{T} from #Xi^{0}");
873   CreateEtaPtHisto2D("EtSimulatedAllHadron","Simulated E_{T} from all hadrons");
874
875
876   CreateEtaPtHisto2D("EtSimulatedLambdaDaughters","Simulated E_{T} from #Lambda Daughters");
877   CreateEtaPtHisto2D("EtSimulatedAntiLambdaDaughters","Simulated E_{T} from #bar{#Lambda} Daughters");
878   CreateEtaPtHisto2D("EtSimulatedK0SDaughters","Simulated E_{T} from K^{0}_{S} Daughters");
879   CreateEtaPtHisto2D("EtSimulatedOmegaDaughters","Simulated E_{T} from #Omega^{-} Daughters");
880   CreateEtaPtHisto2D("EtSimulatedAntiOmegaDaughters","Simulated E_{T} from #Omega^{+} Daughters");
881   CreateEtaPtHisto2D("EtSimulatedXiDaughters","Simulated E_{T} from #Xi^{-} Daughters");
882   CreateEtaPtHisto2D("EtSimulatedAntiXiDaughters","Simulated E_{T} from #Xi^{+} Daughters");
883
884
885   CreateEtaPtHisto2D("EtSimulatedGamma","Simulated E_{T} from #gamma");
886   CreateEtaPtHisto2D("EtSimulatedEta","Simulated E_{T} from #eta");
887   CreateEtaPtHisto2D("EtSimulatedPi0","Simulated E_{T} from #pi^{0}");
888   CreateEtaPtHisto2D("EtSimulatedOmega0","Simulated E_{T} from #omega");
889
890   TString *strTPC = new TString("TPC");
891   TString *strITS = new TString("ITS");
892   TString *strTPCITS = new TString("TPCITS");
893   for(Int_t i=0;i<3;i++){
894     TString *cutName;
895     Float_t maxPtdEdx = 10;
896     Float_t mindEdx = 35;
897     Float_t maxdEdx = 150.0;
898     switch(i){
899     case 0:
900       cutName = strTPC;
901       break;
902     case 1:
903       cutName = strITS;
904       maxPtdEdx = 5;
905       maxdEdx = 500.0;
906       break;
907     case 2:
908       cutName = strTPCITS;
909       break;
910     default:
911       cerr<<"Error:  cannot make histograms!"<<endl;
912       return;
913     }
914
915     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedPiPlus",cutName->Data()),"Reconstructed E_{T} from identified #pi^{+}");
916     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedPiMinus",cutName->Data()),"Reconstructed E_{T} from identified #pi^{-}");
917     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedKPlus",cutName->Data()),"Reconstructed E_{T} from identified K^{+}");
918     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedEMinus",cutName->Data()),"Reconstructed E_{T} from identified e^{-}");
919     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedEPlus",cutName->Data()),"Reconstructed E_{T} from identified e^{+}");
920     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedKMinus",cutName->Data()),"Reconstructed E_{T} from identified K^{-}");
921     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedProton",cutName->Data()),"Reconstructed E_{T} from identified p");
922     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedAntiProton",cutName->Data()),"Reconstructed E_{T} from identified #bar{p}");
923     CreateEtaPtHisto2D(Form("EtNReconstructed%sUnidentified",cutName->Data()),"Number of Reconstructed unidentified particles");
924     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedAssumingPion",cutName->Data()),"Reconstructed E_{T} from unidentified particles assuming pion mass");
925     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentified",cutName->Data()),"Reconstructed E_{T} from unidentified particles using real mass");
926     CreateEtaPtHisto2D(Form("EtReconstructed%sMisidentifiedElectrons",cutName->Data()),"Reconstructed E_{T} from misidentified electrons");
927
928
929     CreateEtaPtHisto2D(Form("EtReconstructed%sPiPlus",cutName->Data()),"Reconstructed E_{T} from #pi^{+}");
930     CreateEtaPtHisto2D(Form("EtReconstructed%sPiMinus",cutName->Data()),"Reconstructed E_{T} from #pi^{-}");
931     CreateEtaPtHisto2D(Form("EtReconstructed%sKPlus",cutName->Data()),"Reconstructed E_{T} from K^{+}");
932     CreateEtaPtHisto2D(Form("EtReconstructed%sKMinus",cutName->Data()),"Reconstructed E_{T} from K^{-}");
933     CreateEtaPtHisto2D(Form("EtReconstructed%sProton",cutName->Data()),"Reconstructed E_{T} from p");
934     CreateEtaPtHisto2D(Form("EtReconstructed%sAntiProton",cutName->Data()),"Reconstructed E_{T} from #bar{p}");
935     CreateEtaPtHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),"Reconstructed E_{T} from charged hadrons");
936     CreateEtaPtHisto2D(Form("EtNReconstructed%sPiPlus",cutName->Data()),"Reconstructed E_{T} from #pi^{+}");
937     CreateEtaPtHisto2D(Form("EtNReconstructed%sPiMinus",cutName->Data()),"Reconstructed E_{T} from #pi^{-}");
938     CreateEtaPtHisto2D(Form("EtNReconstructed%sKPlus",cutName->Data()),"Reconstructed E_{T} from K^{+}");
939     CreateEtaPtHisto2D(Form("EtNReconstructed%sKMinus",cutName->Data()),"Reconstructed E_{T} from K^{-}");
940     CreateEtaPtHisto2D(Form("EtNReconstructed%sProton",cutName->Data()),"Reconstructed E_{T} from p");
941     CreateEtaPtHisto2D(Form("EtNReconstructed%sAntiProton",cutName->Data()),"Reconstructed E_{T} from #bar{p}");
942     CreateEtaPtHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),"Reconstructed E_{T} from charged hadrons");
943
944     CreateEtaPtHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),"Reconstructed E_{T} from charged hadrons assuming they are all pions");
945     CreateEtaPtHisto2D(Form("EtReconstructed%sKPlusAssumingPion",cutName->Data()),"Reconstructed E_{T} from K^{+} assuming #pi mass");
946     CreateEtaPtHisto2D(Form("EtReconstructed%sKMinusAssumingPion",cutName->Data()),"Reconstructed E_{T} from K^{-} assuming #pi mass");
947     CreateEtaPtHisto2D(Form("EtReconstructed%sProtonAssumingPion",cutName->Data()),"Reconstructed E_{T} from p assuming #pi mass");
948     CreateEtaPtHisto2D(Form("EtReconstructed%sAntiProtonAssumingPion",cutName->Data()),"Reconstructed E_{T} from #bar{p} assuming #pi mass");
949
950     CreateEtaPtHisto2D(Form("EtReconstructed%sEPlus",cutName->Data()),"Reconstructed E_{T} from e^{+}");
951     CreateEtaPtHisto2D(Form("EtReconstructed%sEMinus",cutName->Data()),"Reconstructed E_{T} from e^{-}");
952
953
954
955     CreateEtaPtHisto2D(Form("EtReconstructed%sLambdaDaughters",cutName->Data()),"Reconstructed E_{T} from #Lambda Daughters");
956   CreateEtaPtHisto2D(Form("EtReconstructed%sAntiLambdaDaughters",cutName->Data()),"Reconstructed E_{T} from #bar{#Lambda} Daughters");
957   CreateEtaPtHisto2D(Form("EtReconstructed%sK0SDaughters",cutName->Data()),"Reconstructed E_{T} from K^{0}_{S} Daughters");
958   CreateEtaPtHisto2D(Form("EtReconstructed%sOmegaDaughters",cutName->Data()),"Reconstructed E_{T} from #Omega^{-} Daughters");
959   CreateEtaPtHisto2D(Form("EtReconstructed%sAntiOmegaDaughters",cutName->Data()),"Reconstructed E_{T} from #Omega^{+} Daughters");
960   CreateEtaPtHisto2D(Form("EtReconstructed%sXiDaughters",cutName->Data()),"Reconstructed E_{T} from #Xi^{-} Daughters");
961   CreateEtaPtHisto2D(Form("EtReconstructed%sAntiXiDaughters",cutName->Data()),"Reconstructed E_{T} from #Xi^{+} Daughters");
962
963     CreateIntHisto1D(Form("UnidentifiedPIDs%s",cutName->Data()),"PIDs of unidentified particles", "PID", "Number of particles",9, -4,4);
964     CreateHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),"PIDs of misidentified particles", "PID real","PID identified",5, -.5,4.5,5, -.5,4.5);
965     CreateHisto2D(Form("dEdxAll%s",cutName->Data()),"dE/dx for all particles","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
966     CreateHisto2D(Form("dEdxPion%s",cutName->Data()),"dE/dx for #pi^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
967     CreateHisto2D(Form("dEdxKaon%s",cutName->Data()),"dE/dx for K^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
968     CreateHisto2D(Form("dEdxProton%s",cutName->Data()),"dE/dx for p(#bar{p})","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
969     CreateHisto2D(Form("dEdxElectron%s",cutName->Data()),"dE/dx for e^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
970     CreateHisto2D(Form("dEdxUnidentified%s",cutName->Data()),"dE/dx for unidentified particles","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
971   }
972   CreateHisto2D("SimTotEtVsRecoTotEtFullAcceptanceTPC","Simulated total E_{T} vs reconstructed E_{T} with full acceptance for p_{T}>0.15 GeV/c","Simulated total E_{T}","Reconstructed total E_{T} (Full acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
973   CreateHisto2D("SimTotEtVsRecoTotEtEMCALAcceptanceTPC","Simulated total E_{T} vs reconstructed E_{T} with EMCAL acceptance for p_{T}>0.15 GeV/c","Simulated total E_{T}","Reconstructed total E_{T} (EMCAL acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
974   CreateHisto2D("SimTotEtVsRecoTotEtPHOSAcceptanceTPC","Simulated total E_{T} vs reconstructed E_{T} with PHOS acceptance for p_{T}>0.15 GeV/c","Simulated total E_{T}","Reconstructed total E_{T} (PHOS acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
975   CreateHisto2D("SimHadEtVsRecoHadEtFullAcceptanceTPC","Simulated hadronic E_{T} vs reconstructed E_{T} with full acceptance for p_{T}>0.15 GeV/c","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (Full acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
976   CreateHisto2D("SimHadEtVsRecoHadEtEMCALAcceptanceTPC","Simulated hadronic E_{T} vs reconstructed E_{T} with EMCAL acceptance for p_{T}>0.15 GeV/c","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (EMCAL acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
977   CreateHisto2D("SimHadEtVsRecoHadEtPHOSAcceptanceTPC","Simulated hadronic E_{T} vs reconstructed E_{T} with PHOS acceptance for p_{T}>0.15 GeV/c","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (PHOS acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
978   CreateHisto2D("SimTotEtVsRecoTotEtFullAcceptanceITS","Simulated total E_{T} vs reconstructed E_{T} with full acceptance for p_{T}>0.10 GeV/c","Simulated total E_{T}","Reconstructed total E_{T} (Full acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
979   CreateHisto2D("SimTotEtVsRecoTotEtEMCALAcceptanceITS","Simulated total E_{T} vs reconstructed E_{T} with EMCAL acceptance for p_{T}>0.10 GeV/c","Simulated total E_{T}","Reconstructed total E_{T} (EMCAL acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
980   CreateHisto2D("SimTotEtVsRecoTotEtPHOSAcceptanceITS","Simulated total E_{T} vs reconstructed E_{T} with PHOS acceptance for p_{T}>0.10 GeV/c","Simulated total E_{T}","Reconstructed total E_{T} (PHOS acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
981   CreateHisto2D("SimHadEtVsRecoHadEtFullAcceptanceITS","Simulated hadronic E_{T} vs reconstructed E_{T} with full acceptance for p_{T}>0.10 GeV/c","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (Full acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
982   CreateHisto2D("SimHadEtVsRecoHadEtEMCALAcceptanceITS","Simulated hadronic E_{T} vs reconstructed E_{T} with EMCAL acceptance for p_{T}>0.10 GeV/c","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (EMCAL acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
983   CreateHisto2D("SimHadEtVsRecoHadEtPHOSAcceptanceITS","Simulated hadronic E_{T} vs reconstructed E_{T} with PHOS acceptance for p_{T}>0.10 GeV/c","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (PHOS acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
984   CreateHisto2D("SimTotEtVsRecoTotEtFullAcceptanceTPCNoPID","Simulated total E_{T} vs reconstructed E_{T} with full acceptance for p_{T}>0.15 GeV/c, no PID","Simulated total E_{T}","Reconstructed total E_{T} (Full acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
985   CreateHisto2D("SimTotEtVsRecoTotEtEMCALAcceptanceTPCNoPID","Simulated total E_{T} vs reconstructed E_{T} with EMCAL acceptance for p_{T}>0.15 GeV/c, no PID","Simulated total E_{T}","Reconstructed total E_{T} (EMCAL acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
986   CreateHisto2D("SimTotEtVsRecoTotEtPHOSAcceptanceTPCNoPID","Simulated total E_{T} vs reconstructed E_{T} with PHOS acceptance for p_{T}>0.15 GeV/c, no PID","Simulated total E_{T}","Reconstructed total E_{T} (PHOS acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
987   CreateHisto2D("SimHadEtVsRecoHadEtFullAcceptanceTPCNoPID","Simulated hadronic E_{T} vs reconstructed E_{T} with full acceptance for p_{T}>0.15 GeV/c, no PID","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (Full acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
988   CreateHisto2D("SimHadEtVsRecoHadEtEMCALAcceptanceTPCNoPID","Simulated hadronic E_{T} vs reconstructed E_{T} with EMCAL acceptance for p_{T}>0.15 GeV/c, no PID","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (EMCAL acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
989   CreateHisto2D("SimHadEtVsRecoHadEtPHOSAcceptanceTPCNoPID","Simulated hadronic E_{T} vs reconstructed E_{T} with PHOS acceptance for p_{T}>0.15 GeV/c, no PID","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (PHOS acc., p_{T}>0.15 GeV/c)",100,0.0,100.0,100,0.0,100.0);
990   CreateHisto2D("SimTotEtVsRecoTotEtFullAcceptanceITSNoPID","Simulated total E_{T} vs reconstructed E_{T} with full acceptance for p_{T}>0.10 GeV/c, no PID","Simulated total E_{T}","Reconstructed total E_{T} (Full acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
991   CreateHisto2D("SimTotEtVsRecoTotEtEMCALAcceptanceITSNoPID","Simulated total E_{T} vs reconstructed E_{T} with EMCAL acceptance for p_{T}>0.10 GeV/c, no PID","Simulated total E_{T}","Reconstructed total E_{T} (EMCAL acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
992   CreateHisto2D("SimTotEtVsRecoTotEtPHOSAcceptanceITSNoPID","Simulated total E_{T} vs reconstructed E_{T} with PHOS acceptance for p_{T}>0.10 GeV/c, no PID","Simulated total E_{T}","Reconstructed total E_{T} (PHOS acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
993   CreateHisto2D("SimHadEtVsRecoHadEtFullAcceptanceITSNoPID","Simulated hadronic E_{T} vs reconstructed E_{T} with full acceptance for p_{T}>0.10 GeV/c, no PID","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (Full acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
994   CreateHisto2D("SimHadEtVsRecoHadEtEMCALAcceptanceITSNoPID","Simulated hadronic E_{T} vs reconstructed E_{T} with EMCAL acceptance for p_{T}>0.10 GeV/c, no PID","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (EMCAL acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
995   CreateHisto2D("SimHadEtVsRecoHadEtPHOSAcceptanceITSNoPID","Simulated hadronic E_{T} vs reconstructed E_{T} with PHOS acceptance for p_{T}>0.10 GeV/c, no PID","Simulated hadronic E_{T}","Reconstructed hadronic E_{T} (PHOS acc., p_{T}>0.10 GeV/c)",100,0.0,100.0,100,0.0,100.0);
996   CreateIntHisto1D("NEvents","Number of events","number of events","Number of events",1,0,1);
997
998   //CreateHisto1D("MisidentifiedPIDs","PIDs for particles misidentified that are not a #pi,K,p","PID","number of entries",3000,0.5,3000.5);
999 }
1000