]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraBoth.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding new options : asymetric y cut, selection of centrality estimator, plus pid...
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraBoth.cxx
CommitLineData
239a080a 1/**************************************************************************\r
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\r
7 * Permission to use, copy, modify and distribute this software and its *\r
8 * documentation strictly for non-commercial purposes is hereby granted *\r
9 * without fee, provided that the above copyright notice appears in all *\r
10 * copies and that both the copyright notice and this permission notice *\r
11 * appear in the supporting documentation. The authors make no claims *\r
12 * about the suitability of this software for any purpose. It is *\r
13 * provided "as is" without express or implied warranty. *\r
14 **************************************************************************/\r
15\r
16//-----------------------------------------------------------------\r
17// AliAnalysisTaskSpectraBoth class\r
18//-----------------------------------------------------------------\r
19\r
20#include "TChain.h"\r
21#include "TTree.h"\r
22#include "TLegend.h"\r
23#include "TH1F.h"\r
24#include "TH2F.h"\r
25#include "TCanvas.h"\r
26#include "AliAnalysisTask.h"\r
27#include "AliAnalysisManager.h"\r
28#include "AliAODTrack.h"\r
29#include "AliAODMCParticle.h"\r
30#include "AliVParticle.h"\r
31#include "AliAODEvent.h"\r
32#include "AliAODInputHandler.h"\r
33#include "AliAnalysisTaskSpectraBoth.h"\r
34#include "AliAnalysisTaskESDfilter.h"\r
35#include "AliAnalysisDataContainer.h"\r
36#include "AliSpectraBothHistoManager.h"\r
37#include "AliSpectraBothTrackCuts.h"\r
38#include "AliSpectraBothEventCuts.h"\r
39#include "AliCentrality.h"\r
40#include "TProof.h"\r
41#include "AliPID.h"\r
42#include "AliVEvent.h"\r
43#include "AliESDEvent.h"\r
44#include "AliPIDResponse.h"\r
45#include "AliStack.h"\r
46#include "AliSpectraBothPID.h"\r
47#include <TMCProcess.h>\r
48\r
49#include <iostream>\r
50\r
51\r
52\r
53\r
54using namespace AliSpectraNameSpaceBoth;\r
55using namespace std;\r
56\r
57ClassImp(AliAnalysisTaskSpectraBoth)\r
58\r
59//________________________________________________________________________\r
60AliAnalysisTaskSpectraBoth::AliAnalysisTaskSpectraBoth(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fPID(0), fIsMC(0), fNRebin(0),fUseMinSigma(0),fCuts(0),fdotheMCLoopAfterEventCuts(0)\r
61{\r
62 // Default constructor\r
63 \r
64 DefineInput(0, TChain::Class());\r
65 DefineOutput(1, AliSpectraBothHistoManager::Class());\r
66 DefineOutput(2, AliSpectraBothEventCuts::Class());\r
67 DefineOutput(3, AliSpectraBothTrackCuts::Class());\r
68 DefineOutput(4, AliSpectraBothPID::Class());\r
69 fNRebin=0;\r
70 \r
71}\r
72//________________________________________________________________________\r
73//________________________________________________________________________\r
74void AliAnalysisTaskSpectraBoth::UserCreateOutputObjects()\r
75{\r
76 // create output objects\r
77 fHistMan = new AliSpectraBothHistoManager("SpectraHistos",fNRebin);\r
78\r
79 if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
80 if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
81 if (!fPID) AliFatal("PID object should be set in the steering macro");\r
82 fTrackCuts->SetAliESDtrackCuts(fCuts);\r
83 PostData(1, fHistMan );\r
84 PostData(2, fEventCuts);\r
85 PostData(3, fTrackCuts);\r
86 PostData(4, fPID );\r
87\r
88}\r
89//________________________________________________________________________\r
90void AliAnalysisTaskSpectraBoth::UserExec(Option_t *)\r
91{\r
92 // main event loop\r
93 Int_t ifAODEvent=AliSpectraBothTrackCuts::kotherobject;\r
94 fAOD = dynamic_cast<AliVEvent*>(InputEvent());\r
95 \r
96 // AliESDEvent* esdevent=0x0;\r
97 // AliAODEvent* aodevent=0x0;\r
98 \r
99 TString nameoftrack(fAOD->ClassName()); \r
100 if(!nameoftrack.CompareTo("AliESDEvent"))\r
101 {\r
102 ifAODEvent=AliSpectraBothTrackCuts::kESDobject;\r
103 //esdevent=dynamic_cast<AliESDEvent*>(fAOD);\r
104 }\r
105 else if(!nameoftrack.CompareTo("AliAODEvent"))\r
106 {\r
107 ifAODEvent=AliSpectraBothTrackCuts::kAODobject;\r
108 //aodevent=dynamic_cast<AliAODEvent*>(fAOD);\r
109 }\r
110 else\r
111 AliFatal("Not processing AODs or ESDS") ;\r
112 if(fdotheMCLoopAfterEventCuts)\r
113 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
114 \r
115 \r
116 TClonesArray *arrayMC = 0;\r
117 Int_t npar=0;\r
118 AliStack* stack=0x0;\r
119 \r
120 if (fIsMC)\r
121 {\r
122 \r
123 if(ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
124 { \r
125 arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
126 if (!arrayMC) \r
127 {\r
128 AliFatal("Error: MC particles branch not found!\n");\r
129 }\r
130 Int_t nMC = arrayMC->GetEntries();\r
131 for (Int_t iMC = 0; iMC < nMC; iMC++)\r
132 {\r
133 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
134 if(!partMC->Charge()) continue;//Skip neutrals\r
135 //if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax()){//charged hadron are filled inside the eta acceptance\r
136 //Printf("%f %f-%f",partMC->Eta(),fTrackCuts->GetEtaMin(),fTrackCuts->GetEtaMax());\r
137 if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
138 fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
139 \r
140 //rapidity cut\r
b3ea73e1 141 if(partMC->Y() > fTrackCuts->GetYMax()|| partMC->Y() < fTrackCuts->GetYMin() ) \r
142 continue; \r
239a080a 143 if(partMC->IsPhysicalPrimary())\r
144 npar++; \r
145 // check for true PID + and fill P_t histos \r
146 Int_t charge = partMC->Charge() > 0 ? kChPos : kChNeg ;\r
147 Int_t id = fPID->GetParticleSpecie(partMC);\r
148 if(id != kSpUndefined) \r
149 {\r
150 fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
151 }\r
152 }\r
153 }\r
154 if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
155 { \r
156 AliMCEvent* mcEvent = (AliMCEvent*) MCEvent();\r
157 //Printf("MC particles: %d", mcEvent->GetNumberOfTracks()); \r
158 if (!mcEvent) \r
159 {\r
160 AliFatal("Error: MC particles branch not found!\n");\r
161 }\r
162 stack = mcEvent->Stack();\r
163 Int_t nMC = stack->GetNtrack();\r
164 for (Int_t iMC = 0; iMC < nMC; iMC++)\r
165 {\r
166 \r
167 TParticle *partMC = stack->Particle(iMC);\r
168 \r
169 if(!partMC) \r
170 continue; \r
171 \r
172 if(!partMC->GetPDG(0))\r
173 continue;\r
174 if(TMath::Abs(partMC->GetPDG(0)->Charge()/3.0)<0.01) \r
175 continue;//Skip neutrals\r
176 if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
177 fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
b3ea73e1 178 if(partMC->Y() > fTrackCuts->GetYMax() ||partMC->Y() < fTrackCuts->GetYMin() ) \r
239a080a 179 continue;\r
180 if(stack->IsPhysicalPrimary(iMC))\r
181 npar++; \r
182 // check for true PID + and fill P_t histos \r
183 Int_t charge = partMC->GetPDG(0)->Charge()/3.0 > 0 ? kChPos : kChNeg ;\r
184 Int_t id = fPID->GetParticleSpecie(partMC);\r
185 if(id != kSpUndefined) \r
186 {\r
187 fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
188 }\r
189 }\r
190 }\r
191 }\r
192 \r
193 if(!fdotheMCLoopAfterEventCuts)\r
194 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
195\r
196 //main loop on tracks\r
197 Int_t ntracks=0;\r
198 //cout<<fAOD->GetNumberOfTracks()<<endl;\r
199 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) \r
200 {\r
201 AliVTrack* track = dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));\r
202 AliAODTrack* aodtrack=0;\r
203 AliESDtrack* esdtrack=0;\r
204 Float_t dca=-999.;\r
205 if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
206 {\r
207 esdtrack=dynamic_cast<AliESDtrack*>(track);\r
737ad8ac 208 if(!esdtrack)\r
209 continue;\r
239a080a 210 Float_t dcaz=0.0;\r
211 esdtrack->GetImpactParameters(dca,dcaz);\r
212 }\r
213 else if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
214 {\r
737ad8ac 215 aodtrack=dynamic_cast<AliAODTrack*>(track);\r
216 if(!aodtrack)\r
217 continue; \r
239a080a 218 dca=aodtrack->DCA();\r
219 }\r
220 else\r
221 continue;\r
222 if (!fTrackCuts->IsSelected(track,kTRUE)) \r
223 continue;\r
224 ntracks++;\r
225 fPID->FillQAHistos(fHistMan, track, fTrackCuts);\r
226 \r
227 //calculate DCA for AOD track\r
228 if(dca==-999.)\r
229 {// track->DCA() does not work in old AOD production\r
230 Double_t d[2], covd[3];\r
231 AliVTrack* track_clone=(AliVTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
232 Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
233 delete track_clone;\r
234 if(!isDCA)\r
235 d[0]=-999.;\r
236 dca=d[0];\r
237 }\r
238 fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),dca); // PT histo\r
239 // get identity and charge\r
240 Bool_t rec[3]={false,false,false};\r
241 Int_t idRec = fPID->GetParticleSpecie(fHistMan,track, fTrackCuts,rec);\r
242 for(int irec=kSpPion;irec<kNSpecies;irec++)\r
243 {\r
244 \r
245 if(fUseMinSigma)\r
246 {\r
247 if(irec>kSpPion)\r
248 break;\r
249 }\r
250 else\r
251 { \r
252 if(!rec[irec]) \r
253 idRec = kSpUndefined;\r
254 else \r
255 idRec=irec;\r
256 } \r
257 \r
258 Int_t charge = track->Charge() > 0 ? kChPos : kChNeg;\r
259 \r
260 // Fill histograms, only if inside y and nsigma acceptance\r
261 if(idRec != kSpUndefined && fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec))\r
262 fHistMan->GetHistogram2D(kHistPtRecSigma,idRec,charge)->Fill(track->Pt(),dca);\r
263 //can't put a continue because we still have to fill allcharged primaries, done later\r
264 \r
265 /* MC Part */\r
266 if (arrayMC||stack) \r
267 {\r
268 Bool_t isPrimary = kFALSE;\r
269 Bool_t isSecondaryMaterial = kFALSE; \r
270 Bool_t isSecondaryWeak = kFALSE; \r
271 Int_t idGen =kSpUndefined;\r
272 Int_t pdgcode=0;\r
273 if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
274 {\r
275 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
276 if (!partMC) \r
277 { \r
278 AliError("Cannot get MC particle");\r
279 continue; \r
280 }\r
281 // Check if it is primary, secondary from material or secondary from weak decay\r
282 isPrimary = partMC->IsPhysicalPrimary();\r
283 isSecondaryWeak = partMC->IsSecondaryFromWeakDecay();\r
284 isSecondaryMaterial = partMC->IsSecondaryFromMaterial();\r
285 //cout<<"AOD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
286\r
287 if(!isPrimary&&!isSecondaryWeak&&!isSecondaryMaterial) \r
288 {\r
289 AliError("old tagging");\r
290 Int_t mfl=-999,codemoth=-999;\r
291 Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
292 if(indexMoth>=0)\r
293 {//is not fake\r
294 AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
295 codemoth = TMath::Abs(moth->GetPdgCode());\r
296 mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
297 }\r
298 //Int_t uniqueID = partMC->GetUniqueID();\r
299 //cout<<"uniqueID: "<<partMC->GetUniqueID()<<" "<<kPDecay<<endl;\r
300 //cout<<"status: "<<partMC->GetStatus()<<" "<<kPDecay<<endl;\r
301 // if(uniqueID == kPDecay)Printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");\r
302 if(mfl==3) \r
303 isSecondaryWeak = kTRUE; // add if(partMC->GetStatus() & kPDecay)? FIXME\r
304 else \r
305 isSecondaryMaterial = kTRUE;\r
306 }\r
307 //cout<<"AOD 2 tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
308\r
309\r
310 idGen = fPID->GetParticleSpecie(partMC);\r
311 pdgcode=partMC->GetPdgCode(); \r
312 }\r
313 else if (ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
314 {\r
315 TParticle *partMC =stack->Particle(TMath::Abs(track->GetLabel()));\r
316 if (!partMC) \r
317 { \r
318 AliError("Cannot get MC particle");\r
319 continue; \r
320 }\r
321 isPrimary = stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));\r
322 isSecondaryWeak = stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()));\r
323 isSecondaryMaterial = stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()));\r
324 //cout<<"ESD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<endl;\r
325 \r
326 idGen = fPID->GetParticleSpecie(partMC);\r
327 pdgcode=partMC->GetPdgCode(); \r
328 }\r
329 else\r
330 return;\r
331 \r
332 // cout<<isPrimary<<" "<<isSecondaryWeak<<" "<<isSecondaryMaterial<<endl;\r
333 // cout<<" functions "<<partMC->IsPhysicalPrimary()<<" "<<partMC->IsSecondaryFromWeakDecay()<<" "<<partMC->IsSecondaryFromMaterial()<<endl;\r
334 \r
335 if (isPrimary&&irec==kSpPion)\r
94a8ed71 336 fHistMan->GetPtHistogram(kHistPtRecPrimaryAll)->Fill(track->Pt(),dca); // PT histo of primaries\r
239a080a 337 \r
338 //nsigma cut (reconstructed nsigma)\r
339 if(idRec == kSpUndefined) \r
340 continue;\r
341 \r
342 // rapidity cut (reconstructed pt and identity)\r
343 if(!fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec)) continue;\r
344 \r
345 // Get true ID\r
346 \r
347 //if(TMath::Abs(partMC->Y()) > fTrackCuts->GetY() ) continue; // FIXME: do we need a rapidity cut on the generated?\r
348 // Fill histograms for primaries\r
349 \r
350 if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTrue, idGen, charge)->Fill(track->Pt(),dca); \r
351 \r
352 if (isPrimary) \r
353 {\r
354 fHistMan->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),dca); \r
355 if(idGen != kSpUndefined) \r
356 {\r
357 fHistMan->GetHistogram2D(kHistPtRecPrimary, idGen, charge)->Fill(track->Pt(),dca);\r
358 if (idRec == idGen) \r
359 fHistMan->GetHistogram2D(kHistPtRecTruePrimary, idGen, charge)->Fill(track->Pt(),dca); \r
360 }\r
361 }\r
362 //25th Apr - Muons are added to Pions -- FIXME\r
363 if ( pdgcode == 13 && idRec == kSpPion) \r
364 { \r
365 fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),dca); \r
366 if(isPrimary)\r
367 fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),dca); \r
368 }\r
369 if ( pdgcode == -13 && idRec == kSpPion) \r
370 { \r
371 fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),dca); \r
372 if (isPrimary) \r
373 {\r
374 fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),dca); \r
375 }\r
376 }\r
377 \r
378 ///..... END FIXME\r
379 \r
380 // Fill secondaries\r
381 if(isSecondaryWeak ) \r
382 fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),dca);\r
383 if(isSecondaryMaterial) \r
384 fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),dca);\r
385 \r
386 }//end if(arrayMC)\r
387 }\r
388 \r
389 \r
390 } // end loop on tracks\r
391 // cout<< ntracks<<endl;\r
392 fHistMan->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->Fill(npar,ntracks);\r
393 PostData(1, fHistMan );\r
394 PostData(2, fEventCuts);\r
395 PostData(3, fTrackCuts);\r
396 PostData(4, fPID );\r
397}\r
398\r
399//_________________________________________________________________\r
400void AliAnalysisTaskSpectraBoth::Terminate(Option_t *)\r
401{\r
402 // Terminate\r
403}\r