]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2AOD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / AliAnalysisTaskB2AOD.cxx
CommitLineData
9e44412b 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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// Analysis task for B2 (AOD)
17// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
18
19#include <AliAnalysisTaskSE.h>
20#include <AliAnalysisManager.h>
21#include <AliInputEventHandler.h>
22#include <AliExternalTrackParam.h>
23#include <AliAODEvent.h>
24#include <AliAODVertex.h>
25#include <AliAODTrack.h>
26#include <AliAODMCParticle.h>
27#include <TH1D.h>
28#include <TH2D.h>
29#include <TList.h>
30#include <TString.h>
31#include <TVector3.h>
32#include <TProfile.h>
33#include <AliLog.h>
34
35#include "AliLnID.h"
36#include "AliLnHistoMap.h"
37#include "AliLnAODtrackCuts.h"
38#include "AliAnalysisTaskB2AOD.h"
39
40ClassImp(AliAnalysisTaskB2AOD)
41
42AliAnalysisTaskB2AOD::AliAnalysisTaskB2AOD()
43: AliAnalysisTaskSE()
44, fSpecies("Proton")
45, fPartCode(AliPID::kProton)
46, fHeavyIons(0)
47, fSimulation(0)
48, fMultTriggerFired(0)
49, fCentTriggerFired(0)
50, fTriggerFired(0)
51, fGoodVertex(0)
52, fPileUpEvent(0)
53, fV0AND(0)
54, fNoFastOnly(0)
55, fNtrkMultTrigger(0)
56, fMinKNOmult(-10)
57, fMaxKNOmult(100000)
58, fMinCentrality(0)
59, fMaxCentrality(100)
60, fNch(0)
61, fNtrk(0)
62, fMeanNtrk(1)
63, fKNOmult(1)
64, fMinVx(-1)
65, fMaxVx(1)
66, fMinVy(-1)
67, fMaxVy(1)
68, fMinVz(-10)
69, fMaxVz(10)
70, fMinDCAxy(-1)
71, fMaxDCAxy(1)
72, fMinDCAz(-2)
73, fMaxDCAz(2)
74, fMaxNSigma(3)
75, fMinEta(-0.8)
76, fMaxEta(0.8)
77, fMinY(-0.5)
78, fMaxY(0.5)
79, fAODevent(0)
80, fOutputContainer(0)
81, fHistoMap(0)
82, fTrackCuts(0)
83, fLnID(0)
84, fMaxNSigmaITS(3)
85, fMaxNSigmaTPC(3)
86, fMaxNSigmaTOF(3)
87, fMinM2(2.)
88, fMaxM2(6.)
89, fMomentumCorrection(0)
90, fMoCpfx(0)
91
92{
93//
94// Default constructor
95//
96 AliLog::SetGlobalLogLevel(AliLog::kFatal);
97}
98
99AliAnalysisTaskB2AOD::AliAnalysisTaskB2AOD(const char* name)
100: AliAnalysisTaskSE(name)
101, fSpecies("Proton")
102, fPartCode(AliPID::kProton)
103, fHeavyIons(0)
104, fSimulation(0)
105, fMultTriggerFired(0)
106, fCentTriggerFired(0)
107, fTriggerFired(0)
108, fGoodVertex(0)
109, fPileUpEvent(0)
110, fV0AND(0)
111, fNoFastOnly(0)
112, fNtrkMultTrigger(0)
113, fMinKNOmult(-10)
114, fMaxKNOmult(100000)
115, fMinCentrality(-1)
116, fMaxCentrality(100)
117, fNch(0)
118, fNtrk(0)
119, fMeanNtrk(1)
120, fKNOmult(1)
121, fMinVx(-1)
122, fMaxVx(1)
123, fMinVy(-1)
124, fMaxVy(1)
125, fMinVz(-10)
126, fMaxVz(10)
127, fMinDCAxy(-1)
128, fMaxDCAxy(1)
129, fMinDCAz(-2)
130, fMaxDCAz(2)
131, fMaxNSigma(3)
132, fMinEta(-0.8)
133, fMaxEta(0.8)
134, fMinY(-0.5)
135, fMaxY(0.5)
136, fAODevent(0)
137, fOutputContainer(0)
138, fHistoMap(0)
139, fTrackCuts(0)
140, fLnID(0)
141, fMaxNSigmaITS(3)
142, fMaxNSigmaTPC(3)
143, fMaxNSigmaTOF(3)
144, fMinM2(2.)
145, fMaxM2(6.)
146, fMomentumCorrection(0)
147, fMoCpfx(0)
148
149{
150//
151// Constructor
152//
153 DefineOutput(1, TList::Class());
154
155 AliLog::SetGlobalLogLevel(AliLog::kFatal);
156}
157
158void AliAnalysisTaskB2AOD::UserCreateOutputObjects()
159{
160//
161// Create output objects
162//
163 if(fHistoMap == 0) AliFatal("no histogram map"); // should be created somewhere else
164
165 fOutputContainer = new TList();
166 fOutputContainer->SetOwner(kTRUE);
167
168 TObjString* key;
169 TIter iter(fHistoMap->GetMap());
170 while((key = dynamic_cast<TObjString*>(iter.Next()))) fOutputContainer->Add((TH1*)fHistoMap->Get(key));
171
172 PostData(1, fOutputContainer);
173}
174
175AliAnalysisTaskB2AOD::~AliAnalysisTaskB2AOD()
176{
177//
178// Default destructor
179//
180 delete fHistoMap;
181 delete fOutputContainer;
182 delete fTrackCuts;
183 delete fLnID;
184
185 delete fMoCpfx;
186}
187
188void AliAnalysisTaskB2AOD::SetParticleSpecies(const TString& species)
189{
190//
191// set the particle species and the AliPID code
192//
193 fSpecies = species;
194 fPartCode = this->GetPidCode(species);
195}
196
197void AliAnalysisTaskB2AOD::UserExec(Option_t* )
198{
199//
200// Execute analysis for the current event
201//
202 fAODevent = dynamic_cast<AliAODEvent*>(InputEvent());
203
204 if (!fAODevent)
205 {
206 this->Error("Exec", "AOD event not available");
207 return;
208 }
209
210 if(fTrackCuts == 0) AliFatal("track cuts not set");
211
212 if(fLnID == 0) AliFatal("PID not set");
213
214 // multiplicity and centrality
0a918d8d 215 AliAODHeader * header = dynamic_cast<AliAODHeader*>(fAODevent->GetHeader());
216 if(!header) AliFatal("Not a standard AOD");
217
218 fNtrk = (fMaxEta > 0.5) ? header->GetRefMultiplicityComb08() : header->GetRefMultiplicityComb05();
9e44412b 219
220 if(fSimulation) fNch = this->GetChargedMultiplicity(fMaxEta);
221
222 fKNOmult = fNtrk/fMeanNtrk;
223
224 ((TH1D*)fHistoMap->Get(fSpecies + "_Event_Ntrk"))->Fill(fNtrk);
225 ((TH1D*)fHistoMap->Get(fSpecies + "_Event_Zmult"))->Fill(fKNOmult);
226
227 fMultTriggerFired = (fNtrk > 0) && (fKNOmult >= fMinKNOmult) && (fKNOmult < fMaxKNOmult);
228
229 if(fHeavyIons)
230 {
0a918d8d 231 Double_t centrality = ((AliVAODHeader*)fAODevent->GetHeader())->GetCentrality();
9e44412b 232
233 fCentTriggerFired = (centrality >= fMinCentrality) && (centrality < fMaxCentrality);
234 }
235
236 // trigger
237
238 AliInputEventHandler* eventH = dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
239 if(eventH == 0)
240 {
241 this->Error("Exec", "could not get AliInputEventHandler");
242 return;
243 }
244
245 UInt_t triggerBits = eventH->IsEventSelected();
246
247 if(fHeavyIons)
248 {
249 fTriggerFired = ( this->IsMB(triggerBits) && fCentTriggerFired );
250 }
251 else
252 {
253 fTriggerFired = this->IsMB(triggerBits);
254 if(fNoFastOnly) fTriggerFired = ( fTriggerFired && !this->IsFastOnly(triggerBits) );
255 if(fV0AND) fTriggerFired = ( fTriggerFired && this->IsV0AND(triggerBits) );
256 if(fNtrkMultTrigger) fTriggerFired = ( fTriggerFired && fMultTriggerFired );
257 }
258
259 // vertex
260
261 fGoodVertex = kFALSE;
262
263 const AliAODVertex* vtx = fAODevent->GetPrimaryVertex(); // best primary vertex
264
265 if(vtx->GetNContributors()>0)
266 {
267 fGoodVertex=kTRUE;
268
269 ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_YX"))->Fill(vtx->GetX(), vtx->GetY());
270 ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_YZ"))->Fill(vtx->GetZ(), vtx->GetY());
271 ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_XZ"))->Fill(vtx->GetZ(), vtx->GetX());
272 }
273
274 if( (vtx->GetX() <= fMinVx) || (vtx->GetX() >= fMaxVx) ) fGoodVertex=kFALSE;
275 if( (vtx->GetY() <= fMinVy) || (vtx->GetY() >= fMaxVy) ) fGoodVertex=kFALSE;
276 if( (vtx->GetZ() <= fMinVz) || (vtx->GetZ() >= fMaxVz) ) fGoodVertex=kFALSE;
277
278 // pile up
279
280 fPileUpEvent = fAODevent->IsPileupFromSPDInMultBins();
281
282 // event stats
283
284 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(0); // number of events
285
286 if(fTriggerFired)
287 {
288 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(1); // triggering events
289
290 if(vtx->GetNContributors()>0)
291 {
292 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(3); // valid vertex within a triggering event
293
294 if((vtx->GetZ() > fMinVz) && (vtx->GetZ() < fMaxVz))
295 {
296 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(4); // Vz
297
298 if( (vtx->GetX() > fMinVx) && (vtx->GetX() < fMaxVx)
299 && (vtx->GetY() > fMinVy) && (vtx->GetY() < fMaxVy))
300 {
301 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(5); // Vx, Vy
302
303 if(fGoodVertex)
304 {
305 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(6); // VertexerZ
306
307 if(!fPileUpEvent)
308 {
309 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(7); // no pile-up
310 }
311 }
312 }
313 }
314 }
315 }
316
317 // Particles and tracks for this event
318
319 if(fSimulation) this->GetParticles();
320
321 this->GetTracks();
322
323 // Post the data (input/output slots #0 already used by the base class)
324 PostData(1, fOutputContainer);
325}
326
327Int_t AliAnalysisTaskB2AOD::GetParticles()
328{
329//
330// Get particles from current event
331//
332 Int_t nParticles = 0;
333
334 TList* lst = fAODevent->GetList();
335 TClonesArray* lop = dynamic_cast<TClonesArray*>(lst->FindObject(AliAODMCParticle::StdBranchName()));
336
337 if (!lop)
338 {
339 AliDebug(AliLog::kWarning, "list of particles not available");
340 return 0;
341 }
342
343 for (Int_t i = 0; i < lop->GetEntries(); ++i)
344 {
345 AliAODMCParticle* iParticle = dynamic_cast<AliAODMCParticle*>(lop->At(i));
346
347 if(!iParticle) continue;
348
349 Int_t pid = fLnID->GetPID(iParticle);
350
351 if(pid != fPartCode) continue;
352
353 // physical primaries
354
355 if(!iParticle->IsPhysicalPrimary()) continue;
356
357 TString particle = fSpecies;
358 if(iParticle->Charge() < 0) particle.Prepend("Anti");
359
360 Double_t genP = iParticle->P();
361 Double_t genPt = iParticle->Pt();
362 Double_t genY = iParticle->Y();
363 Double_t genPhi = iParticle->Phi();
364 Double_t genEta = iParticle->Eta();
365
366 // multiplicity and centrality
367
368 if( fHeavyIons && !fCentTriggerFired) continue;
369 if( fNtrkMultTrigger && !fMultTriggerFired) continue;
370
371 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_P"))->Fill(genP);
372 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Pt"))->Fill(genPt);
373 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Y"))->Fill(genY);
374 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Phi"))->Fill(genPhi);
375 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Eta"))->Fill(genEta);
376 ((TH2D*)fHistoMap->Get(particle + "_Gen_Prim_PtY"))->Fill(genY, genPt);
377 ((TH2D*)fHistoMap->Get(particle + "_Gen_Prim_EtaY"))->Fill(genY, genEta);
378
379 // is within phase space?
380
381 if(TMath::Abs(genY) >= fMaxY) continue;
382
383 ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_P"))->Fill(genP);
384 ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_Pt"))->Fill(genPt);
385
386 // is from a triggering event? (same as rec.)
387
388 if(!fTriggerFired)
389 {
390 ((TH1D*)fHistoMap->Get(particle + "_Gen_NoTrig_Prim_Pt"))->Fill(genPt);
391 continue;
392 }
393
394 ((TH1D*)fHistoMap->Get(particle + "_Gen_Trig_Prim_Pt"))->Fill(genPt);
395
396 // is from a triggering event with good vertex?
397
398 if(fAODevent->GetPrimaryVertex()->GetNContributors() <= 0)
399 {
400 ((TH1D*)fHistoMap->Get(particle + "_Gen_NoVtx_Prim_Pt"))->Fill(genPt);
401 }
402
403 if(!fGoodVertex) continue;
404 if(fPileUpEvent) continue;
405
406 ((TH1D*)fHistoMap->Get(particle + "_Gen_Nch"))->Fill(fNch);
407 ((TH1D*)fHistoMap->Get(particle + "_Gen_Vtx_Prim_Pt"))->Fill(genPt);
408
409 // is within the geometrical acceptance?
410
411 if(!fTrackCuts->IsWithinGeoAcceptance(iParticle)) continue;
412
413 ((TH2D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_PtY"))->Fill(genY,genPt);
414 ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_P"))->Fill(genP);
415 ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Pt"))->Fill(genPt);
416 ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Phi"))->Fill(genPhi);
417 ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Y"))->Fill(genY);
418 }
419
420 return nParticles;
421}
422
423Int_t AliAnalysisTaskB2AOD::GetTracks()
424{
425//
426// Get tracks from current event
427//
428 using namespace std;
429
430 Int_t nTracks = 0;
431
432 // trigger, vertex and pile-up
433
434 if(!fTriggerFired) return 0;
435 if(!fGoodVertex) return 0;
436 if(fPileUpEvent) return 0;
437
438 // this is a 'good' event
439
440 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(2); // analyzed events
441
442 ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_Event_Ntrk"))->Fill(fNtrk);
443 ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_Event_Zmult"))->Fill(fKNOmult);
444
445 if(fSimulation)
446 {
447 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Event_Nch_Ntrk"))->Fill(fNtrk, fNch);
448 }
449
450 const AliAODVertex* vtx = fAODevent->GetPrimaryVertex();
451
452 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_YX"))->Fill(vtx->GetX(), vtx->GetY());
453 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_YZ"))->Fill(vtx->GetZ(), vtx->GetY());
454 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_XZ"))->Fill(vtx->GetZ(), vtx->GetX());
455
456 // track loop
457 for(Int_t i = 0; i < fAODevent->GetNumberOfTracks(); ++i)
458 {
459 AliAODTrack* iTrack = fAODevent->GetTrack(i);
460
461 if(!iTrack) continue;
462
463 TString particle = fSpecies;
464 if(iTrack->Charge() < 0) particle.Prepend("Anti");
465
466 // track parameters at the DCA to the primary vertex
467
468 AliExternalTrackParam trkparam;
469 trkparam.CopyFromVTrack(iTrack);
470
471 if(trkparam.GetX() > 3.) continue; // only valid for propagation inside the beam pipe
472
473 // impact parameters
474 Double_t b[2];
475 Double_t bCov[3];
476 if(!trkparam.PropagateToDCA(fAODevent->GetPrimaryVertex(), fAODevent->GetMagneticField(), 10000., b, bCov)) continue;
477
478 Double_t dcaxy = b[0];
479 Double_t dcaz = b[1];
480
481 if(iTrack->Charge() < 0) // in case of asymmetry
482 {
483 dcaxy = -dcaxy;
484 dcaz = -dcaz;
485 }
486
487 Double_t nSigmaVtx = fTrackCuts->GetNSigmaToVertex(b, bCov);
488
489 // momentum
368c671a 490 Double_t pDCA[3] = { iTrack->Px(), iTrack->Py(), iTrack->Pz() };
491 //trkparam.GetPxPyPz(pDCA);
9e44412b 492
493 Double_t p = TMath::Sqrt(pDCA[0]*pDCA[0] + pDCA[1]*pDCA[1] + pDCA[2]*pDCA[2]);
494 Double_t pt = TMath::Sqrt(pDCA[0]*pDCA[0] + pDCA[1]*pDCA[1]);
495 Double_t pz = pDCA[2];
496
497 // track cuts
498
499 TVector3 pVect(pDCA[0],pDCA[1],pDCA[2]);
500 Double_t eta = pVect.Eta();
501 Double_t phi = this->GetPhi(pDCA);
502
503 ((TH2D*)fHistoMap->Get(particle + "_Before_Phi_Eta"))->Fill(eta,phi);
504
505 if(!fTrackCuts->AcceptTrack(iTrack,b, bCov)) continue; // with next track
506 if(!fTrackCuts->IsWithinGeoAcceptance(pDCA)) continue; // with next track
507
508 // end track cuts
509
510 ++nTracks;
511
512 ((TH2D*)fHistoMap->Get(particle + "_After_Phi_Eta"))->Fill(eta, phi);
513
514 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_DCAxy"))->Fill(dcaxy);
515 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_DCAz"))->Fill(dcaz);
516 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_NSigma"))->Fill(nSigmaVtx);
517 //((TH1D*)fHistoMap->Get(particle + "_TrackCuts_ITSchi2PerCls"))->Fill(this->GetITSchi2PerCluster(iTrack));
518 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCncls"))->Fill(iTrack->GetTPCNcls());
519 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCxRowsOverF"))->Fill(fTrackCuts->GetNTPCXRowsOverFindable(iTrack));
520 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCxRows"))->Fill(iTrack->GetTPCNCrossedRows());
521 //((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCchi2PerCls"))->Fill(iTrack->GetTPCchi2()/iTrack->GetTPCNcls());
522 //((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCchi2Global"))->Fill(iTrack->GetChi2TPCConstrainedVsGlobal(fAODevent->GetPrimaryVertex()));
523
524 // detector signals
525
526 Double_t pITS = p;
527 Double_t pTPC = iTrack->GetTPCmomentum();
528 Double_t pTOF = p;
529 Double_t dEdxITS = iTrack->GetITSsignal();
530 Double_t dEdxTPC = iTrack->GetTPCsignal();
531 Int_t nPointsITS = this->GetITSnPointsPID(iTrack);
532 Int_t nPointsTPC = iTrack->GetTPCsignalN();
533 Double_t beta = 0;
534 Double_t mass = 0;
535 Double_t m2 = 0;
536
537 ((TH2D*)fHistoMap->Get(particle + "_ITS_dEdx_P"))->Fill(pITS, dEdxITS);
538 ((TH2D*)fHistoMap->Get(particle + "_TPC_dEdx_P"))->Fill(pTPC, dEdxTPC);
539
540 if(fTrackCuts->TOFmatch())
541 {
542 beta = this->GetBeta(iTrack);
543 m2 = this->GetMassSquared(p, beta);
544 mass = TMath::Sqrt(TMath::Abs(m2));
545
546 ((TH2D*)fHistoMap->Get(particle + "_TOF_Beta_P"))->Fill(pTOF, beta);
547 ((TH2D*)fHistoMap->Get(particle + "_TOF_Mass_P"))->Fill(pTOF, mass);
548 }
549
368c671a 550 // get the pid
551 Int_t pid = fLnID->GetPID( fPartCode, pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta, fMaxNSigmaITS, fMaxNSigmaTPC, fMaxNSigmaTOF);
552
553 Int_t offset = AliPID::kDeuteron - 5;
554
555 // for bayes iteration
556 if(pid != -1)
557 {
558 Int_t iBin = (pid > AliPID::kProton) ? pid - offset : pid;
559 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats_PID"))->Fill(iBin);
560 }
561
562 // fix momentum
563
564 if(fPartCode > AliPID::kTriton)
565 {
566 p *= 2.;
567 pt *= 2.;
568 pz *= 2.;
569 }
570
571 if(fMomentumCorrection)
572 {
573 pt += this->GetMomentumCorrection(pt);
574 p = TMath::Sqrt(pt*pt + pz*pz);
575
576 if(fTrackCuts->TOFmatch())
577 {
578 m2 = this->GetMassSquared(p, beta);
579 mass = TMath::Sqrt(TMath::Abs(m2));
580 }
581 }
582
9e44412b 583 // for pid and efficiency
584
585 Double_t simPt = 0;
586 Double_t simPhi = 0;
587 Double_t simY = 0;
588
589 Int_t simpid = -1;
590 AliAODMCParticle* iParticle = 0;
591
592 if(fSimulation)
593 {
594 iParticle = this->GetParticle(iTrack);
595
596 if(iParticle == 0) continue;
597
598 simPt = iParticle->Pt();
599 simPhi = iParticle->Phi();
600 simY = iParticle->Y();
601
602 simpid = fLnID->GetPID(iParticle);
603
604 if(simpid == fPartCode)
605 {
606 TString simparticle = fSpecies;
607 if(iParticle->Charge()<0) simparticle.Prepend("Anti");
608
609 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_PtY"))->Fill(simY, simPt);
610
611 if(TMath::Abs(simY) < fMaxY)
612 {
613 ((TH2D*)fHistoMap->Get(simparticle + "_Response_Matrix"))->Fill(pt, simPt);
614
615 if(iParticle->IsPhysicalPrimary())
616 {
617 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Ntrk"))->Fill(fNtrk);
618 }
619
620 // for pid
621 if(!this->IsFakeTrack(iTrack))
622 {
623 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Pt"))->Fill(simPt);
624
625 if(iParticle->IsPhysicalPrimary()) // the efficiency is calculated on the primaries
626 {
627 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Pt"))->Fill(simPt);
628 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Y"))->Fill(simY);
629 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Phi"))->Fill(simPhi);
630 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Rec_Pt"))->Fill(pt);
631
632 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_DCAxy_Pt"))->Fill(simPt,dcaxy);
633
634 ((TH2D*)fHistoMap->Get(simparticle + "_Prim_Response_Matrix"))->Fill(pt, simPt);
635
636 ((TH2D*)fHistoMap->Get(simparticle + "_Prim_DiffPt_RecPt"))->Fill(pt,simPt-pt);
637
638 if( fTrackCuts->TOFmatch() )
639 {
640 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_P"))->Fill(pTOF, m2);
641 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_Pt"))->Fill(pt, m2);
642 }
643 }
644 else if(iParticle->IsSecondaryFromWeakDecay())
645 {
646 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_Pt"))->Fill(simPt);
647
648 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_DCAxy_Pt"))->Fill(simPt,dcaxy);
649 }
650 else
651 {
652 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Mat_Pt"))->Fill(simPt);
653
654 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Mat_DCAxy_Pt"))->Fill(simPt,dcaxy);
655 }
656 }
657 else // fake tracks
658 {
659 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Pt"))->Fill(simPt);
660 if(iParticle->IsPhysicalPrimary())
661 {
662 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Prim_Pt"))->Fill(simPt);
663 }
664 else if(iParticle->IsSecondaryFromWeakDecay())
665 {
666 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Fdwn_Pt"))->Fill(simPt);
667 }
668 else
669 {
670 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Mat_Pt"))->Fill(simPt);
671 }
672 }
673 }
674 }
368c671a 675
676 // pid table for prior probabilities (only Bayes)
677
9e44412b 678 Int_t sim = (simpid > AliPID::kProton) ? simpid - offset : simpid;
679 Int_t rec = (pid > AliPID::kProton) ? pid - offset : pid;
680
681 if((sim > -1) && (rec > -1)) // pid performance
682 {
683 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(sim, rec);
684 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(9, rec);
685 }
686 if(sim > -1)
687 {
688 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(sim, 9);
689 }
690 }
691
368c671a 692 // candidate tracks
9e44412b 693
694 if(pid != fPartCode) continue;
695
9e44412b 696 Bool_t goodPid = 0;
697
698 if(fSimulation)
699 {
700 goodPid = ( simpid == pid );
701 }
702
703 ((TH2D*)fHistoMap->Get(particle + "_PID_Ntrk_pTPC"))->Fill(pTPC,fNtrk);
704 ((TH2D*)fHistoMap->Get(particle + "_PID_Zmult_pTPC"))->Fill(pTPC,fKNOmult);
705
706 // pid performance
707 ((TH2D*)fHistoMap->Get(particle + "_PID_ITSdEdx_P"))->Fill(pITS, dEdxITS);
708 ((TH2D*)fHistoMap->Get(particle + "_PID_TPCdEdx_P"))->Fill(pTPC, dEdxTPC);
709
710 if(fTrackCuts->TOFmatch())
711 {
712 ((TH2D*)fHistoMap->Get(particle + "_PID_Beta_P"))->Fill(pTOF, beta);
713 ((TH2D*)fHistoMap->Get(particle + "_PID_Mass_P"))->Fill(pTOF, mass);
714 if(fSimulation && goodPid)
715 {
716 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Mass"))->Fill(mass);
717 }
718 }
719
9e44412b 720 Double_t y = this->GetRapidity(p, pz, AliPID::ParticleMass(fPartCode));
721
722 ((TH1D*)fHistoMap->Get(particle + "_PID_Y"))->Fill(y);
723 ((TH2D*)fHistoMap->Get(particle + "_PID_Pt_Y"))->Fill(y, pt);
724
368c671a 725 // results in |y| < fMaxY
726
9e44412b 727 if(TMath::Abs(y) >= fMaxY) continue;
728
729 ((TH1D*)fHistoMap->Get(particle + "_PID_Pt"))->Fill(pt);
730 ((TH1D*)fHistoMap->Get(particle + "_PID_Phi"))->Fill(phi);
731
732 if(iTrack->IsOn(AliAODTrack::kTRDin))
733 {
734 ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_Pt"))->Fill(pt);
735
736 if(iTrack->IsOn(AliAODTrack::kTRDout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_TRDout_Pt"))->Fill(pt);
737 if(iTrack->IsOn(AliAODTrack::kTOFout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_TOFout_Pt"))->Fill(pt);
738 }
739
740 if(iTrack->IsOn(AliAODTrack::kTOFin))
741 {
742 ((TH1D*)fHistoMap->Get(particle + "_PID_TOFin_Pt"))->Fill(pt);
743
744 if(iTrack->IsOn(AliAODTrack::kTOFout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TOFin_TOFout_Pt"))->Fill(pt);
745 }
746
747 if(fTrackCuts->TOFmatch())
748 {
749 Double_t dm2 = this->GetM2Difference(beta, p, AliPID::ParticleMass(fPartCode));
750 Double_t t = iTrack->GetTOFsignal()*1.e-3; // ns
751 Double_t dt = t - this->GetExpectedTime(iTrack, AliPID::ParticleMass(fPartCode))*1.e-3;
752
753 ((TH2D*)fHistoMap->Get(particle + "_PID_M2_Pt"))->Fill(pt, m2);
754 ((TH2D*)fHistoMap->Get(particle + "_PID_DM2_Pt"))->Fill(pt, dm2);
755 ((TH2D*)fHistoMap->Get(particle + "_PID_Time_Pt"))->Fill(pt, t);
756 ((TH2D*)fHistoMap->Get(particle + "_PID_DTime_Pt"))->Fill(pt, dt);
757 ((TH1D*)fHistoMap->Get(particle + "_PID_TOFmatch_Pt"))->Fill(pt);
758 }
759
760 // secondaries
761
762 Bool_t m2match = kTRUE;
763
764 if( fTrackCuts->TOFmatch() && (fLnID->GetPidProcedure() > AliLnID::kMaxLikelihood))
765 {
766 if ((m2 < fMinM2) || (m2 >= fMaxM2)) m2match = kFALSE;
767 }
768
769 if(m2match)
770 {
771 ((TH2D*)fHistoMap->Get(particle + "_PID_DCAxy_Pt"))->Fill(pt, dcaxy);
772 ((TH2D*)fHistoMap->Get(particle + "_PID_DCAz_Pt"))->Fill(pt, dcaz);
773 ((TH2D*)fHistoMap->Get(particle + "_PID_NSigma_Pt"))->Fill(pt, nSigmaVtx);
774 }
775
776 if(fSimulation && goodPid)
777 {
778 // for unfolding and pid contamination
779 if(!this->IsFakeTrack(iTrack))
780 {
781 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Pt"))->Fill(simPt);
782
783 if( fTrackCuts->TOFmatch() )
784 {
785 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_M2_Pt"))->Fill(pt, m2);
786 }
787
788 if(iParticle->IsPhysicalPrimary())
789 {
790 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Prim_Pt"))->Fill(simPt);
791 }
792
793 if(m2match)
794 {
795 if(iParticle->IsPhysicalPrimary())
796 {
797 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_DCAxy_Pt"))->Fill(pt, dcaxy);
798 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_DCAz_Pt"))->Fill(pt, dcaz);
799 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_NSigma_Pt"))->Fill(pt, nSigmaVtx);
800 }
801 else if(iParticle->IsSecondaryFromWeakDecay())
802 {
803 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_DCAxy_Pt"))->Fill(pt, dcaxy);
804 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_DCAz_Pt"))->Fill(pt, dcaz);
805 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_NSigma_Pt"))->Fill(pt, nSigmaVtx);
806 }
807 else // from materials
808 {
809 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_DCAxy_Pt"))->Fill(pt, dcaxy);
810 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_DCAz_Pt"))->Fill(pt, dcaz);
811 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_NSigma_Pt"))->Fill(pt, nSigmaVtx);
812 }
813 }
814 }
815 else // fake tracks
816 {
817 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Pt"))->Fill(simPt);
818
819 if(m2match)
820 {
821 if(iParticle->IsPhysicalPrimary())
822 {
823 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Prim_DCAxy_Pt"))->Fill(pt, dcaxy);
824 }
825 else if(iParticle->IsSecondaryFromWeakDecay())
826 {
827 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt"))->Fill(pt, dcaxy);
828 }
829 else // from materials
830 {
831 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Mat_DCAxy_Pt"))->Fill(pt, dcaxy);
832 }
833 }
834 }
835 }
836 }
837
838 return nTracks;
839}
840
841void AliAnalysisTaskB2AOD::Terminate(Option_t* )
842{
843 // The Terminate() function is the last function to be called during
844 // a query. It always runs on the client, it can be used to present
845 // the results graphically or save the results to file.
846
847}
848
849Bool_t AliAnalysisTaskB2AOD::IsV0AND(UInt_t /*triggerBits*/) const
850{
851//
852// signals in both V0A and V0C
853//
854 //return ( (triggerBits&AliVEvent::kINT7) == AliVEvent::kINT7 );
855
856 AliAODVZERO* vzero = fAODevent->GetVZEROData();
857 if(vzero == 0) return kFALSE;
858
859 if(!vzero->TestBit(AliAODVZERO::kDecisionFilled)) return kFALSE;
860
861 return ( (vzero->GetV0ADecision() == AliVVZERO::kV0BB) && (vzero->GetV0CDecision() == AliVVZERO::kV0BB) );
862}
863
864Bool_t AliAnalysisTaskB2AOD::IsFastOnly(UInt_t triggerBits) const
865{
866//
867// kFastOnly trigger
868//
869 return ( (triggerBits&AliVEvent::kFastOnly) == AliVEvent::kFastOnly );
870}
871
872Bool_t AliAnalysisTaskB2AOD::IsMB(UInt_t triggerBits) const
873{
874//
875// MB event
876//
877 return ( (triggerBits&AliVEvent::kMB) == AliVEvent::kMB );
878}
879
880AliAODMCParticle* AliAnalysisTaskB2AOD::GetParticle(const AliAODTrack* trk) const
881{
882//
883// Particle that left the track
884//
885 TList* lst = fAODevent->GetList();
886 TClonesArray* lop = dynamic_cast<TClonesArray*>(lst->FindObject(AliAODMCParticle::StdBranchName()));
887
888 if (!lop) return 0;
889
890 Int_t label = TMath::Abs(trk->GetLabel()); // if negative then it shares points from other tracks
891 if( label >= lop->GetEntries() ) return 0;
892
893 return dynamic_cast<AliAODMCParticle*>(lop->At(label));
894}
895
896Bool_t AliAnalysisTaskB2AOD::IsFakeTrack(const AliAODTrack* trk) const
897{
898//
899// Check if the track shares some clusters with different particles
900// (definition changed to label=0? )
901//
902 return ( trk->GetLabel() < 0 );
903}
904
905Double_t AliAnalysisTaskB2AOD::GetPhi(Double_t p[3]) const
906{
907//
908// Azimuthal angle [0,2pi)
909//
910 Double_t px = p[0];
911 Double_t py = p[1];
912
913 return TMath::Pi()+TMath::ATan2(-py, -px);
914}
915
916Double_t AliAnalysisTaskB2AOD::GetBeta(const AliAODTrack* trk) const
917{
918//
919// Velocity
920//
921 Double_t t = trk->GetTOFsignal(); // ps
922 Double_t l = fTrackCuts->GetIntegratedLength(trk, fPartCode); // cm
923
924 if(t <= 0) return 1.e+6; // 1M times the speed of light ;)
925
926 return (l/t)/2.99792458e-2;
927}
928
929Double_t AliAnalysisTaskB2AOD::GetExpectedTime(const AliAODTrack* trk, Double_t m) const
930{
931//
932// Expected time workaround for the given mass hypothesis (ps)
933//
934 Double_t p = (fPartCode > AliPID::kTriton) ? 2.*trk->P() : trk->P();
935 Double_t beta = p/TMath::Sqrt(p*p + m*m);
936 Double_t l = fTrackCuts->GetIntegratedLength(trk, fPartCode);
937
938 return l/beta/2.99792458e-2;
939}
940
941Double_t AliAnalysisTaskB2AOD::GetMassSquared(Double_t p, Double_t beta) const
942{
943//
944// Mass squared
945//
946 return p*p*(1./(beta*beta) - 1.);
947}
948
949Double_t AliAnalysisTaskB2AOD::GetM2Difference(Double_t beta, Double_t p, Double_t m) const
950{
951//
952// Mass squared difference
953//
954 Double_t expBeta2 = p*p/(p*p+m*m);
955
956 return p*p*(1./(beta*beta)-1./expBeta2);
957}
958
959Int_t AliAnalysisTaskB2AOD::GetChargedMultiplicity(Double_t etaMax) const
960{
961//
962// Charged particle multiplicity using ALICE physical primary definition
963//
964 TList* lst = fAODevent->GetList();
965 TClonesArray* lop = dynamic_cast<TClonesArray*>(lst->FindObject(AliAODMCParticle::StdBranchName()));
966
967 if (!lop) return 0;
968
969 Int_t nch = 0;
970
971 for (Int_t i = 0; i < lop->GetEntries(); ++i)
972 {
973 AliAODMCParticle* iParticle = dynamic_cast<AliAODMCParticle*>(lop->At(i));
974
975 if(!iParticle) continue;
976
977 if(!iParticle->IsPhysicalPrimary()) continue;
978
979 if(TMath::Abs(iParticle->Eta()) >= etaMax) continue;
980
981 if(iParticle->Charge()==0) continue;
982
983 ++nch;
984 }
985
986 return nch;
987}
988
989Int_t AliAnalysisTaskB2AOD::GetITSnClusters(const AliAODTrack* trk) const
990{
991//
992// ITS number of clusters
993//
994 UChar_t map = trk->GetITSClusterMap();
995
996 Int_t npoints=0;
997 for(Int_t j=0; j<6; j++)
998 {
999 if(map&(1<<j)) ++npoints;
1000 }
1001
1002 return npoints;
1003}
1004
1005Int_t AliAnalysisTaskB2AOD::GetITSnPointsPID(const AliAODTrack* trk) const
1006{
1007//
1008// ITS number of points for PID
1009//
1010 UChar_t map = trk->GetITSClusterMap();
1011
1012 Int_t npoints = 0;
1013 for(Int_t j=2; j<6; j++)
1014 {
1015 if(map&(1<<j)) ++npoints;
1016 }
1017
1018 return npoints;
1019}
1020
1021Int_t AliAnalysisTaskB2AOD::GetPidCode(const TString& species) const
1022{
1023//
1024// Return AliPID code of the given species
1025//
1026 TString name = species;
1027 name.ToLower();
1028
1029 if(name == "electron") return AliPID::kElectron;
1030 if(name == "muon") return AliPID::kMuon;
1031 if(name == "pion") return AliPID::kPion;
1032 if(name == "kaon") return AliPID::kKaon;
1033 if(name == "proton") return AliPID::kProton;
1034 if(name == "deuteron") return AliPID::kDeuteron;
1035 if(name == "triton") return AliPID::kTriton;
1036 if(name == "he3") return AliPID::kHe3;
1037 if(name == "alpha") return AliPID::kAlpha;
1038
1039 return -1;
1040}
1041
1042Double_t AliAnalysisTaskB2AOD::GetMomentumCorrection(Double_t ptrec) const
1043{
1044//
1045// momentum correction for low pt
1046//
1047 if(fMoCpfx == 0) return 0;
1048
1049 return fMoCpfx->Interpolate(ptrec);
1050}
1051
1052Double_t AliAnalysisTaskB2AOD::GetRapidity(Double_t p, Double_t pz, Double_t m) const
1053{
1054//
1055// Rapidity
1056//
1057 Double_t e = TMath::Sqrt(p*p + m*m);
1058
1059 if(e <= pz) return 1.e+16;
1060
1061 return 0.5*TMath::Log( (e+pz)/(e-pz) );
1062}