flag for double counting
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraAllChAOD.cxx
CommitLineData
c683985a 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
48using namespace AliHelperPIDNameSpace;
49using namespace std;
50
51ClassImp(AliAnalysisTaskSpectraAllChAOD)
52
53//________________________________________________________________________
54AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),
55 fAOD(0x0),
56 fTrackCuts(0x0),
57 fEventCuts(0x0),
58 fHelperPID(0x0),
59 fIsMC(0),
bc07eb4b 60 fDoDoubleCounting(0),
c683985a 61 fCharge(0),
62 fVZEROside(0),
63 fOutput(0x0),
64 fnCentBins(20),
65 fnQvecBins(40)
66{
67 // Default constructor
68 DefineInput(0, TChain::Class());
69 DefineOutput(1, TList::Class());
70 DefineOutput(2, AliSpectraAODEventCuts::Class());
71 DefineOutput(3, AliSpectraAODTrackCuts::Class());
72 DefineOutput(4, AliHelperPID::Class());
73}
74
75//________________________________________________________________________
76void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
77{
78 // create output objects
79 fOutput=new TList();
80 fOutput->SetOwner();
81 fOutput->SetName("fOutput");
82
83 if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
84 if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
85 if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");
86
87 // binning common to all the THn
88 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.};
89 const Int_t nptBins=26;
90
91 //dimensions of THnSparse for tracks
92 const Int_t nvartrk=8;
93 // pt cent Q vec IDrec IDgen isph iswd y
94 Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 2, 2, 2};
f99a88b2 95 Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, -0.5, -0.5, -0.5};
96 Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., 8., 3.5, 2.5, 1.5, 1.5, 0.5};
c683985a 97 THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
98 NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
99 NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
100 NSparseHistTrk->SetBinEdges(0,ptBins);
101 NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
102 NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
103 NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
104 NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
105 NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
106 NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
107 NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
108 NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
109 NSparseHistTrk->GetAxis(5)->SetTitle("isph");
110 NSparseHistTrk->GetAxis(5)->SetName("isph");
111 NSparseHistTrk->GetAxis(6)->SetTitle("iswd");
112 NSparseHistTrk->GetAxis(6)->SetName("iswd");
113 NSparseHistTrk->GetAxis(7)->SetTitle("y");
114 NSparseHistTrk->GetAxis(7)->SetName("y");
115 fOutput->Add(NSparseHistTrk);
116
117 //dimensions of THnSparse for stack
118 const Int_t nvarst=5;
119 // pt cent IDgen isph y
120 Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2};
121 Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5};
122 Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5};
123 THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
124 NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
125 NSparseHistSt->SetBinEdges(0,ptBins);
126 NSparseHistSt->GetAxis(0)->SetName("pT_rec");
127 NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
128 NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
129 NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
130 NSparseHistSt->GetAxis(2)->SetName("ID_gen");
131 NSparseHistSt->GetAxis(3)->SetTitle("isph");
132 NSparseHistSt->GetAxis(3)->SetName("isph");
133 NSparseHistSt->GetAxis(4)->SetTitle("y");
134 NSparseHistSt->GetAxis(4)->SetName("y");
135 fOutput->Add(NSparseHistSt);
136
137 //dimensions of THnSparse for the normalization
138 const Int_t nvarev=2;
139 // cent Q vec
140 Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins};
141 Double_t xminHistRealEv[nvarev] = { 0., 0.};
142 Double_t xmaxHistRealEv[nvarev] = { 100., 10.};
143 THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
144 NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
145 NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
146 NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
147 NSparseHistEv->GetAxis(1)->SetName("Q_vec");
148 fOutput->Add(NSparseHistEv);
149
150 PostData(1, fOutput );
151 PostData(2, fEventCuts);
152 PostData(3, fTrackCuts);
153 PostData(4, fHelperPID);
154}
155
156//________________________________________________________________________
157
158void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
159{
160 //Printf("An event");
161 // main event loop
162 fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
163 if (!fAOD) {
164 AliWarning("ERROR: AliAODEvent not available \n");
165 return;
166 }
167
168 if (strcmp(fAOD->ClassName(), "AliAODEvent"))
169 {
170 AliFatal("Not processing AODs");
171 }
172
173 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
174
175 //Default TPC priors
176 if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
177
178 Double_t Qvec=0.;//in case of MC we save space in the memory
179 if(!fIsMC){
180 if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
181 else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
182 }
183
184 Double_t Cent=fEventCuts->GetCent();
185
186 // First do MC to fill up the MC particle array
187 TClonesArray *arrayMC = 0;
188 if (fIsMC)
189 {
190 arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
191 if (!arrayMC) {
192 AliFatal("Error: MC particles branch not found!\n");
193 }
194 Int_t nMC = arrayMC->GetEntries();
195 for (Int_t iMC = 0; iMC < nMC; iMC++)
196 {
197 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
198 if(!partMC->Charge()) continue;//Skip neutrals
660a66d4 199 if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
c683985a 200
201 if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!
202
203 //pt cent IDgen isph y
204 Double_t varSt[5];
205 varSt[0]=partMC->Pt();
206 varSt[1]=Cent;
207 varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
208 varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
209 varSt[4]=partMC->Y();
210 ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
211
212 //Printf("a particle");
213
214 }
215 }
216
217 //main loop on tracks
218 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
219 AliAODTrack* track = fAOD->GetTrack(iTracks);
220 if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
221 if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
222 Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
223 Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
224 Int_t IDgen=kSpUndefined;//set if MC
225 Int_t isph=-999;
226 Int_t iswd=-999;
227
228 if (arrayMC) {
229 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
230 if (!partMC) {
231 AliError("Cannot get MC particle");
232 continue;
233 }
234 IDgen=fHelperPID->GetParticleSpecies(partMC);
235 isph=partMC->IsPhysicalPrimary();
236 iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions
237 }
238
239 //pt cent Q vec IDrec IDgen isph iswd y
240 Double_t varTrk[8];
241 varTrk[0]=track->Pt();
242 varTrk[1]=Cent;
243 varTrk[2]=Qvec;
244 varTrk[3]=(Double_t)IDrec;
245 varTrk[4]=(Double_t)IDgen;
246 varTrk[5]=(Double_t)isph;
247 varTrk[6]=(Double_t)iswd;
248 varTrk[7]=y;
249 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
250
251 //for nsigma PID fill double counting of ID
bc07eb4b 252 if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
c683985a 253 Bool_t *HasDC;
254 HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
255 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
256 if(HasDC[ipart]==kTRUE){
257 varTrk[3]=(Double_t)ipart;
258 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
259 }
260 }
261 }
262
f99a88b2 263 //fill all charged (4)
264 varTrk[3]=4.;
c683985a 265 ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
266
267 //Printf("a track");
268
269 } // end loop on tracks
270
271 Double_t varEv[2];
272 varEv[0]=Cent;
273 varEv[1]=Qvec;
274 ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
275
276 PostData(1, fOutput );
277 PostData(2, fEventCuts);
278 PostData(3, fTrackCuts);
279 PostData(4, fHelperPID);
280}
281
282//_________________________________________________________________
283void AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
284{
285 // Terminate
286}