1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-----------------------------------------------------------------
17 // AliAnalysisTaskSpectraAllChAOD class
18 //-----------------------------------------------------------------
25 #include "THnSparse.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"
42 #include "AliVEvent.h"
44 #include <TMCProcess.h>
48 using namespace AliHelperPIDNameSpace;
51 ClassImp(AliAnalysisTaskSpectraAllChAOD)
53 //________________________________________________________________________
54 AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),
77 fDoCentrSystCentrality(0)
79 // Default constructor
80 DefineInput(0, TChain::Class());
81 DefineOutput(1, TList::Class());
82 DefineOutput(2, AliSpectraAODEventCuts::Class());
83 DefineOutput(3, AliSpectraAODTrackCuts::Class());
84 DefineOutput(4, AliHelperPID::Class());
87 //________________________________________________________________________
88 void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
90 // create output objects
93 fOutput->SetName("fOutput");
95 if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
96 if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
97 if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");
99 // binning common to all the THn
100 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.};
101 const Int_t nptBins=34;
102 //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.};
103 //const Int_t nptBins=44;
105 //dimensions of THnSparse for tracks
106 const Int_t nvartrk=8;
107 // pt cent Q vec IDrec IDgen isph y DCA
108 Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 3, 2, fnDCABins};
109 Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, 0.5, -0.5, fDCAmin};
110 Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., fQvecUpperLim, 3.5, 2.5, 3.5, 0.5, fDCAmax};
111 THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
112 NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
113 NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
114 NSparseHistTrk->SetBinEdges(0,ptBins);
115 NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
116 NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
117 NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
118 NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
119 NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
120 NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
121 NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
122 NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
123 NSparseHistTrk->GetAxis(5)->SetTitle("isph");
124 NSparseHistTrk->GetAxis(5)->SetName("isph");
125 NSparseHistTrk->GetAxis(6)->SetTitle("y");
126 NSparseHistTrk->GetAxis(6)->SetName("y");
127 NSparseHistTrk->GetAxis(7)->SetTitle("dca");
128 NSparseHistTrk->GetAxis(7)->SetName("dca");
129 fOutput->Add(NSparseHistTrk);
131 //dimensions of THnSparse for stack
132 const Int_t nvarst=7;
133 // pt cent IDgen isph y etaselected Qvec gen
134 Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2, 1, fnQvecBins};
135 Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5, 0.5, 0.};
136 Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5, 1.5, fQvecUpperLim};
137 THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
138 NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
139 NSparseHistSt->SetBinEdges(0,ptBins);
140 NSparseHistSt->GetAxis(0)->SetName("pT_rec");
141 NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
142 NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
143 NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
144 NSparseHistSt->GetAxis(2)->SetName("ID_gen");
145 NSparseHistSt->GetAxis(3)->SetTitle("isph");
146 NSparseHistSt->GetAxis(3)->SetName("isph");
147 NSparseHistSt->GetAxis(4)->SetTitle("y");
148 NSparseHistSt->GetAxis(4)->SetName("y");
149 NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
150 NSparseHistSt->GetAxis(5)->SetName("etaselected");
151 NSparseHistSt->GetAxis(6)->SetTitle("Q vec gen");
152 NSparseHistSt->GetAxis(6)->SetName("Q_vec_gen");
153 fOutput->Add(NSparseHistSt);
155 //dimensions of THnSparse for the normalization
156 const Int_t nvarev=4;
157 // cent Q vec Nch Qvec gen
158 Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, fnNchBins, fnQvecBins};
159 Double_t xminHistRealEv[nvarev] = { 0., 0., 0., 0.};
160 Double_t xmaxHistRealEv[nvarev] = { 100., fQvecUpperLim, 2000., fQvecUpperLim};
161 THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
162 NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
163 NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
164 NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
165 NSparseHistEv->GetAxis(1)->SetName("Q_vec");
166 NSparseHistEv->GetAxis(2)->SetTitle("N charged");
167 NSparseHistEv->GetAxis(2)->SetName("N_ch");
168 NSparseHistEv->GetAxis(3)->SetTitle("Q vec gen");
169 NSparseHistEv->GetAxis(3)->SetName("Q_vec_gen");
170 fOutput->Add(NSparseHistEv);
172 PostData(1, fOutput );
173 PostData(2, fEventCuts);
174 PostData(3, fTrackCuts);
175 PostData(4, fHelperPID);
178 //________________________________________________________________________
180 void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
182 //Printf("An event");
184 fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
186 AliWarning("ERROR: AliAODEvent not available \n");
190 if (strcmp(fAOD->ClassName(), "AliAODEvent"))
192 AliFatal("Not processing AODs");
195 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
198 if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME we should modify the task to change priors
201 if(fIsQvecCalibMode){
202 if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
203 else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
204 else if (fVZEROside==2)Qvec=fEventCuts->GetqTPC();
206 else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
208 Double_t QvecMC = 0.;
210 if(fIsQvecCalibMode){
211 QvecMC = fEventCuts->CalculateQVectorMC(fVZEROside, fQgenType);
213 else QvecMC = fEventCuts->GetQvecPercentileMC(fVZEROside, fQgenType);
216 Double_t Cent=(fDoCentrSystCentrality)?1.01*fEventCuts->GetCent():fEventCuts->GetCent();
218 // First do MC to fill up the MC particle array
219 TClonesArray *arrayMC = 0;
222 arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
224 AliFatal("Error: MC particles branch not found!\n");
226 Int_t nMC = arrayMC->GetEntries();
227 for (Int_t iMC = 0; iMC < nMC; iMC++)
229 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
230 if(!partMC->Charge()) continue;//Skip neutrals
231 if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
233 //flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
234 Double_t etaselected=-1.;
235 if(partMC->Eta()>=fTrackCuts->GetEtaMin() && partMC->Eta()<=fTrackCuts->GetEtaMax())etaselected=1.;
237 //pt cent IDgen isph y
239 varSt[0]=partMC->Pt();
241 varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
242 varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
243 varSt[4]=partMC->Y();
244 varSt[5]=etaselected;
246 ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
248 //Printf("a particle");
253 //main loop on tracks
257 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
258 AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
259 if(!track) AliFatal("Not a standard AOD");
260 if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
261 if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
262 if(!fFillOnlyEvents){
263 Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
264 Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
265 Int_t IDgen=kSpUndefined;//set if MC
270 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
272 AliError("Cannot get MC particle");
275 IDgen=fHelperPID->GetParticleSpecies(partMC);
276 isph=partMC->IsPhysicalPrimary();
277 //iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014
279 if(fIsAOD160){// enabled for new ADO160 only
280 if(partMC->IsSecondaryFromWeakDecay()) isph=2.;
281 if(partMC->IsSecondaryFromMaterial()) isph=3.;
287 Double_t dcaxy = -999.;
288 Double_t dcaz = -999.;
291 if(GetDCA(track,p)){ dcaxy=p[0]; dcaz=p[1];}
293 if(dcaz >= fDCAzCut) continue;
295 //pt cent Q vec IDrec IDgen isph y
297 varTrk[0]=track->Pt();
299 if(fIsMC && fQvecGen) varTrk[2]=QvecMC;
301 varTrk[3]=(Double_t)IDrec;
302 varTrk[4]=(Double_t)IDgen;
303 varTrk[5]=(Double_t)isph;
306 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
308 //for nsigma PID fill double counting of ID
309 if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
311 HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
312 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
313 if(HasDC[ipart]==kTRUE){
314 varTrk[3]=(Double_t)ipart;
315 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
320 //fill all charged (3)
323 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
324 }//end if fFillOnlyEvents
329 } // end loop on tracks
336 ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
338 PostData(1, fOutput );
339 PostData(2, fEventCuts);
340 PostData(3, fTrackCuts);
341 PostData(4, fHelperPID);
344 //_________________________________________________________________
345 Bool_t AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
347 //AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
348 //FIXME should update EventCuts?
349 //FIXME add track->GetXYZ(p) method
351 double xyz[3],cov[3];
353 if (!trk->GetXYZ(xyz)) { // dca is not stored
354 AliExternalTrackParam etp;
355 etp.CopyFromVTrack(trk);
356 AliVEvent* ev = (AliVEvent*)trk->GetEvent();
357 if (!ev) {/*printf("Event is not connected to the track\n");*/ return kFALSE;}
358 if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
366 //_________________________________________________________________
367 void AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)