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