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