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