Taking into account that only 1 or 2 values may be present for the
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonFeedDownAnalysis.cxx
1 #include <Riostream.h>
2 #include <TFile.h>
3 #include <TSystem.h>
4 #include <TF1.h>
5 #include <TH2D.h>
6 #include <TH1D.h>
7 #include <TH1I.h>
8 #include <TParticle.h>
9
10 #include <AliExternalTrackParam.h>
11 #include <AliAODEvent.h>
12 #include <AliESDEvent.h>
13 //#include <AliLog.h>
14 #include <AliPID.h>
15 #include <AliStack.h>
16 #include <AliCFContainer.h>
17 #include <AliCFEffGrid.h>
18 #include <AliCFDataGrid.h>
19 //#include <AliESDVertex.h>
20 class AliLog;
21 class AliESDVertex;
22
23 #include "AliProtonFeedDownAnalysis.h"
24 #include "AliProtonAnalysisBase.h"
25
26 ClassImp(AliProtonFeedDownAnalysis)
27
28 //____________________________________________________________________//
29 AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis() : 
30   TObject(), fProtonAnalysisBase(0),
31   fNBinsY(0), fMinY(0), fMaxY(0),
32   fNBinsPt(0), fMinPt(0), fMaxPt(0),
33   fProtonContainer(0), fAntiProtonContainer(0), fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)
34   {
35   //Default constructor
36  }
37 //____________________________________________________________________//
38 AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) : 
39   TObject(), fProtonAnalysisBase(0),
40   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
41   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
42   fProtonContainer(0), fAntiProtonContainer(0),fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)
43   {
44         //Default constructor
45         
46         //setting up the containers
47         Int_t iBin[2];
48         iBin[0] = nbinsY;
49         iBin[1] = nbinsPt;
50         Double_t *binLimY = new Double_t[nbinsY+1];
51         Double_t *binLimPt = new Double_t[nbinsPt+1];
52         //values for bin lower bounds
53         for(Int_t i = 0; i <= nbinsY; i++) 
54                 binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
55         for(Int_t i = 0; i <= nbinsPt; i++) 
56                 binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
57         
58         fProtonContainer = new AliCFContainer("containerProtons","container for protons",4,2,iBin);
59         fProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
60         fProtonContainer->SetBinLimits(1,binLimPt); //pT
61         fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",4,2,iBin);
62         fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
63         fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
64         
65 }
66 //____________________________________________________________________//
67 AliProtonFeedDownAnalysis::~AliProtonFeedDownAnalysis() 
68 {
69         //Default destructor
70         if(fProtonAnalysisBase) delete fProtonAnalysisBase;
71         if(fProtonContainer) delete fProtonContainer;
72         if(fAntiProtonContainer) delete fAntiProtonContainer;
73         if(fweightfunction) delete fweightfunction; 
74         if(fLambda) delete fLambda;
75         if(fLambdaweighted) delete fLambdaweighted;
76         if(fAntiLambda) delete fAntiLambda;
77         if(fAntiLambdaweighted) delete fAntiLambdaweighted;
78 }
79 //____________________________________________________________________//
80 void AliProtonFeedDownAnalysis::InitAnalysisHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) 
81 {
82         //Initializes the histograms
83         fNBinsY = nbinsY;
84         fMinY = fLowY;
85         fMaxY = fHighY;
86         fNBinsPt = nbinsPt;
87         fMinPt = fLowPt;
88         fMaxPt = fHighPt;
89         
90         
91         //setting up the containers
92         Int_t iBin[2];
93         iBin[0] = nbinsY;
94         iBin[1] = nbinsPt;
95         Double_t *binLimY = new Double_t[nbinsY+1];
96         Double_t *binLimPt = new Double_t[nbinsPt+1];
97         //values for bin lower bounds
98         for(Int_t i = 0; i <= nbinsY; i++) 
99                 binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
100         for(Int_t i = 0; i <= nbinsPt; i++) 
101                 binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
102         
103         fProtonContainer = new AliCFContainer("containerProtons","container for protons",kNSteps,2,iBin);
104         fProtonContainer->SetBinLimits(0,binLimY); //rapidity
105         fProtonContainer->SetBinLimits(1,binLimPt); //pT
106         fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",kNSteps,2,iBin);
107         fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
108         fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
109         fLambda=new TH1F("Lambda","Lambda",35,0.5,4.0);
110         fLambdaweighted=new TH1F("Lambdaweighted","Lambdaweighted",35,0.5,4.0);
111         fAntiLambda=new TH1F("AntiLambda","AntiLambda",35,0.5,4.0);
112         fAntiLambdaweighted=new TH1F("AntiLambdaweighted","AntiLambdaweighted",35,0.5,4.0);
113 }
114 //_________________________________________________________________________//
115 void AliProtonFeedDownAnalysis::Analyze(AliESDEvent *esd, const AliESDVertex *vertex,AliStack *stack)
116 {       
117         Int_t nTracks = 0;
118         Int_t nIdentifiedProtons = 0, nIdentifiedAntiProtons = 0;
119         Int_t nSurvivedProtons = 0, nSurvivedAntiProtons = 0;
120         
121         //fHistEvents->Fill(0); //number of analyzed events
122         Double_t containerInput[2] ;
123         Double_t gPt = 0.0, gP = 0.0;
124         nTracks = esd->GetNumberOfTracks();
125         Float_t weight;
126         Int_t trackflag;
127         for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) 
128         {
129                 
130                 AliESDtrack* track = esd->GetTrack(iTracks);
131                 AliESDtrack trackTPC;
132                 Int_t label= track->GetLabel();
133         /*      Int_t trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else 
134                 if (trackflag!=0)
135                         weight=GetWeightforProton(label,stack); 
136                 else
137                         weight=1.0;     */
138         
139                 if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid)) 
140                 {
141                         AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
142                         if(!tpcTrack)
143                                  continue;
144                         gPt = tpcTrack->Pt();
145                         gP = tpcTrack->P();
146                         if(fProtonAnalysisBase->IsProton(track)) 
147                         {
148                                 trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else 
149                                 if (trackflag!=0)
150                                         weight=GetWeightforProton(label,stack); 
151                                 else
152                                         weight=1.0;     
153                                 if(tpcTrack->Charge() > 0) 
154                                 {
155                                         nIdentifiedProtons += 1;
156                                         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) 
157                                                 continue;//track cuts
158                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) 
159                                                 continue; //track outside the analyzed y-Pt
160                                         nSurvivedProtons += 1;
161                                         if(fProtonAnalysisBase->GetEtaMode()) 
162                                         {
163                                                 containerInput[0] = tpcTrack->Eta();
164                                         }
165                                         else 
166                                         {
167                                                 //fill the container
168                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
169                                         }
170                                         containerInput[1] = gPt;
171                                         fProtonContainer->Fill(containerInput,kSelected,weight);   
172                                         //Feed-down check
173                                         if (trackflag==1)
174                                                 fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight); 
175                                         
176                                 }//protons
177                                 else if(tpcTrack->Charge() < 0) 
178                                 {
179                                         nIdentifiedAntiProtons += 1;
180                                         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) 
181                                                 continue;//track cuts
182                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) 
183                                                 continue; //track outside the analyzed y-Pt
184                                         nSurvivedAntiProtons += 1;
185                                         if(fProtonAnalysisBase->GetEtaMode()) 
186                                         {
187                                                 containerInput[0] = tpcTrack->Eta();
188                                         }
189                                         else 
190                                         {
191                                                 //fill the container
192                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
193                                         }
194                                         containerInput[1] = gPt;
195                                         fAntiProtonContainer->Fill(containerInput,kSelected,weight);
196                                         //Feed-down check                                                       
197                                         if(trackflag==-1)
198                                                 fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);                                  
199                                 }//antiprotons   
200                         }//proton check
201                 }//TPC only tracks
202                 else if(fProtonAnalysisBase->GetAnalysisMode() == AliProtonAnalysisBase::kGlobal) 
203                 {
204                         gPt = track->Pt();
205                         gP = track->P();
206                         if(fProtonAnalysisBase->IsProton(track)) 
207                         {
208                                 trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else 
209                                 if (trackflag!=0)
210                                         weight=GetWeightforProton(label,stack); 
211                                 else
212                                         weight=1.0;     
213                                 if(track->Charge() > 0) 
214                                 {
215                                         nIdentifiedProtons += 1;
216                                         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) 
217                                                 continue;//track cuts
218                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) 
219                                                 continue; //track outside the analyzed y-Pt
220                                         nSurvivedProtons += 1;
221                                         if(fProtonAnalysisBase->GetEtaMode()) 
222                                         {
223                                                 containerInput[0] = track->Eta();
224                                         }
225                                         else 
226                                         {
227                                                 //fill the container
228                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
229                                         }
230                                         containerInput[1] = gPt;
231                                         fProtonContainer->Fill(containerInput,kSelected,weight);  
232                                          //Feed-down check
233                                         if (trackflag==1)
234                                                 fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight); 
235                                 }//protons
236                                 else if(track->Charge() < 0) 
237                                 {
238                                         nIdentifiedAntiProtons += 1;
239                                         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) 
240                                                 continue;//track cuts
241                                         if(!fProtonAnalysisBase->IsInPhaseSpace(track))
242                                                 continue; //track outside the analyzed y-Pt
243                                         nSurvivedAntiProtons += 1;
244                                         if(fProtonAnalysisBase->GetEtaMode()) 
245                                         {
246                                                 containerInput[0] = track->Eta();
247                                         }
248                                         else 
249                                         {
250                                                 //fill the container
251                                                 containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
252                                         }
253                                         containerInput[1] = gPt;
254                                         fAntiProtonContainer->Fill(containerInput,kSelected,weight);   
255                                         if(trackflag==-1)
256                                                 fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
257                                 }//antiprotons
258                         }//proton check 
259                 }//combined tracking
260         }//track loop 
261         
262         if(fProtonAnalysisBase->GetDebugMode())
263                 Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
264
265 }
266 //_________________________________________________________________________//
267 void AliProtonFeedDownAnalysis::Analyze(AliAODEvent *fAOD)
268 {
269   // Analysis from AOD: to be implemented...
270   // in the near future.
271   fAOD->Print();
272
273 }
274 //_________________________________________________________________________//
275 void AliProtonFeedDownAnalysis::Analyze(AliStack *stack)
276 {
277          Double_t containerInput[2] ;
278          Float_t weight;
279         for (Int_t ipart=0; ipart<stack->GetNtrack(); ipart++) 
280         { 
281                 TParticle *particle  = stack->Particle(ipart);
282                 Int_t code=particle->GetPdgCode();
283                 if (code==3122)
284                 {
285                         fLambda->Fill(particle->Pt());
286                         fLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));                       
287                 }
288                 if (code==-3122)
289                 {
290                         fAntiLambda->Fill(particle->Pt());
291                         fAntiLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));                   
292                 }
293                 if (TMath::Abs(code)==2212)
294                 {
295                         Int_t trackflag=LambdaIsMother(ipart,stack);//1 mother lambda -1 mother anti lambda 0 mother something else 
296                         if (trackflag!=0)
297                                 weight=GetWeightforProton(ipart,stack);
298                         else
299                                 weight=1.0;     
300                         if(particle->GetUniqueID() == 13) //recject hadronic inter.
301                                 continue; 
302                         if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt)) 
303                                 continue;
304                         if(fProtonAnalysisBase->GetEtaMode()) 
305                         {
306                                 if((particle->Eta()> fMaxY)||(particle->Eta()< fMinY))
307                                         continue; 
308                         }       
309                         else
310                         {       
311                                 if((fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz())< fMinY)) 
312                                         continue;
313                         }               
314                         if(fProtonAnalysisBase->GetEtaMode()) 
315                         {
316                                 containerInput[0] =particle->Eta();
317                         }
318                         else 
319                         {
320                                 containerInput[0] = fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz());
321                         }
322                         //containerInput[0] = particle->Eta() ;
323                         containerInput[1] = (Float_t)particle->Pt();
324                         if (particle->GetPdgCode()==2212)
325                         {
326                                 fProtonContainer->Fill(containerInput,kAll,weight);
327                                 if(ipart<stack->GetNprimary())
328                                 {
329                                         fProtonContainer->Fill(containerInput,kPrimary,weight);
330                                 }
331                                 else
332                                 {
333                                         if(trackflag==1)
334                                                 fProtonContainer->Fill(containerInput,kFromLambda,weight);
335                                 }
336                         }
337                         if (particle->GetPdgCode()==-2212)
338                         {
339                                 fAntiProtonContainer->Fill(containerInput,kAll,weight);
340                                 if(ipart<stack->GetNprimary())
341                                 {
342                                         fAntiProtonContainer->Fill(containerInput,kPrimary,weight);
343                                 }
344                                 else
345                                 {                                                               
346                                         if(trackflag==-1)
347                                                 fAntiProtonContainer->Fill(containerInput,kFromLambda,weight);
348                                 }                       
349                         }
350                 }
351         }   
352 }
353 //______________________________________________________________________________________________
354 Int_t AliProtonFeedDownAnalysis::LambdaIsMother(Int_t number, AliStack *stack)
355 {
356          if(number<0)
357                 return 0;
358         TParticle *particle  = stack->Particle(number);
359         Int_t motherPDG=0;
360         Int_t lmother=-1;
361         lmother=particle->GetFirstMother();
362         if (lmother<0)          
363                 return 0;
364         TParticle *mparticle=stack->Particle(lmother);
365         motherPDG=mparticle->GetPdgCode();                                                                              
366         if(motherPDG==3122)
367                 return 1;
368         if(motherPDG==-3122)    
369                 return -1;
370         return 0;       
371
372 //___________________________________________________________________________________________
373 Float_t AliProtonFeedDownAnalysis::GetWeightforProton(Int_t number,AliStack *stack)
374 {
375          if(number<0)
376                 return 1.0;
377         TParticle *particle  = stack->Particle(number);
378         Int_t lmother=particle->GetFirstMother();
379         TParticle *mparticle=stack->Particle(lmother);
380         return  fweightfunction->Eval(mparticle->Pt());
381 }
382 //______________________________________________________________________________________________
383 Float_t AliProtonFeedDownAnalysis::GetWeightforLambda(Float_t pt)
384 {
385         return  fweightfunction->Eval(pt);
386 }