]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
option for TPC crossed rows
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / AliAnalysisTaskB2.cxx
CommitLineData
6b91e8c0 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
17// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
18
19#include <Riostream.h>
20#include <TTree.h>
21#include <TChain.h>
22#include <TString.h>
23
24#include <AliLog.h>
25
26#include <AliAnalysisTask.h>
27#include <AliAnalysisManager.h>
28
29#include <AliESD.h>
30#include <AliESDEvent.h>
31#include <AliESDInputHandler.h>
32#include <AliESDtrackCuts.h>
33#include <AliExternalTrackParam.h>
34
35#include <AliMCEvent.h>
36#include <AliMCEventHandler.h>
37#include <AliMCVertex.h>
38#include <AliStack.h>
39#include <TParticle.h>
40#include <TParticlePDG.h>
41#include <TMCProcess.h>
42
43#include <AliTriggerAnalysis.h> // for offline signals
44#include <AliCentrality.h>
45#include <AliESDUtils.h>
46
47#include <AliESDpid.h>
48
49#include <TH1D.h>
50#include <TH2D.h>
68a84436 51#include <TList.h>
2b4e5f2c 52#include <TProfile.h>
6b91e8c0 53
54#include "AliLnID.h"
55#include "AliLnHistoMap.h"
56#include "AliAnalysisTaskB2.h"
57
58ClassImp(AliAnalysisTaskB2)
59
60AliAnalysisTaskB2::AliAnalysisTaskB2()
61: AliAnalysisTask()
62, fSpecies("Proton")
63, fPartCode(AliPID::kProton)
6b91e8c0 64, fHeavyIons(0)
65, fSimulation(0)
dc8b5f72 66, fMultTriggerFired(0)
67, fCentTriggerFired(0)
6b91e8c0 68, fTriggerFired(0)
69, fGoodVertex(0)
70, fPileUpEvent(0)
6b91e8c0 71, fV0AND(0)
72, fNoFastOnly(0)
dc8b5f72 73, fNtrkMultTrigger(0)
6b91e8c0 74, fMinKNOmult(-10)
75, fMaxKNOmult(100000)
76, fMinCentrality(0)
77, fMaxCentrality(100)
ec499c6d 78, fNch(0)
6b91e8c0 79, fNtrk(0)
80, fMeanNtrk(1)
81, fKNOmult(1)
6b91e8c0 82, fMinVx(-1)
83, fMaxVx(1)
84, fMinVy(-1)
85, fMaxVy(1)
86, fMinVz(-10)
87, fMaxVz(10)
6b91e8c0 88, fMinDCAxy(-1)
89, fMaxDCAxy(1)
90, fMinDCAz(-2)
91, fMaxDCAz(2)
92, fMaxNSigma(3)
6b91e8c0 93, fMinEta(-0.8)
94, fMaxEta(0.8)
95, fMinY(-0.5)
96, fMaxY(0.5)
6b91e8c0 97, fMCevent(0)
98, fESDevent(0)
68a84436 99, fOutputContainer(0)
6b91e8c0 100, fHistoMap(0)
101, fESDtrackCuts(0)
102, fLnID(0)
6b91e8c0 103, fMaxNSigmaITS(3)
104, fMaxNSigmaTPC(3)
105, fMaxNSigmaTOF(3)
6b91e8c0 106, fTrigAna(0)
107, fESDpid(0)
108, fIsPidOwner(0)
109, fTimeZeroType(AliESDpid::kTOF_T0)
110, fTOFmatch(0)
1cb68411 111, fMinM2(2.)
112, fMaxM2(6.)
2b4e5f2c 113, fMomentumCorrection(0)
114, fMoCpfx(0)
6b91e8c0 115
116{
117//
118// Default constructor
119//
120 AliLog::SetGlobalLogLevel(AliLog::kFatal);
121
122 fTrigAna = new AliTriggerAnalysis();
123}
124
125AliAnalysisTaskB2::AliAnalysisTaskB2(const char* name)
126: AliAnalysisTask(name,"")
127, fSpecies("Proton")
128, fPartCode(AliPID::kProton)
6b91e8c0 129, fHeavyIons(0)
130, fSimulation(0)
dc8b5f72 131, fMultTriggerFired(0)
132, fCentTriggerFired(0)
6b91e8c0 133, fTriggerFired(0)
134, fGoodVertex(0)
135, fPileUpEvent(0)
6b91e8c0 136, fV0AND(0)
137, fNoFastOnly(0)
dc8b5f72 138, fNtrkMultTrigger(0)
6b91e8c0 139, fMinKNOmult(-10)
140, fMaxKNOmult(100000)
141, fMinCentrality(-1)
142, fMaxCentrality(100)
b9e8e52e 143, fNch(0)
6b91e8c0 144, fNtrk(0)
145, fMeanNtrk(1)
146, fKNOmult(1)
6b91e8c0 147, fMinVx(-1)
148, fMaxVx(1)
149, fMinVy(-1)
150, fMaxVy(1)
151, fMinVz(-10)
152, fMaxVz(10)
6b91e8c0 153, fMinDCAxy(-1)
154, fMaxDCAxy(1)
155, fMinDCAz(-2)
156, fMaxDCAz(2)
157, fMaxNSigma(3)
6b91e8c0 158, fMinEta(-0.8)
159, fMaxEta(0.8)
160, fMinY(-0.5)
161, fMaxY(0.5)
6b91e8c0 162, fMCevent(0)
163, fESDevent(0)
68a84436 164, fOutputContainer(0)
6b91e8c0 165, fHistoMap(0)
166, fESDtrackCuts(0)
167, fLnID(0)
6b91e8c0 168, fMaxNSigmaITS(3)
169, fMaxNSigmaTPC(3)
170, fMaxNSigmaTOF(3)
6b91e8c0 171, fTrigAna(0)
172, fESDpid(0)
173, fIsPidOwner(0)
174, fTimeZeroType(AliESDpid::kTOF_T0)
175, fTOFmatch(0)
1cb68411 176, fMinM2(2.)
177, fMaxM2(6.)
2b4e5f2c 178, fMomentumCorrection(0)
179, fMoCpfx(0)
6b91e8c0 180
181{
182//
183// Constructor
184//
185 DefineInput(0, TChain::Class());
b2c24401 186 DefineOutput(0, TTree::Class());
187 DefineOutput(1, TList::Class());
6b91e8c0 188
189 //kFatal, kError, kWarning, kInfo, kDebug, kMaxType
190 AliLog::SetGlobalLogLevel(AliLog::kFatal);
191
192 fTrigAna = new AliTriggerAnalysis();
193}
194
195void AliAnalysisTaskB2::ConnectInputData(Option_t *)
196{
197//
198// Connect input data
199//
6b91e8c0 200 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
201 if(!tree)
202 {
f20bd9a4 203 this->Error("ConnectInputData", "could not read chain from input slot 0");
6b91e8c0 204 return;
205 }
206
207 // Get the pointer to the existing analysis manager via the static access method.
208 AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
209 if (!mgr)
210 {
f20bd9a4 211 this->Error("ConnectInputData", "could not get analysis manager");
6b91e8c0 212 return;
213 }
214
215 // cast type AliVEventHandler
216 AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(mgr->GetInputEventHandler());
217 if (!esdH)
218 {
f20bd9a4 219 this->Error("ConnectInputData", "could not get ESDInputHandler");
6b91e8c0 220 return;
221 }
222
223 // Get pointer to esd event from input handler
224 fESDevent = esdH->GetEvent();
225
226 // PID object for TOF
227 fESDpid = esdH->GetESDpid();
228
229 if(!fESDpid)
230 { //in case of no Tender attached
231 fESDpid = new AliESDpid();
232 fIsPidOwner = kTRUE;
233 }
6b91e8c0 234}
235
236void AliAnalysisTaskB2::CreateOutputObjects()
237{
238//
239// Create output objects
240//
241 if(fHistoMap == 0) exit(1); // should be created somewhere else
242
68a84436 243 fOutputContainer = new TList();
05798e84 244 fOutputContainer->SetOwner(kTRUE);
68a84436 245
246 TObjString* key;
247 TIter iter(fHistoMap->GetMap());
248 while((key = (TObjString*)iter.Next())) fOutputContainer->Add((TH1*)fHistoMap->Get(key));
249
b2c24401 250 PostData(1, fOutputContainer);
6b91e8c0 251}
252
253AliAnalysisTaskB2::~AliAnalysisTaskB2()
254{
255//
256// Default destructor
257//
258 delete fHistoMap;
68a84436 259 delete fOutputContainer;
6b91e8c0 260 delete fESDtrackCuts;
261 delete fLnID;
262
263 delete fTrigAna;
264
265 if(fIsPidOwner) delete fESDpid;
2b4e5f2c 266
267 delete fMoCpfx;
6b91e8c0 268}
269
270void AliAnalysisTaskB2::SetParticleSpecies(const TString& species)
271{
272//
273// set the particle species and the AliPID code
274//
275 fSpecies = species;
276 fPartCode = this->GetPidCode(species);
277}
278
279void AliAnalysisTaskB2::Exec(Option_t* )
280{
281//
282// Execute analysis for the current event
283//
f20bd9a4 284 if(fESDevent == 0)
6b91e8c0 285 {
f20bd9a4 286 this->Error("Exec", "AliESDEvent not available");
6b91e8c0 287 return;
288 }
289
290 if(fESDtrackCuts == 0)
291 {
f20bd9a4 292 this->Error("Exec", "ESD track cuts not set");
6b91e8c0 293 return;
294 }
295
296 if(fLnID == 0)
297 {
f20bd9a4 298 this->Error("Exec", "PID not set");
6b91e8c0 299 return;
300 }
301
87ed2a3a 302 if(fSimulation)
303 {
304 AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
305
306 if(mcH == 0) return;
307
308 fMCevent = mcH->MCEvent();
309
310 if(fMCevent == 0) return;
311 }
312
6b91e8c0 313 // --------- multiplicity and centrality ------------------
314
68a84436 315 fNtrk = AliESDtrackCuts::GetReferenceMultiplicity(fESDevent, AliESDtrackCuts::kTrackletsITSTPC, fMaxEta);
316 if(fSimulation) fNch = this->GetChargedMultiplicity(fMaxEta);
6b91e8c0 317
6b91e8c0 318 fKNOmult = fNtrk/fMeanNtrk;
319
1a9524cf 320 ((TH1D*)fHistoMap->Get(fSpecies + "_Event_Ntrk"))->Fill(fNtrk);
321 ((TH1D*)fHistoMap->Get(fSpecies + "_Event_Zmult"))->Fill(fKNOmult);
322
dc8b5f72 323 fMultTriggerFired = (fNtrk > 0) && (fKNOmult >= fMinKNOmult) && (fKNOmult < fMaxKNOmult);
6b91e8c0 324
325 if(fHeavyIons)
326 {
327 AliCentrality* esdCent = fESDevent->GetCentrality();
328
329 Float_t centrality = esdCent->GetCentralityPercentile("V0M");
330
dc8b5f72 331 fCentTriggerFired = (centrality >= fMinCentrality) && (centrality < fMaxCentrality);
6b91e8c0 332
333 Float_t v0ScaMult;
334 Float_t v0Mult = AliESDUtils::GetCorrV0(fESDevent,v0ScaMult);
335
336 ((TH1D*)fHistoMap->Get(fSpecies + "_V0_Mult"))->Fill(v0Mult);
337 ((TH1D*)fHistoMap->Get(fSpecies + "_V0_Scaled_Mult"))->Fill(v0ScaMult);
338 }
339
340 // ----------------- trigger --------------------
341
2e693567 342 AliInputEventHandler* eventH = dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
343 if(eventH == 0)
344 {
f20bd9a4 345 this->Error("Exec", "could not get AliInputEventHandler");
2e693567 346 return;
347 }
348
349 UInt_t triggerBits = eventH->IsEventSelected();
6b91e8c0 350
351 if(fHeavyIons)
352 {
dc8b5f72 353 fTriggerFired = ( this->IsMB(triggerBits) && fCentTriggerFired );
6b91e8c0 354 }
355 else
356 {
357 fTriggerFired = this->IsMB(triggerBits);
dc8b5f72 358 if(fNoFastOnly) fTriggerFired = ( fTriggerFired && !this->IsFastOnly(triggerBits) );
359 if(fV0AND) fTriggerFired = ( fTriggerFired && this->IsV0AND() );
360 if(fNtrkMultTrigger) fTriggerFired = ( fTriggerFired && fMultTriggerFired );
6b91e8c0 361 }
362
363 // --------------------- vertex --------------
364
dc8b5f72 365 fGoodVertex = kFALSE;
366
6b91e8c0 367 const AliESDVertex* vtx = fESDevent->GetPrimaryVertex(); // best primary vertex
368
369 if(vtx->GetStatus())
370 {
371 fGoodVertex=kTRUE;
372
373 ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_YX"))->Fill(vtx->GetX(), vtx->GetY());
374 ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_YZ"))->Fill(vtx->GetZ(), vtx->GetY());
375 ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_XZ"))->Fill(vtx->GetZ(), vtx->GetX());
376 }
377
378 if(fESDevent->GetPrimaryVertex()->IsFromVertexerZ())
379 {
380 if(fESDevent->GetPrimaryVertex()->GetDispersion() > 0.02) fGoodVertex=kFALSE;
381 if(fESDevent->GetPrimaryVertex()->GetZRes() > 0.25 ) fGoodVertex=kFALSE;
382 }
383
384 if( (vtx->GetX() <= fMinVx) || (vtx->GetX() >= fMaxVx) ) fGoodVertex=kFALSE;
385 if( (vtx->GetY() <= fMinVy) || (vtx->GetY() >= fMaxVy) ) fGoodVertex=kFALSE;
386 if( (vtx->GetZ() <= fMinVz) || (vtx->GetZ() >= fMaxVz) ) fGoodVertex=kFALSE;
387
388 // -------------------- pile up ------------------
389
dc8b5f72 390 fPileUpEvent = fESDevent->IsPileupFromSPDInMultBins();
6b91e8c0 391
392 // ---------------- event stats ----------------
393
394 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(0); // number of events
395
396 if(fTriggerFired)
397 {
f4d6dd11 398 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(1); // triggering events
6b91e8c0 399
400 if(vtx->GetStatus())
401 {
f4d6dd11 402 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(3); // valid vertex within a triggering event
6b91e8c0 403
404 if((vtx->GetZ() > fMinVz) && (vtx->GetZ() < fMaxVz))
405 {
406 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(4); // Vz
407
408 if( (vtx->GetX() > fMinVx) && (vtx->GetX() < fMaxVx)
409 && (vtx->GetY() > fMinVy) && (vtx->GetY() < fMaxVy))
410 {
411 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(5); // Vx, Vy
412
413 if(fGoodVertex)
414 {
415 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(6); // VertexerZ
416
417 if(!fPileUpEvent)
418 {
419 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(7); // no pile-up
420 }
421 }
422 }
423 }
424 }
425 }
426
427 // ------------- Particles and tracks for this event ---------------
428
429 Int_t nParticles = 0;
430 if(fSimulation)
431 {
432 nParticles = this->GetParticles();
433 }
434
435 this->GetTracks();
436
437 // Post the data (input/output slots #0 already used by the base class)
b2c24401 438 PostData(1, fOutputContainer);
6b91e8c0 439}
440
441Int_t AliAnalysisTaskB2::GetParticles()
442{
443//
444// Get particles from current event
445//
77dac0a6 446 Int_t nParticles = 0;
6b91e8c0 447
448 AliStack* stack = fMCevent->Stack();
449 if (!stack)
450 {
451 AliDebug(AliLog::kWarning, "stack not available");
452 return 0;
453 }
454
455 for (Int_t i = 0; i < fMCevent->GetNumberOfTracks(); ++i)
456 {
457 TParticle* iParticle = stack->Particle(i);
458
459 if(!iParticle) continue;
460
461 Int_t pid = fLnID->GetPID(iParticle);
462
463 if(pid != fPartCode) continue;
464
465 // ------- physical primaries ------------
466
467 if(!stack->IsPhysicalPrimary(i)) continue;
468
469 TString particle = fSpecies;
470 if(iParticle->GetPDG()->Charge() < 0) particle.Prepend("Anti");
471
472 Double_t genP = iParticle->P();
473 Double_t genPt = iParticle->Pt();
474 Double_t genY = iParticle->Y();
475 Double_t genPhi = iParticle->Phi();
476 Double_t genEta = iParticle->Eta();
477
ec499c6d 478 // ------- multiplicity and centrality -------------
479
dc8b5f72 480 if( fHeavyIons && !fCentTriggerFired) continue;
481 if( fNtrkMultTrigger && !fMultTriggerFired) continue;
6b91e8c0 482
483 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_P"))->Fill(genP);
484 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Pt"))->Fill(genPt);
485 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Y"))->Fill(genY);
486 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Phi"))->Fill(genPhi);
487 ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Eta"))->Fill(genEta);
488 ((TH2D*)fHistoMap->Get(particle + "_Gen_Prim_PtY"))->Fill(genY, genPt);
489 ((TH2D*)fHistoMap->Get(particle + "_Gen_Prim_EtaY"))->Fill(genY, genEta);
490
6b91e8c0 491 // ------ is within phase space? ----------
492
493 if(TMath::Abs(genY) >= fMaxY) continue;
494
ec499c6d 495 ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_P"))->Fill(genP);
6b91e8c0 496 ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_Pt"))->Fill(genPt);
497
f4d6dd11 498 // ------- is from a triggering event? (same as rec.) --------
6b91e8c0 499
500 if(!fTriggerFired)
501 {
502 ((TH1D*)fHistoMap->Get(particle + "_Gen_NoTrig_Prim_Pt"))->Fill(genPt);
503 continue;
504 }
505
506 ((TH1D*)fHistoMap->Get(particle + "_Gen_Trig_Prim_Pt"))->Fill(genPt);
507
f4d6dd11 508 // ------- is from a triggering event with good vertex? -------------
6b91e8c0 509
510 if(!fESDevent->GetPrimaryVertex()->GetStatus())
511 {
512 ((TH1D*)fHistoMap->Get(particle + "_Gen_NoVtx_Prim_Pt"))->Fill(genPt);
513 }
514
515 if(!fGoodVertex) continue;
516 if(fPileUpEvent) continue;
517
ec499c6d 518 ((TH1D*)fHistoMap->Get(particle + "_Gen_Nch"))->Fill(fNch);
6b91e8c0 519 ((TH1D*)fHistoMap->Get(particle + "_Gen_Vtx_Prim_Pt"))->Fill(genPt);
520
521 // ------ is within the geometrical acceptance? ------------
522
f4d6dd11 523 if(TMath::Abs(genEta) >= fMaxEta) continue;
6b91e8c0 524
525 ((TH2D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_PtY"))->Fill(genY,genPt);
526 ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_P"))->Fill(genP);
527 ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Pt"))->Fill(genPt);
528 ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Phi"))->Fill(genPhi);
529 ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Y"))->Fill(genY);
530 }
531
532 return nParticles;
533}
534
535Int_t AliAnalysisTaskB2::GetTracks()
536{
537//
538// Get tracks from current event
539//
540 using namespace std;
541
542 Int_t nTracks = 0;
543
544 // --------- trigger, vertex and pile-up -----------
545
546 if(!fTriggerFired) return 0;
547 if(!fGoodVertex) return 0;
548 if(fPileUpEvent) return 0;
549
550 // ------------ this is a 'good' event -------------
551
552 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(2); // analyzed events
553
554 ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_Event_Ntrk"))->Fill(fNtrk);
1a9524cf 555 ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_Event_Zmult"))->Fill(fKNOmult);
6b91e8c0 556
557 if(fSimulation)
558 {
ec499c6d 559 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Event_Nch_Ntrk"))->Fill(fNtrk, fNch);
6b91e8c0 560 }
561
562 const AliESDVertex* vtx = fESDevent->GetPrimaryVertex();
563
564 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_YX"))->Fill(vtx->GetX(), vtx->GetY());
565 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_YZ"))->Fill(vtx->GetZ(), vtx->GetY());
566 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_XZ"))->Fill(vtx->GetZ(), vtx->GetX());
567
568 if(fHeavyIons)
569 {
570 Float_t v0ScaMult;
571 Float_t v0Mult = AliESDUtils::GetCorrV0(fESDevent,v0ScaMult);
572
573 ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_V0_Mult"))->Fill(v0Mult);
574 ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_V0_Scaled_Mult"))->Fill(v0ScaMult);
575 }
576
577 // ------- track loop --------
578
579 for(Int_t i = 0; i < fESDevent->GetNumberOfTracks(); ++i)
580 {
581 AliESDtrack* iTrack = fESDevent->GetTrack(i);
582
583 if(!iTrack) continue;
584
ec499c6d 585 Float_t sign = 1;
586 TString particle = fSpecies;
587 if(iTrack->GetSign()<0 )
588 {
589 particle.Prepend("Anti");
590 sign = -1;
591 }
592
6b91e8c0 593 // -------- track cuts ------------------------
594
052278a6 595 Double_t eta = iTrack->Eta();
596 Double_t phi = this->GetPhi(iTrack);
6b91e8c0 597
052278a6 598 ((TH2D*)fHistoMap->Get(particle + "_Before_Phi_Eta"))->Fill(eta,phi);
6b91e8c0 599
600 if(!fESDtrackCuts->AcceptTrack(iTrack)) continue; // with next track
601 if(fTOFmatch && !this->AcceptTOFtrack(iTrack)) continue; // with next track
602
603 // --------------------- end track cuts -----------------------
604
052278a6 605 ((TH2D*)fHistoMap->Get(particle + "_After_Phi_Eta"))->Fill(eta, phi);
6b91e8c0 606
607 ++nTracks;
608
609 // charge
610
611 Double_t z = 1;
612 if(fPartCode>AliPID::kTriton) z = 2;
5b88d60e 613
6b91e8c0 614 // impact parameters
615
616 Float_t dcaxy, dcaz;
617 iTrack->GetImpactParameters(dcaxy, dcaz);
618
619 dcaxy *= sign;
620 dcaz *= sign;
621
622 Double_t nSigmaVtx = fESDtrackCuts->GetSigmaToVertex(iTrack);
623
624 // momentum
625
2b4e5f2c 626 Double_t p = iTrack->GetP()*z; // p at DCA
77dac0a6 627 Double_t pt = iTrack->Pt()*z; // pt at DCA
6b91e8c0 628 Double_t y = this->GetRapidity(iTrack, fPartCode);
629 Double_t pITS = this->GetITSmomentum(iTrack);
630 Double_t pTPC = iTrack->GetTPCmomentum();
631 Double_t pTOF = this->GetTOFmomentum(iTrack);
632 Double_t dEdxITS = iTrack->GetITSsignal();
633 Double_t dEdxTPC = iTrack->GetTPCsignal();
634 Int_t nPointsITS = this->GetITSnPointsPID(iTrack);
635 Int_t nPointsTPC = iTrack->GetTPCsignalN();
636
2b4e5f2c 637 if(fMomentumCorrection)
638 {
639 pt += this->GetMomentumCorrection(pt);
640 p = TMath::Sqrt(pt*pt + iTrack->Pz()*iTrack->Pz());
641 y = this->GetRapidity(p, iTrack->Pz(), fPartCode);
642 }
643
6b91e8c0 644 Double_t beta = 0;
645 Double_t mass = 0;
646 Double_t m2 = 0;
5b88d60e 647 Double_t dm2 = -100;
067aa49d 648 Double_t t = 0;
649 Double_t dt = -1000;
6b91e8c0 650
651 Double_t simPt = 0;
652 Double_t simPhi = 0;
653 Double_t simY = 0;
654
655 // --------- track cuts ------------
656
ec499c6d 657 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_DCAxy"))->Fill(dcaxy);
658 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_DCAz"))->Fill(dcaz);
659 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_NSigma"))->Fill(nSigmaVtx);
6b91e8c0 660
ec499c6d 661 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_ITSchi2PerCls"))->Fill(this->GetITSchi2PerCluster(iTrack));
6b91e8c0 662
ec499c6d 663 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCncls"))->Fill(iTrack->GetTPCNcls());
b3f10f7e 664 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCxRowsOverF"))->Fill(static_cast<Double_t>(iTrack->GetTPCCrossedRows())/static_cast<Double_t>(iTrack->GetTPCNclsF()));
ec499c6d 665 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCxRows"))->Fill(iTrack->GetTPCCrossedRows());
666 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCchi2PerCls"))->Fill(iTrack->GetTPCchi2()/iTrack->GetTPCNcls());
667 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCchi2Global"))->Fill(iTrack->GetChi2TPCConstrainedVsGlobal(fESDevent->GetPrimaryVertex()));
6b91e8c0 668
669 // -------------
670
ec499c6d 671 ((TH2D*)fHistoMap->Get(particle + "_ITS_dEdx_P"))->Fill(pITS, dEdxITS);
672 ((TH2D*)fHistoMap->Get(particle + "_TPC_dEdx_P"))->Fill(pTPC, dEdxTPC);
6b91e8c0 673
674 if(this->TOFmatch(iTrack))
675 {
676 beta = this->GetBeta(iTrack);
067aa49d 677 m2 = this->GetMassSquared(iTrack);
6b91e8c0 678 mass = TMath::Sqrt(TMath::Abs(m2));
067aa49d 679 dm2 = this->GetM2Difference(beta, pTOF, AliPID::ParticleMass(fPartCode));
680 t = this->GetTimeOfFlight(iTrack)*1.e-3; // ns
681 dt = t - 1.e-3*this->GetExpectedTime(iTrack, AliPID::ParticleMass(fPartCode));
6b91e8c0 682
ec499c6d 683 ((TH2D*)fHistoMap->Get(particle + "_TOF_Beta_P"))->Fill(pTOF, beta);
684 ((TH2D*)fHistoMap->Get(particle + "_TOF_Mass_P"))->Fill(pTOF, mass);
6b91e8c0 685 }
686
687 // -----------------------------------------------
688 // for pid and efficiency
689 // -----------------------------------------------
690
691 Int_t simpid = -1;
692 TParticle* iParticle = 0;
693
694 if(fSimulation)
695 {
696 iParticle = this->GetParticle(iTrack);
697
698 if(iParticle == 0) continue;
699
700 simPt = iParticle->Pt();
701 simPhi = iParticle->Phi();
702 simY = iParticle->Y();
703
704 simpid = fLnID->GetPID(iParticle);
705
706 if(simpid == fPartCode)
707 {
708 TString simparticle = fSpecies;
709 if(this->GetSign(iParticle)<0) simparticle.Prepend("Anti");
710
711 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_PtY"))->Fill(simY, simPt);
712
713 if(TMath::Abs(y) < fMaxY)
714 {
715 ((TH2D*)fHistoMap->Get(simparticle + "_Response_Matrix"))->Fill(pt, simPt);
716
ec499c6d 717 if(this->IsPhysicalPrimary(iParticle))
718 {
719 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Ntrk"))->Fill(fNtrk);
720 }
721
6b91e8c0 722 // for pid
723 if(!this->IsFakeTrack(iTrack))
724 {
725 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Pt"))->Fill(simPt);
726
727 if(this->IsPhysicalPrimary(iParticle)) // the efficiency is calculated on the primaries
728 {
729 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Pt"))->Fill(simPt);
730 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Y"))->Fill(simY);
731 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Phi"))->Fill(simPhi);
052278a6 732 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Rec_Pt"))->Fill(pt);
6b91e8c0 733
734 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_DCAxy_Pt"))->Fill(simPt,dcaxy);
735
736 ((TH2D*)fHistoMap->Get(simparticle + "_Prim_Response_Matrix"))->Fill(pt, simPt);
737
738 ((TH2D*)fHistoMap->Get(simparticle + "_Prim_DiffPt_RecPt"))->Fill(pt,simPt-pt);
739
740 if( this->TOFmatch(iTrack) )
741 {
742 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_P"))->Fill(pTOF, m2);
743 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_Pt"))->Fill(pt, m2);
744 }
1cb68411 745 }
746 else if(this->IsFromWeakDecay(iParticle))
747 {
748 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_Pt"))->Fill(simPt);
6b91e8c0 749
1cb68411 750 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_DCAxy_Pt"))->Fill(simPt,dcaxy);
6b91e8c0 751 }
1cb68411 752 else // if(this->IsFromMaterial(iParticle)
6b91e8c0 753 {
1cb68411 754 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Mat_Pt"))->Fill(simPt);
755
756 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Mat_DCAxy_Pt"))->Fill(simPt,dcaxy);
6b91e8c0 757 }
758 }
759 else // fake tracks
760 {
761 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Pt"))->Fill(simPt);
762 if(this->IsPhysicalPrimary(iParticle))
763 {
764 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Prim_Pt"))->Fill(simPt);
765 }
766 else if(this->IsFromWeakDecay(iParticle))
767 {
768 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Fdwn_Pt"))->Fill(simPt);
769 }
1cb68411 770 else
6b91e8c0 771 {
772 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Mat_Pt"))->Fill(simPt);
773 }
774 }
775 }
776 }
777 }
778
779 // get the pid
780 Int_t pid = fLnID->GetPID( fPartCode, pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta, fMaxNSigmaITS, fMaxNSigmaTPC, fMaxNSigmaTOF);
781
782 // pid table for prior probabilities (only Bayes)
783
784 Int_t offset = AliPID::kDeuteron - 5;
785
786 if( fSimulation )
787 {
788 Int_t sim = simpid > AliPID::kProton ? simpid - offset : simpid;
789 Int_t rec = pid > AliPID::kProton ? pid - offset : pid;
790
791 if((sim > -1) && (rec > -1)) // pid performance
792 {
793 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(sim, rec);
794 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(9, rec);
795 }
796 if(sim > -1)
797 {
798 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(sim, 9);
799 }
800 }
801
802 // for bayes iteration
803 if(pid != -1)
804 {
805 Int_t iBin = pid > AliPID::kProton ? pid - offset : pid;
806 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats_PID"))->Fill(iBin);
807 }
808
809 if(pid != fPartCode) continue;
810
811 Bool_t goodPid = 0;
812
813 if(fSimulation)
814 {
815 goodPid = ( simpid == pid );
816 }
817
067aa49d 818 ((TH2D*)fHistoMap->Get(particle + "_PID_Ntrk_pTPC"))->Fill(pTPC,fNtrk);
819 ((TH2D*)fHistoMap->Get(particle + "_PID_Zmult_pTPC"))->Fill(pTPC,fKNOmult);
ec499c6d 820
6b91e8c0 821 // pid performance
2e693567 822 ((TH2D*)fHistoMap->Get(particle + "_PID_ITSdEdx_P"))->Fill(pITS, dEdxITS);
77dac0a6 823 ((TH2D*)fHistoMap->Get(particle + "_PID_TPCdEdx_P"))->Fill(pTPC, dEdxTPC);
6b91e8c0 824
825 if(this->TOFmatch(iTrack))
826 {
77dac0a6 827 ((TH2D*)fHistoMap->Get(particle + "_PID_Beta_P"))->Fill(pTOF, beta);
828 ((TH2D*)fHistoMap->Get(particle + "_PID_Mass_P"))->Fill(pTOF, mass);
6b91e8c0 829 if(fSimulation && goodPid)
830 {
831 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Mass"))->Fill(mass);
832 }
833 }
834
835 // ---------------------------------
836 // results in |y| < fMaxY
837 // ---------------------------------
838
839 ((TH1D*)fHistoMap->Get(particle + "_PID_Y"))->Fill(y);
840 ((TH2D*)fHistoMap->Get(particle + "_PID_Pt_Y"))->Fill(y, pt);
841
842 if(TMath::Abs(y) >= fMaxY) continue;
843
844 ((TH1D*)fHistoMap->Get(particle + "_PID_Pt"))->Fill(pt);
845 ((TH1D*)fHistoMap->Get(particle + "_PID_Phi"))->Fill(phi);
846
7421dabb 847 if(iTrack->IsOn(AliESDtrack::kTRDin))
848 {
849 ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_Pt"))->Fill(pt);
850
851 if(iTrack->IsOn(AliESDtrack::kTRDout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_TRDout_Pt"))->Fill(pt);
852 if(iTrack->IsOn(AliESDtrack::kTOFout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_TOFout_Pt"))->Fill(pt);
853 }
854
855 if(iTrack->IsOn(AliESDtrack::kTOFin))
856 {
857 ((TH1D*)fHistoMap->Get(particle + "_PID_TOFin_Pt"))->Fill(pt);
858
859 if(iTrack->IsOn(AliESDtrack::kTOFout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TOFin_TOFout_Pt"))->Fill(pt);
860 }
861
6b91e8c0 862 if( this->TOFmatch(iTrack) )
863 {
864 ((TH2D*)fHistoMap->Get(particle + "_PID_M2_Pt"))->Fill(pt, m2);
5b88d60e 865 ((TH2D*)fHistoMap->Get(particle + "_PID_DM2_Pt"))->Fill(pt, dm2);
067aa49d 866 ((TH2D*)fHistoMap->Get(particle + "_PID_Time_Pt"))->Fill(pt, t);
867 ((TH2D*)fHistoMap->Get(particle + "_PID_DTime_Pt"))->Fill(pt, dt);
7421dabb 868 ((TH1D*)fHistoMap->Get(particle + "_PID_TOFmatch_Pt"))->Fill(pt);
6b91e8c0 869 }
870
871 // secondaries
1cb68411 872
873 Bool_t m2match = kTRUE;
874
875 if( fTOFmatch && (fLnID->GetPidProcedure() > AliLnID::kMaxLikelihood))
6b91e8c0 876 {
067aa49d 877 if( pt < 1.8 )
878 {
879 if ((m2 < fMinM2) || (m2 >= fMaxM2)) m2match = kFALSE;
880 }
881 else
882 {
883 if((dt < -0.6) || (dt >= 1)) m2match = kFALSE;
884 }
6b91e8c0 885 }
1cb68411 886
887 if(m2match)
6b91e8c0 888 {
889 ((TH2D*)fHistoMap->Get(particle + "_PID_DCAxy_Pt"))->Fill(pt, dcaxy);
890 ((TH2D*)fHistoMap->Get(particle + "_PID_DCAz_Pt"))->Fill(pt, dcaz);
891 ((TH2D*)fHistoMap->Get(particle + "_PID_NSigma_Pt"))->Fill(pt, nSigmaVtx);
892 }
893
894 if(fSimulation && goodPid)
895 {
896 // for unfolding and pid contamination
897 if(!this->IsFakeTrack(iTrack))
898 {
899 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Pt"))->Fill(simPt);
900
901 if( this->TOFmatch(iTrack) )
902 {
13dc5112 903 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_M2_Pt"))->Fill(pt, m2);
6b91e8c0 904 }
905
906 if(this->IsPhysicalPrimary(iParticle))
907 {
908 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Prim_Pt"))->Fill(simPt);
6b91e8c0 909 }
1cb68411 910
911 if(m2match)
6b91e8c0 912 {
1cb68411 913 if(this->IsPhysicalPrimary(iParticle))
914 {
915 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_DCAxy_Pt"))->Fill(pt, dcaxy);
916 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_DCAz_Pt"))->Fill(pt, dcaz);
917 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_NSigma_Pt"))->Fill(pt, nSigmaVtx);
918 }
919 else if(this->IsFromWeakDecay(iParticle))
6b91e8c0 920 {
13dc5112 921 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_DCAxy_Pt"))->Fill(pt, dcaxy);
922 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_DCAz_Pt"))->Fill(pt, dcaz);
923 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_NSigma_Pt"))->Fill(pt, nSigmaVtx);
6b91e8c0 924 }
1cb68411 925 else // from materials
6b91e8c0 926 {
13dc5112 927 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_DCAxy_Pt"))->Fill(pt, dcaxy);
928 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_DCAz_Pt"))->Fill(pt, dcaz);
929 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_NSigma_Pt"))->Fill(pt, nSigmaVtx);
6b91e8c0 930 }
931 }
932 }
933 else // fake tracks
934 {
935 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Pt"))->Fill(simPt);
936
1cb68411 937 if(m2match)
6b91e8c0 938 {
1cb68411 939 if(this->IsPhysicalPrimary(iParticle))
940 {
941 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Prim_DCAxy_Pt"))->Fill(pt, dcaxy);
942 }
943 else if(this->IsFromWeakDecay(iParticle))
944 {
945 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt"))->Fill(pt, dcaxy);
946 }
947 else // from materials
948 {
949 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Mat_DCAxy_Pt"))->Fill(pt, dcaxy);
950 }
6b91e8c0 951 }
952 }
953 }
954 }
955
956 return nTracks;
957}
958
959void AliAnalysisTaskB2::Terminate(Option_t* )
960{
961 // The Terminate() function is the last function to be called during
962 // a query. It always runs on the client, it can be used to present
963 // the results graphically or save the results to file.
964
965}
966
967Bool_t AliAnalysisTaskB2::IsV0AND() const
968{
969//
970// signals in both V0A and V0C
971//
972 return ( fTrigAna->IsOfflineTriggerFired(fESDevent, AliTriggerAnalysis::kV0A) &&
973 fTrigAna->IsOfflineTriggerFired(fESDevent, AliTriggerAnalysis::kV0C) );
974}
975
976Bool_t AliAnalysisTaskB2::IsFastOnly(UInt_t triggerBits) const
977{
978//
979// kFastOnly trigger
980//
981 return ( (triggerBits&AliVEvent::kFastOnly) == AliVEvent::kFastOnly );
982}
983
984Bool_t AliAnalysisTaskB2::IsMB(UInt_t triggerBits) const
985{
986//
987// MB event
988//
989 return ( (triggerBits&AliVEvent::kMB) == AliVEvent::kMB );
990}
991
992Bool_t AliAnalysisTaskB2::TOFmatch(const AliESDtrack* trk) const
993{
994//
995// Check for TOF match signal
996//
997 return ( trk->IsOn(AliESDtrack::kTOFout) && trk->IsOn(AliESDtrack::kTIME) );
998}
999
1000Bool_t AliAnalysisTaskB2::AcceptTOFtrack(const AliESDtrack* trk) const
1001{
1002//
1003// Additional checks for TOF match signal
1004//
1005 if( !this->TOFmatch(trk) ) return kFALSE;
1006 if( trk->GetIntegratedLength() < 350) return kFALSE;
1007 if( trk->GetTOFsignal() < 1e-6) return kFALSE;
1008
1009 return kTRUE;
1010}
1011
1012TParticle* AliAnalysisTaskB2::GetParticle(const AliESDtrack* trk) const
1013{
1014//
1015// Particle that left the track
1016//
1017 AliStack* stack = fMCevent->Stack();
1018
1019 Int_t label = TMath::Abs(trk->GetLabel()); // if negative then it shares points from other tracks
1020 if( label >= fMCevent->GetNumberOfTracks() ) return 0;
1021
1022 return stack->Particle(label);
1023}
1024
1025Bool_t AliAnalysisTaskB2::IsFakeTrack(const AliESDtrack* trk) const
1026{
1027//
1028// Check if the track shares some clusters with different particles
0f539a2b 1029// (definition changed to label=0? )
6b91e8c0 1030//
1031 return ( trk->GetLabel() < 0 );
1032}
1033
1034Bool_t AliAnalysisTaskB2::IsPhysicalPrimary(const TParticle* prt) const
1035{
1036//
1037// Check if the particle is physical primary
1038//
1039 AliStack* stack = fMCevent->Stack();
1040 Int_t index = stack->Particles()->IndexOf(prt);
1041
1042 return stack->IsPhysicalPrimary(index);
1043}
1044
1045Bool_t AliAnalysisTaskB2::IsFromMaterial(const TParticle* prt) const
1046{
1047//
1048// Check if the particle is originated at the materials
1049//
1050 AliStack* stack = fMCevent->Stack();
1051 Int_t index = stack->Particles()->IndexOf(prt);
1052
1053 return stack->IsSecondaryFromMaterial(index);
1054}
1055
1056Bool_t AliAnalysisTaskB2::IsFromWeakDecay(const TParticle* prt) const
1057{
1058//
1059// Check if the particle comes from a weak decay
1060//
1061 AliStack* stack = fMCevent->Stack();
1062 Int_t index = stack->Particles()->IndexOf(prt);
1063
1064 return stack->IsSecondaryFromWeakDecay(index);
1065}
1066
1067Double_t AliAnalysisTaskB2::GetSign(TParticle* prt) const
1068{
1069//
1070// Sign of the particle
1071//
1072 TParticlePDG* pdg = prt->GetPDG();
1073
1074 if(pdg != 0) return pdg->Charge();
1075
1076 return 0;
1077}
1078
1079Double_t AliAnalysisTaskB2::GetPhi(const AliESDtrack* trk) const
1080{
1081//
1082// Azimuthal angle [0,2pi) using the pt at the DCA point
1083//
1084 Double_t px = trk->Px();
1085 Double_t py = trk->Py();
1086
1087 return TMath::Pi()+TMath::ATan2(-py, -px);
1088}
1089
1090Double_t AliAnalysisTaskB2::GetTheta(const AliESDtrack* trk) const
1091{
1092//
1093// Polar angle using the pt at the DCA point
1094//
1095 Double_t p = trk->GetP();
1096 Double_t pz = trk->Pz();
1097
1098 return (pz == 0) ? TMath::PiOver2() : TMath::ACos(pz/p);
1099}
1100
1101Double_t AliAnalysisTaskB2::GetRapidity(const AliESDtrack* trk, Int_t pid) const
1102{
1103//
1104// Rapidity assuming the given particle hypothesis
1105// and using the momentum at the DCA
1106//
1107 Double_t m = AliPID::ParticleMass(pid);
77dac0a6 1108 Double_t p = (pid>AliPID::kTriton) ? 2.*trk->GetP() : trk->GetP();
6b91e8c0 1109 Double_t e = TMath::Sqrt(p*p + m*m);
77dac0a6 1110 Double_t pz = (pid>AliPID::kTriton) ? 2.*trk->Pz() : trk->Pz();
6b91e8c0 1111
1112 if(e <= pz) return 1.e+16;
1113
1114 return 0.5*TMath::Log( (e+pz)/(e-pz) );
1115}
1116
1117Double_t AliAnalysisTaskB2::GetITSmomentum(const AliESDtrack* trk) const
1118{
1119//
1120// Momentum for ITS pid
1121//
1122 Double_t pDCA = trk->GetP();
1123 Double_t pTPC = trk->GetTPCmomentum();
1124
1125 return (pDCA+pTPC)/2.;
1126}
1127
1128Double_t AliAnalysisTaskB2::GetTOFmomentum(const AliESDtrack* trk) const
1129{
1130//
1131// Momentum for TOF pid
1132//
067aa49d 1133 Double_t pIn = trk->GetTPCmomentum();
6b91e8c0 1134
1135 const AliExternalTrackParam* param = trk->GetOuterParam();
1136 Double_t pOut = param ? param->GetP() : trk->GetP();
1137
067aa49d 1138 return (pIn+pOut)/2.;
6b91e8c0 1139}
1140
1141Double_t AliAnalysisTaskB2::GetBeta(const AliESDtrack* trk) const
1142{
1143//
1144// Velocity
1145//
1146 Double_t t = this->GetTimeOfFlight(trk); // ps
1147 Double_t l = trk->GetIntegratedLength(); // cm
1148
1149 if(t <= 0) return 1.e6; // 1M times the speed of light ;)
1150
1151 return (l/t)/2.99792458e-2;
1152}
1153
067aa49d 1154Double_t AliAnalysisTaskB2::GetExpectedTime(const AliESDtrack* trk, Double_t m) const
6b91e8c0 1155{
1156//
067aa49d 1157// Expected time (ps) for the given mass hypothesis
1158//
1159 Double_t p = (fPartCode>AliPID::kTriton) ? 2.*this->GetTOFmomentum(trk) : this->GetTOFmomentum(trk);
1160 Double_t beta = p/TMath::Sqrt(p*p + m*m);
1161 Double_t l = trk->GetIntegratedLength();
1162
1163 return l/beta/2.99792458e-2;
1164}
1165
1166Double_t AliAnalysisTaskB2::GetMassSquared(const AliESDtrack* trk) const
1167{
1168//
1169// Mass squared
6b91e8c0 1170//
77dac0a6 1171 Double_t p = (fPartCode>AliPID::kTriton) ? 2.*this->GetTOFmomentum(trk) : this->GetTOFmomentum(trk);
6b91e8c0 1172 Double_t beta = this->GetBeta(trk);
1173
1174 return p*p*(1./(beta*beta) - 1.);
1175}
1176
067aa49d 1177Double_t AliAnalysisTaskB2::GetM2Difference(Double_t beta, Double_t momentum, Double_t m) const
1178{
1179//
1180// Mass squared difference
1181//
1182 Double_t p = (fPartCode>AliPID::kTriton) ? 2.*momentum : momentum;
1183 Double_t expBeta2 = p*p/(p*p+m*m);
1184
1185 return p*p*(1./(beta*beta)-1./expBeta2);
1186}
1187
77dac0a6 1188Int_t AliAnalysisTaskB2::GetChargedMultiplicity(Double_t etaMax) const
6b91e8c0 1189{
1190//
1191// Charged particle multiplicity using ALICE physical primary definition
1192//
1193 AliStack* stack = fMCevent->Stack();
1194
77dac0a6 1195 Int_t nch = 0;
6b91e8c0 1196 //for (Int_t i=0; i < stack->GetNprimary(); ++i)
1197 for (Int_t i=0; i < stack->GetNtrack(); ++i)
1198 {
1199 if(!stack->IsPhysicalPrimary(i)) continue;
1200
1201 TParticle* iParticle = stack->Particle(i);
1202
1203 if(TMath::Abs(iParticle->Eta()) >= etaMax) continue;
1204 //if (iParticle->Pt() < -1) continue;
1205
1206 TParticlePDG* iPDG = iParticle->GetPDG(); // There are some particles with no PDG
1207 if (iPDG && iPDG->Charge() == 0) continue;
1208
1209 ++nch;
1210 }
1211
1212 return nch;
1213}
1214
1215Double_t AliAnalysisTaskB2::GetTimeOfFlight(const AliESDtrack* trk) const
1216{
1217//
1218// Time of flight associated to the track.
1219// Adapted from ANALYSIS/AliAnalysisTaskESDfilter.cxx
1220//
1221 if(!fESDevent->GetTOFHeader())
1222 { //protection in case the pass2 LHC10b,c,d have been processed without tender.
1223 Float_t t0spread[10];
1224 Float_t intrinsicTOFres=100; //ps ok for LHC10b,c,d pass2!!
1225 for (Int_t j=0; j<10; j++) t0spread[j] = (TMath::Sqrt(fESDevent->GetSigma2DiamondZ()))/0.03; //0.03 to convert from cm to ps
1226 fESDpid->GetTOFResponse().SetT0resolution(t0spread);
1227 fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
1228
1229 fESDpid->SetTOFResponse(fESDevent, (AliESDpid::EStartTimeType_t)fTimeZeroType);
1230 }
1231
1232 if(fESDevent->GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(fESDevent, (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production strating form LHC10e without Tender.
1233
1234 Double_t timeZero = fESDpid->GetTOFResponse().GetStartTime(trk->P());
1235
1236 return trk->GetTOFsignal()-timeZero;
1237}
1238
1239Int_t AliAnalysisTaskB2::GetITSnClusters(const AliESDtrack* trk) const
1240{
1241//
1242// ITS number of clusters
1243//
1244 UChar_t map = trk->GetITSClusterMap();
1245
1246 Int_t npoints=0;
1247 for(Int_t j=0; j<6; j++)
1248 {
1249 if(map&(1<<j)) ++npoints;
1250 }
1251
1252 return npoints;
1253}
1254
1255Double_t AliAnalysisTaskB2::GetITSchi2PerCluster(const AliESDtrack* trk) const
1256{
1257//
1258// ITS chi2 per number of clusters
1259//
1260 Double_t chi2 = trk->GetITSchi2();
1261 Int_t ncls = this->GetITSnClusters(trk);
1262
1263 Double_t chi2ncls = (ncls==0) ? 1.e10 : chi2/ncls;
1264
1265 return chi2ncls;
1266}
1267
1268Int_t AliAnalysisTaskB2::GetITSnPointsPID(const AliESDtrack* trk) const
1269{
1270//
1271// ITS number of points for PID
1272//
1273 UChar_t map = trk->GetITSClusterMap();
1274
1275 Int_t npoints = 0;
1276 for(Int_t j=2; j<6; j++)
1277 {
1278 if(map&(1<<j)) ++npoints;
1279 }
1280
1281 return npoints;
1282}
1283
1284Int_t AliAnalysisTaskB2::GetPidCode(const TString& species) const
1285{
1286//
1287// Return AliPID code of the given species
1288//
1289 TString name = species;
1290 name.ToLower();
1291
1292 if(name == "electron") return AliPID::kElectron;
1293 if(name == "muon") return AliPID::kMuon;
1294 if(name == "pion") return AliPID::kPion;
1295 if(name == "kaon") return AliPID::kKaon;
1296 if(name == "proton") return AliPID::kProton;
1297 if(name == "deuteron") return AliPID::kDeuteron;
1298 if(name == "triton") return AliPID::kTriton;
1299 if(name == "he3") return AliPID::kHe3;
1300 if(name == "alpha") return AliPID::kAlpha;
1301
1302 return -1;
1303}
2b4e5f2c 1304
1305Double_t AliAnalysisTaskB2::GetMomentumCorrection(Double_t ptrec) const
1306{
1307//
1308// momentum correction for low pt
1309//
1310 if(fMoCpfx == 0) return 0;
1311
1312 return fMoCpfx->Interpolate(ptrec);
1313}
1314
1315Double_t AliAnalysisTaskB2::GetRapidity(Double_t p, Double_t pz, Int_t pid) const
1316{
1317//
1318// Rapidity (for momentum correction)
1319//
1320 Double_t m = AliPID::ParticleMass(pid);
1321 Double_t e = TMath::Sqrt(p*p + m*m);
1322
1323 if(e <= pz) return 1.e+16;
1324
1325 return 0.5*TMath::Log( (e+pz)/(e-pz) );
1326}