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),
69 // Default constructor
70 DefineInput(0, TChain::Class());
71 DefineOutput(1, TList::Class());
72 DefineOutput(2, AliSpectraAODEventCuts::Class());
73 DefineOutput(3, AliSpectraAODTrackCuts::Class());
74 DefineOutput(4, AliHelperPID::Class());
77 //________________________________________________________________________
78 void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
80 // create output objects
83 fOutput->SetName("fOutput");
85 if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
86 if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
87 if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");
89 // binning common to all the THn
90 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.};
91 const Int_t nptBins=26;
93 //dimensions of THnSparse for tracks
94 const Int_t nvartrk=8;
95 // pt cent Q vec IDrec IDgen isph iswd y
96 Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 2, 2, 2};
97 Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, -0.5, -0.5, -0.5};
98 Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., 100., 3.5, 2.5, 1.5, 1.5, 0.5};
99 THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
100 NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
101 NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
102 NSparseHistTrk->SetBinEdges(0,ptBins);
103 NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
104 NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
105 NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
106 NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
107 NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
108 NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
109 NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
110 NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
111 NSparseHistTrk->GetAxis(5)->SetTitle("isph");
112 NSparseHistTrk->GetAxis(5)->SetName("isph");
113 NSparseHistTrk->GetAxis(6)->SetTitle("iswd");
114 NSparseHistTrk->GetAxis(6)->SetName("iswd");
115 NSparseHistTrk->GetAxis(7)->SetTitle("y");
116 NSparseHistTrk->GetAxis(7)->SetName("y");
117 fOutput->Add(NSparseHistTrk);
119 //dimensions of THnSparse for stack
120 const Int_t nvarst=5;
121 // pt cent IDgen isph y
122 Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2};
123 Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5};
124 Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5};
125 THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
126 NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
127 NSparseHistSt->SetBinEdges(0,ptBins);
128 NSparseHistSt->GetAxis(0)->SetName("pT_rec");
129 NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
130 NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
131 NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
132 NSparseHistSt->GetAxis(2)->SetName("ID_gen");
133 NSparseHistSt->GetAxis(3)->SetTitle("isph");
134 NSparseHistSt->GetAxis(3)->SetName("isph");
135 NSparseHistSt->GetAxis(4)->SetTitle("y");
136 NSparseHistSt->GetAxis(4)->SetName("y");
137 fOutput->Add(NSparseHistSt);
139 //dimensions of THnSparse for the normalization
140 const Int_t nvarev=3;
142 Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, fnNchBins};
143 Double_t xminHistRealEv[nvarev] = { 0., 0., 0.};
144 Double_t xmaxHistRealEv[nvarev] = { 100., 100., 2000.};
145 THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
146 NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
147 NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
148 NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
149 NSparseHistEv->GetAxis(1)->SetName("Q_vec");
150 NSparseHistEv->GetAxis(2)->SetTitle("N charged");
151 NSparseHistEv->GetAxis(2)->SetName("N_ch");
152 fOutput->Add(NSparseHistEv);
154 PostData(1, fOutput );
155 PostData(2, fEventCuts);
156 PostData(3, fTrackCuts);
157 PostData(4, fHelperPID);
160 //________________________________________________________________________
162 void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
164 //Printf("An event");
166 fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
168 AliWarning("ERROR: AliAODEvent not available \n");
172 if (strcmp(fAOD->ClassName(), "AliAODEvent"))
174 AliFatal("Not processing AODs");
177 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
180 if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
182 Double_t Qvec=0.;//in case of MC we save space in the memory
184 // if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
185 // else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
186 Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
189 Double_t Cent=fEventCuts->GetCent();
191 // First do MC to fill up the MC particle array
192 TClonesArray *arrayMC = 0;
195 arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
197 AliFatal("Error: MC particles branch not found!\n");
199 Int_t nMC = arrayMC->GetEntries();
200 for (Int_t iMC = 0; iMC < nMC; iMC++)
202 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
203 if(!partMC->Charge()) continue;//Skip neutrals
204 if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
206 //if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!
208 //pt cent IDgen isph y
210 varSt[0]=partMC->Pt();
212 varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
213 varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
214 varSt[4]=partMC->Y();
215 ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
217 //Printf("a particle");
222 //main loop on tracks
226 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
227 AliAODTrack* track = fAOD->GetTrack(iTracks);
228 if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
229 if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
230 if(!fFillOnlyEvents){
231 Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
232 Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
233 Int_t IDgen=kSpUndefined;//set if MC
238 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
240 AliError("Cannot get MC particle");
243 IDgen=fHelperPID->GetParticleSpecies(partMC);
244 isph=partMC->IsPhysicalPrimary();
245 iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions
248 //pt cent Q vec IDrec IDgen isph iswd y
250 varTrk[0]=track->Pt();
253 varTrk[3]=(Double_t)IDrec;
254 varTrk[4]=(Double_t)IDgen;
255 varTrk[5]=(Double_t)isph;
256 varTrk[6]=(Double_t)iswd;
258 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
260 //for nsigma PID fill double counting of ID
261 if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
263 HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
264 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
265 if(HasDC[ipart]==kTRUE){
266 varTrk[3]=(Double_t)ipart;
267 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
272 //fill all charged (3)
275 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
276 }//end if fFillOnlyEvents
281 } // end loop on tracks
287 ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
289 PostData(1, fOutput );
290 PostData(2, fEventCuts);
291 PostData(3, fTrackCuts);
292 PostData(4, fHelperPID);
295 //_________________________________________________________________
296 void AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)