6e3e26e7da6d846a4e8811cde6f541306100749d
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraAllChAOD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //-----------------------------------------------------------------
17 //         AliAnalysisTaskSpectraAllChAOD class
18 //-----------------------------------------------------------------
19
20 #include "TChain.h"
21 #include "TTree.h"
22 #include "TLegend.h"
23 #include "TH1F.h"
24 #include "TH2F.h"
25 #include "THnSparse.h"
26 #include "TCanvas.h"
27 #include "AliAnalysisTask.h"
28 #include "AliAODTrack.h"
29 #include "AliAODMCParticle.h"
30 #include "AliVParticle.h"
31 #include "AliAODEvent.h"
32 #include "AliAODInputHandler.h"
33 #include "AliAnalysisTaskSpectraAllChAOD.h"
34 #include "AliAnalysisTaskESDfilter.h"
35 #include "AliAnalysisDataContainer.h"
36 #include "AliSpectraAODTrackCuts.h"
37 #include "AliSpectraAODEventCuts.h"
38 #include "AliHelperPID.h"
39 #include "AliPIDCombined.h"
40 #include "AliCentrality.h"
41 #include "TProof.h"
42 #include "AliVEvent.h"
43 #include "AliStack.h"
44 #include <TMCProcess.h>
45
46 #include <iostream>
47
48 using namespace AliHelperPIDNameSpace;
49 using namespace std;
50
51 ClassImp(AliAnalysisTaskSpectraAllChAOD)
52
53 //________________________________________________________________________
54 AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),  
55   fAOD(0x0),
56   fTrackCuts(0x0),
57   fEventCuts(0x0),
58   fHelperPID(0x0),
59   fIsMC(0),
60   fDoDoubleCounting(0),
61   fFillOnlyEvents(0),
62   fCharge(0),
63   fVZEROside(0),
64   fOutput(0x0),
65   fnCentBins(20),
66   fnQvecBins(100),
67   fnNchBins(200),
68   fIsQvecCalibMode(0),
69   fQvecUpperLim(100),
70   fIsAOD160(1),
71   fnDCABins(60),
72   fDCAmin(-3),
73   fDCAmax(3),
74   fDCAzCut(999999.)
75 {
76   // Default constructor
77   DefineInput(0, TChain::Class());
78   DefineOutput(1, TList::Class());
79   DefineOutput(2, AliSpectraAODEventCuts::Class());
80   DefineOutput(3, AliSpectraAODTrackCuts::Class());
81   DefineOutput(4, AliHelperPID::Class());
82 }
83
84 //________________________________________________________________________
85 void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
86 {
87   // create output objects
88   fOutput=new TList();
89   fOutput->SetOwner();
90   fOutput->SetName("fOutput");
91   
92   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
93   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
94   if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");
95   
96   // binning common to all the THn
97   const Double_t ptBins[] = {0.20,0.30,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.,20.,25.,30.,35.,40.,50.,75.,100.};
98   const Int_t nptBins=34;
99   //const Double_t ptBins[] = {0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.,20.,25.,30.,35.,40.,50.,75.,100.};
100   //const Int_t nptBins=44;
101   
102   //dimensions of THnSparse for tracks
103   const Int_t nvartrk=8;
104   //                                             pt          cent          Q vec       IDrec          IDgen        isph         y             DCA 
105   Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,     fnQvecBins,        4,             3,          3,        2,        fnDCABins};
106   Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,              0.,       -.5,          -0.5,        0.5,      -0.5,          fDCAmin};
107   Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,     fQvecUpperLim,       3.5,      2.5,        3.5,      0.5,         fDCAmax};    
108   THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
109   NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
110   NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
111   NSparseHistTrk->SetBinEdges(0,ptBins);
112   NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
113   NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
114   NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
115   NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
116   NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
117   NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
118   NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
119   NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
120   NSparseHistTrk->GetAxis(5)->SetTitle("isph");
121   NSparseHistTrk->GetAxis(5)->SetName("isph");
122   NSparseHistTrk->GetAxis(6)->SetTitle("y");
123   NSparseHistTrk->GetAxis(6)->SetName("y");
124   NSparseHistTrk->GetAxis(7)->SetTitle("dca");
125   NSparseHistTrk->GetAxis(7)->SetName("dca");
126   fOutput->Add(NSparseHistTrk);
127   
128   //dimensions of THnSparse for stack
129   const Int_t nvarst=6;
130   //                                             pt          cent    IDgen        isph        y    etaselected
131   Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2,             1};
132   Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5,            0.5};
133   Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5,           1.5};
134   THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
135   NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
136   NSparseHistSt->SetBinEdges(0,ptBins);
137   NSparseHistSt->GetAxis(0)->SetName("pT_rec");
138   NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
139   NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
140   NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
141   NSparseHistSt->GetAxis(2)->SetName("ID_gen");
142   NSparseHistSt->GetAxis(3)->SetTitle("isph");
143   NSparseHistSt->GetAxis(3)->SetName("isph");
144   NSparseHistSt->GetAxis(4)->SetTitle("y");
145   NSparseHistSt->GetAxis(4)->SetName("y");
146   NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
147   NSparseHistSt->GetAxis(5)->SetName("etaselected");
148   fOutput->Add(NSparseHistSt);
149   
150   //dimensions of THnSparse for the normalization
151   const Int_t nvarev=3;
152   //                                             cent             Q vec                Nch
153   Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           fnNchBins};
154   Double_t xminHistRealEv[nvarev] = {           0.,               0.,                   0.};
155   Double_t xmaxHistRealEv[nvarev] = {       100.,  fQvecUpperLim,               2000.};
156   THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
157   NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
158   NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
159   NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
160   NSparseHistEv->GetAxis(1)->SetName("Q_vec");
161   NSparseHistEv->GetAxis(2)->SetTitle("N charged");
162   NSparseHistEv->GetAxis(2)->SetName("N_ch");
163   fOutput->Add(NSparseHistEv);
164   
165   PostData(1, fOutput  );
166   PostData(2, fEventCuts);
167   PostData(3, fTrackCuts);
168   PostData(4, fHelperPID);
169 }
170
171 //________________________________________________________________________
172
173 void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
174 {
175   //Printf("An event");
176   // main event loop
177   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
178   if (!fAOD) {
179     AliWarning("ERROR: AliAODEvent not available \n");
180     return;
181   }
182   
183   if (strcmp(fAOD->ClassName(), "AliAODEvent"))
184     {
185       AliFatal("Not processing AODs");
186     }
187   
188   if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
189
190   //Default TPC priors
191   //if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
192   
193   Double_t Qvec=0.;//in case of MC we save space in the memory
194   if(!fIsMC){
195     if(fIsQvecCalibMode){
196       if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
197       else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
198     }
199     else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
200   }
201   
202   Double_t Cent=fEventCuts->GetCent();
203   
204   // First do MC to fill up the MC particle array
205   TClonesArray *arrayMC = 0;
206   if (fIsMC)
207     {
208       arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
209       if (!arrayMC) {
210         AliFatal("Error: MC particles branch not found!\n");
211       }
212       Int_t nMC = arrayMC->GetEntries();
213       for (Int_t iMC = 0; iMC < nMC; iMC++)
214         {
215           AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
216           if(!partMC->Charge()) continue;//Skip neutrals
217           if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
218           
219           //flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
220           Double_t etaselected=-1.;
221           if(partMC->Eta()>=fTrackCuts->GetEtaMin() && partMC->Eta()<=fTrackCuts->GetEtaMax())etaselected=1.;
222           
223           //pt     cent    IDgen        isph        y
224           Double_t varSt[6];
225           varSt[0]=partMC->Pt();
226           varSt[1]=Cent;
227           varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
228           varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
229           varSt[4]=partMC->Y();
230           varSt[5]=etaselected;
231           ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
232           
233           //Printf("a particle");
234           
235         }
236     }
237   
238   //main loop on tracks
239   
240   Int_t Nch = 0.;
241   
242   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
243     AliAODTrack* track = fAOD->GetTrack(iTracks);
244     if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge 
245     if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
246     if(!fFillOnlyEvents){
247       Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
248       Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
249       Int_t IDgen=kSpUndefined;//set if MC
250       Int_t isph=-999;
251 //       Int_t iswd=-999;
252       
253       if (arrayMC) {
254         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
255         if (!partMC) { 
256           AliError("Cannot get MC particle");
257           continue; 
258         }
259         IDgen=fHelperPID->GetParticleSpecies(partMC);
260         isph=partMC->IsPhysicalPrimary();
261         //iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014
262         
263         if(fIsAOD160){// enabled for new ADO160 only
264           if(partMC->IsSecondaryFromWeakDecay()) isph=2.;
265           if(partMC->IsSecondaryFromMaterial()) isph=3.;
266
267         }
268       }
269       
270       /*** DCA ***/
271       Double_t dcaxy = -999.;
272       Double_t dcaz = -999.;
273       
274       Double_t p[2]; 
275       if(GetDCA(track,p)){ dcaxy=p[0]; dcaz=p[1];}
276       
277       if(dcaz >= fDCAzCut) continue;
278       
279     //pt     cent    Q vec     IDrec     IDgen       isph      y
280       Double_t varTrk[8];
281       varTrk[0]=track->Pt();
282       varTrk[1]=Cent;
283       varTrk[2]=Qvec;
284       varTrk[3]=(Double_t)IDrec;
285       varTrk[4]=(Double_t)IDgen;
286       varTrk[5]=(Double_t)isph;
287       varTrk[6]=y;
288       varTrk[7]=dcaxy;
289       ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
290       
291       //for nsigma PID fill double counting of ID
292       if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
293         Bool_t *HasDC;
294         HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
295         for(Int_t ipart=0;ipart<kNSpecies;ipart++){
296           if(HasDC[ipart]==kTRUE){
297             varTrk[3]=(Double_t)ipart;
298             ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
299           }
300         }
301       }
302       
303       //fill all charged (3)
304       varTrk[3]=3.;
305       varTrk[4]=3.;
306       ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
307     }//end if fFillOnlyEvents
308     
309     //Printf("a track");
310       
311     Nch++;
312   } // end loop on tracks
313   
314   Double_t varEv[3];
315   varEv[0]=Cent;
316   varEv[1]=Qvec;
317   varEv[2]=Nch;
318   ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
319   
320   PostData(1, fOutput  );
321   PostData(2, fEventCuts);
322   PostData(3, fTrackCuts);
323   PostData(4, fHelperPID);
324 }
325
326 //_________________________________________________________________
327 Bool_t  AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
328   
329   //AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
330   //FIXME should update EventCuts?
331   //FIXME add track->GetXYZ(p) method
332   
333   double xyz[3],cov[3];
334   
335   if (!trk->GetXYZ(xyz)) { // dca is not stored
336     AliExternalTrackParam etp;
337     etp.CopyFromVTrack(trk);
338     AliVEvent* ev = (AliVEvent*)trk->GetEvent();
339     if (!ev) {/*printf("Event is not connected to the track\n");*/ return kFALSE;}
340     if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
341   }
342   p[0] = xyz[0];
343   p[1] = xyz[1];
344   return kTRUE;
345
346 }
347
348 //_________________________________________________________________
349 void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
350 {
351   // Terminate
352 }