PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx

   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>
  51 #include <TList.h>
  52
  53 #include "AliLnID.h"
  54 #include "AliLnHistoMap.h"
  55 #include "AliAnalysisTaskB2.h"
  56
  57 ClassImp(AliAnalysisTaskB2)
  58
  59 AliAnalysisTaskB2::AliAnalysisTaskB2()
  60 : AliAnalysisTask()
  61 , fSpecies("Proton")
  62 , fPartCode(AliPID::kProton)
  63 , fHeavyIons(0)
  64 , fSimulation(0)
  65 , fMultTrigger(0)
  66 , fCentTrigger(0)
  67 , fTriggerFired(0)
  68 , fGoodVertex(0)
  69 , fPileUpEvent(0)
  70 , fV0AND(0)
  71 , fNoFastOnly(0)
  72 , fMinKNOmult(-10)
  73 , fMaxKNOmult(100000)
  74 , fMinCentrality(0)
  75 , fMaxCentrality(100)
  76 , fNch(0)
  77 , fNtrk(0)
  78 , fMeanNtrk(1)
  79 , fKNOmult(1)
  80 , fMinVx(-1)
  81 , fMaxVx(1)
  82 , fMinVy(-1)
  83 , fMaxVy(1)
  84 , fMinVz(-10)
  85 , fMaxVz(10)
  86 , fMinDCAxy(-1)
  87 , fMaxDCAxy(1)
  88 , fMinDCAz(-2)
  89 , fMaxDCAz(2)
  90 , fMaxNSigma(3)
  91 , fMinEta(-0.8)
  92 , fMaxEta(0.8)
  93 , fMinY(-0.5)
  94 , fMaxY(0.5)
  95 , fMCevent(0)
  96 , fESDevent(0)
  97 , fOutputContainer(0)
  98 , fHistoMap(0)
  99 , fESDtrackCuts(0)
 100 , fLnID(0)
 101 , fMaxNSigmaITS(3)
 102 , fMaxNSigmaTPC(3)
 103 , fMaxNSigmaTOF(3)
 104 , fTrigAna(0)
 105 , fESDpid(0)
 106 , fIsPidOwner(0)
 107 , fTimeZeroType(AliESDpid::kTOF_T0)
 108 , fTOFmatch(0)
 109 , fMinM2(2.)
 110 , fMaxM2(6.)
 111
 112 {
 113 //
 114 // Default constructor
 115 //
 116         AliLog::SetGlobalLogLevel(AliLog::kFatal);
 117
 118         fTrigAna = new AliTriggerAnalysis();
 119 }
 120
 121 AliAnalysisTaskB2::AliAnalysisTaskB2(const char* name)
 122 : AliAnalysisTask(name,"")
 123 , fSpecies("Proton")
 124 , fPartCode(AliPID::kProton)
 125 , fHeavyIons(0)
 126 , fSimulation(0)
 127 , fMultTrigger(0)
 128 , fCentTrigger(0)
 129 , fTriggerFired(0)
 130 , fGoodVertex(0)
 131 , fPileUpEvent(0)
 132 , fV0AND(0)
 133 , fNoFastOnly(0)
 134 , fMinKNOmult(-10)
 135 , fMaxKNOmult(100000)
 136 , fMinCentrality(-1)
 137 , fMaxCentrality(100)
 138 , fNch(0)
 139 , fNtrk(0)
 140 , fMeanNtrk(1)
 141 , fKNOmult(1)
 142 , fMinVx(-1)
 143 , fMaxVx(1)
 144 , fMinVy(-1)
 145 , fMaxVy(1)
 146 , fMinVz(-10)
 147 , fMaxVz(10)
 148 , fMinDCAxy(-1)
 149 , fMaxDCAxy(1)
 150 , fMinDCAz(-2)
 151 , fMaxDCAz(2)
 152 , fMaxNSigma(3)
 153 , fMinEta(-0.8)
 154 , fMaxEta(0.8)
 155 , fMinY(-0.5)
 156 , fMaxY(0.5)
 157 , fMCevent(0)
 158 , fESDevent(0)
 159 , fOutputContainer(0)
 160 , fHistoMap(0)
 161 , fESDtrackCuts(0)
 162 , fLnID(0)
 163 , fMaxNSigmaITS(3)
 164 , fMaxNSigmaTPC(3)
 165 , fMaxNSigmaTOF(3)
 166 , fTrigAna(0)
 167 , fESDpid(0)
 168 , fIsPidOwner(0)
 169 , fTimeZeroType(AliESDpid::kTOF_T0)
 170 , fTOFmatch(0)
 171 , fMinM2(2.)
 172 , fMaxM2(6.)
 173
 174 {
 175 //
 176 // Constructor
 177 //
 178         DefineInput(0, TChain::Class());
 179         DefineOutput(0, TTree::Class());
 180         DefineOutput(1, TList::Class());
 181
 182         //kFatal, kError, kWarning, kInfo, kDebug, kMaxType
 183         AliLog::SetGlobalLogLevel(AliLog::kFatal);
 184
 185         fTrigAna = new AliTriggerAnalysis();
 186 }
 187
 188 void AliAnalysisTaskB2::ConnectInputData(Option_t *)
 189 {
 190 //
 191 // Connect input data
 192 //
 193         TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 194         if(!tree)
 195         {
 196                 this->Error("ConnectInputData", "could not read chain from input slot 0");
 197                 return;
 198         }
 199
 200         // Get the pointer to the existing analysis manager via the static access method.
 201         AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
 202         if (!mgr)
 203         {
 204                 this->Error("ConnectInputData", "could not get analysis manager");
 205                 return;
 206         }
 207
 208         // cast type AliVEventHandler
 209         AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(mgr->GetInputEventHandler());
 210         if (!esdH)
 211         {
 212                 this->Error("ConnectInputData", "could not get ESDInputHandler");
 213                 return;
 214         }
 215
 216         // Get pointer to esd event from input handler
 217         fESDevent = esdH->GetEvent();
 218
 219         // PID object for TOF
 220         fESDpid = esdH->GetESDpid();
 221
 222         if(!fESDpid)
 223         { //in case of no Tender attached
 224                 fESDpid = new AliESDpid();
 225                 fIsPidOwner = kTRUE;
 226         }
 227
 228         if(!fSimulation) return;
 229
 230         AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (mgr->GetMCtruthEventHandler());
 231         if (!mcH)
 232         {
 233                 this->Error("ConnectInputData", "could not get AliMCEventHandler");
 234                 return;
 235         }
 236
 237         fMCevent = mcH->MCEvent();
 238         if (!fMCevent)
 239         {
 240                 this->Error("ConnectInputData", "could not get MC fLnEvent");
 241                 return;
 242         }
 243 }
 244
 245 void AliAnalysisTaskB2::CreateOutputObjects()
 246 {
 247 //
 248 // Create output objects
 249 //
 250         if(fHistoMap == 0) exit(1); // should be created somewhere else
 251
 252         fOutputContainer = new TList();
 253         fOutputContainer->SetOwner(kTRUE);
 254
 255         TObjString* key;
 256         TIter iter(fHistoMap->GetMap());
 257         while((key = (TObjString*)iter.Next())) fOutputContainer->Add((TH1*)fHistoMap->Get(key));
 258
 259         PostData(1, fOutputContainer);
 260 }
 261
 262 AliAnalysisTaskB2::~AliAnalysisTaskB2()
 263 {
 264 //
 265 // Default destructor
 266 //
 267         delete fHistoMap;
 268         delete fOutputContainer;
 269         delete fESDtrackCuts;
 270         delete fLnID;
 271
 272         delete fTrigAna;
 273
 274         if(fIsPidOwner) delete fESDpid;
 275 }
 276
 277 void AliAnalysisTaskB2::SetParticleSpecies(const TString& species)
 278 {
 279 //
 280 // set the particle species and the AliPID code
 281 //
 282         fSpecies = species;
 283         fPartCode = this->GetPidCode(species);
 284 }
 285
 286 void AliAnalysisTaskB2::Exec(Option_t* )
 287 {
 288 //
 289 // Execute analysis for the current event
 290 //
 291         if(fESDevent == 0)
 292         {
 293                 this->Error("Exec", "AliESDEvent not available");
 294                 return;
 295         }
 296
 297         if(fESDtrackCuts == 0)
 298         {
 299                 this->Error("Exec", "ESD track cuts not set");
 300                 return;
 301         }
 302
 303         if(fLnID == 0)
 304         {
 305                 this->Error("Exec", "PID not set");
 306                 return;
 307         }
 308
 309         fMultTrigger  = kFALSE;
 310         fCentTrigger  = kFALSE;
 311         fTriggerFired = kFALSE;
 312         fGoodVertex   = kFALSE;
 313         fPileUpEvent  = kFALSE;
 314
 315         // --------- multiplicity and centrality ------------------
 316
 317         fNtrk = AliESDtrackCuts::GetReferenceMultiplicity(fESDevent, AliESDtrackCuts::kTrackletsITSTPC, fMaxEta);
 318         if(fSimulation) fNch = this->GetChargedMultiplicity(fMaxEta);
 319
 320         ((TH1D*)fHistoMap->Get(fSpecies + "_Event_Ntrk"))->Fill(fNtrk);
 321
 322         fKNOmult = fNtrk/fMeanNtrk;
 323
 324         if( (fKNOmult >= fMinKNOmult) && (fKNOmult < fMaxKNOmult) ) fMultTrigger = kTRUE;
 325
 326         if(fHeavyIons)
 327         {
 328                 AliCentrality* esdCent = fESDevent->GetCentrality();
 329
 330                 Float_t centrality = esdCent->GetCentralityPercentile("V0M");
 331
 332                 if((centrality >= fMinCentrality) && (centrality < fMaxCentrality)) fCentTrigger = kTRUE;
 333
 334                 Float_t v0ScaMult;
 335                 Float_t v0Mult = AliESDUtils::GetCorrV0(fESDevent,v0ScaMult);
 336
 337                 ((TH1D*)fHistoMap->Get(fSpecies + "_V0_Mult"))->Fill(v0Mult);
 338                 ((TH1D*)fHistoMap->Get(fSpecies + "_V0_Scaled_Mult"))->Fill(v0ScaMult);
 339         }
 340
 341         // ----------------- trigger --------------------
 342
 343         AliInputEventHandler* eventH = dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 344         if(eventH == 0)
 345         {
 346                 this->Error("Exec", "could not get AliInputEventHandler");
 347                 return;
 348         }
 349
 350         UInt_t triggerBits = eventH->IsEventSelected();
 351
 352         if(fHeavyIons)
 353         {
 354                 fTriggerFired = ( this->IsMB(triggerBits) && fCentTrigger );
 355         }
 356         else
 357         {
 358                 fTriggerFired = this->IsMB(triggerBits);
 359                 if(fNoFastOnly) fTriggerFired = ( fTriggerFired && !this->IsFastOnly(triggerBits) );
 360                 if(fV0AND) fTriggerFired = ( fTriggerFired && this->IsV0AND() );
 361
 362                 fTriggerFired = ( fTriggerFired && fMultTrigger );
 363         }
 364
 365         // --------------------- vertex --------------
 366
 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
 390         if(fESDevent->IsPileupFromSPDInMultBins()) fPileUpEvent = kTRUE;
 391
 392         // ---------------- event stats ----------------
 393
 394         ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(0); // number of events
 395
 396         if(fTriggerFired)
 397         {
 398                 ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(1); // triggering events
 399
 400                 if(vtx->GetStatus())
 401                 {
 402                         ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(3); // valid vertex within a triggering event
 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)
 438         PostData(1, fOutputContainer);
 439 }
 440
 441 Int_t AliAnalysisTaskB2::GetParticles()
 442 {
 443 //
 444 // Get particles from current event
 445 //
 446         Int_t nParticles = 0;
 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
 478                 // ------- multiplicity and centrality -------------
 479
 480                 if( fHeavyIons && !fCentTrigger) continue;
 481                 if(!fHeavyIons && !fMultTrigger) continue;
 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
 491                 // ------ is within phase space? ----------
 492
 493                 if(TMath::Abs(genY) >= fMaxY) continue;
 494
 495                 ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_P"))->Fill(genP);
 496                 ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_Pt"))->Fill(genPt);
 497
 498                 // ------- is from a triggering event? (same as rec.) --------
 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
 508                 // ------- is from a triggering event with good vertex? -------------
 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
 518                 ((TH1D*)fHistoMap->Get(particle + "_Gen_Nch"))->Fill(fNch);
 519                 ((TH1D*)fHistoMap->Get(particle + "_Gen_Vtx_Prim_Pt"))->Fill(genPt);
 520
 521                 // ------ is within the geometrical acceptance? ------------
 522
 523                 if(TMath::Abs(genEta) >= fMaxEta) continue;
 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
 535 Int_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);
 555
 556         if(fSimulation)
 557         {
 558                 ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Event_Nch_Ntrk"))->Fill(fNtrk, fNch);
 559         }
 560
 561         const AliESDVertex* vtx = fESDevent->GetPrimaryVertex();
 562
 563         ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_YX"))->Fill(vtx->GetX(), vtx->GetY());
 564         ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_YZ"))->Fill(vtx->GetZ(), vtx->GetY());
 565         ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_XZ"))->Fill(vtx->GetZ(), vtx->GetX());
 566
 567         if(fHeavyIons)
 568         {
 569                 Float_t v0ScaMult;
 570                 Float_t v0Mult = AliESDUtils::GetCorrV0(fESDevent,v0ScaMult);
 571
 572                 ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_V0_Mult"))->Fill(v0Mult);
 573                 ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_V0_Scaled_Mult"))->Fill(v0ScaMult);
 574         }
 575
 576         // ------- track loop --------
 577
 578         for(Int_t i = 0; i < fESDevent->GetNumberOfTracks(); ++i)
 579         {
 580                 AliESDtrack* iTrack = fESDevent->GetTrack(i);
 581
 582                 if(!iTrack) continue;
 583
 584                 Float_t sign = 1;
 585                 TString particle = fSpecies;
 586                 if(iTrack->GetSign()<0 )
 587                 {
 588                         particle.Prepend("Anti");
 589                         sign = -1;
 590                 }
 591
 592                 // -------- track cuts ------------------------
 593
 594                 Double_t eta = iTrack->Eta();
 595                 Double_t phi = this->GetPhi(iTrack);
 596
 597                 ((TH2D*)fHistoMap->Get(particle + "_Before_Phi_Eta"))->Fill(eta,phi);
 598
 599                 if(!fESDtrackCuts->AcceptTrack(iTrack)) continue;  // with next track
 600                 if(fTOFmatch && !this->AcceptTOFtrack(iTrack)) continue; // with next track
 601
 602                 // --------------------- end track cuts -----------------------
 603
 604                 ((TH2D*)fHistoMap->Get(particle + "_After_Phi_Eta"))->Fill(eta, phi);
 605
 606                 ++nTracks;
 607
 608                 // charge
 609
 610                 Double_t z = 1;
 611                 if(fPartCode>AliPID::kTriton)  z = 2;
 612                 // impact parameters
 613
 614                 Float_t dcaxy, dcaz;
 615                 iTrack->GetImpactParameters(dcaxy, dcaz);
 616
 617                 dcaxy *= sign;
 618                 dcaz  *= sign;
 619
 620                 Double_t nSigmaVtx = fESDtrackCuts->GetSigmaToVertex(iTrack);
 621
 622                 // momentum
 623
 624                 Double_t pt      = iTrack->Pt()*z; // pt at DCA
 625                 Double_t y       = this->GetRapidity(iTrack, fPartCode);
 626                 Double_t pITS    = this->GetITSmomentum(iTrack);
 627                 Double_t pTPC    = iTrack->GetTPCmomentum();
 628                 Double_t pTOF    = this->GetTOFmomentum(iTrack);
 629                 Double_t dEdxITS = iTrack->GetITSsignal();
 630                 Double_t dEdxTPC = iTrack->GetTPCsignal();
 631                 Int_t nPointsITS = this->GetITSnPointsPID(iTrack);
 632                 Int_t nPointsTPC = iTrack->GetTPCsignalN();
 633
 634                 Double_t beta = 0;
 635                 Double_t mass = 0;
 636                 Double_t m2   = 0;
 637
 638                 Double_t simPt  = 0;
 639                 Double_t simPhi = 0;
 640                 Double_t simY   = 0;
 641
 642                 // --------- track cuts ------------
 643
 644                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_DCAxy"))->Fill(dcaxy);
 645                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_DCAz"))->Fill(dcaz);
 646                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_NSigma"))->Fill(nSigmaVtx);
 647
 648                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_ITSchi2PerCls"))->Fill(this->GetITSchi2PerCluster(iTrack));
 649
 650                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCncls"))->Fill(iTrack->GetTPCNcls());
 651                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCclsOverF"))->Fill((Double_t)iTrack->GetTPCNcls()/(Double_t)iTrack->GetTPCNclsF());
 652                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCxRows"))->Fill(iTrack->GetTPCCrossedRows());
 653                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCchi2PerCls"))->Fill(iTrack->GetTPCchi2()/iTrack->GetTPCNcls());
 654                 ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCchi2Global"))->Fill(iTrack->GetChi2TPCConstrainedVsGlobal(fESDevent->GetPrimaryVertex()));
 655
 656                 // -------------
 657
 658                 ((TH2D*)fHistoMap->Get(particle + "_ITS_dEdx_P"))->Fill(pITS, dEdxITS);
 659                 ((TH2D*)fHistoMap->Get(particle + "_TPC_dEdx_P"))->Fill(pTPC, dEdxTPC);
 660
 661                 if(this->TOFmatch(iTrack))
 662                 {
 663                         beta = this->GetBeta(iTrack);
 664                         m2   = this->GetMassSquare(iTrack);
 665                         mass = TMath::Sqrt(TMath::Abs(m2));
 666
 667                         ((TH2D*)fHistoMap->Get(particle + "_TOF_Beta_P"))->Fill(pTOF, beta);
 668                         ((TH2D*)fHistoMap->Get(particle + "_TOF_Mass_P"))->Fill(pTOF, mass);
 669                 }
 670
 671                 // -----------------------------------------------
 672                 //  for pid and efficiency
 673                 // -----------------------------------------------
 674
 675                 Int_t simpid = -1;
 676                 TParticle* iParticle = 0;
 677
 678                 if(fSimulation)
 679                 {
 680                         iParticle = this->GetParticle(iTrack);
 681
 682                         if(iParticle == 0) continue;
 683
 684                         simPt  = iParticle->Pt();
 685                         simPhi = iParticle->Phi();
 686                         simY   = iParticle->Y();
 687
 688                         simpid = fLnID->GetPID(iParticle);
 689
 690                         if(simpid == fPartCode)
 691                         {
 692                                 TString simparticle = fSpecies;
 693                                 if(this->GetSign(iParticle)<0) simparticle.Prepend("Anti");
 694
 695                                 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_PtY"))->Fill(simY, simPt);
 696
 697                                 if(TMath::Abs(y) < fMaxY)
 698                                 {
 699                                         ((TH2D*)fHistoMap->Get(simparticle + "_Response_Matrix"))->Fill(pt, simPt);
 700
 701                                         if(this->IsPhysicalPrimary(iParticle))
 702                                         {
 703                                                 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Ntrk"))->Fill(fNtrk);
 704                                         }
 705
 706                                         // for pid
 707                                         if(!this->IsFakeTrack(iTrack))
 708                                         {
 709                                                 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Pt"))->Fill(simPt);
 710
 711                                                 if(this->IsPhysicalPrimary(iParticle)) // the efficiency is calculated on the primaries
 712                                                 {
 713                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Pt"))->Fill(simPt);
 714                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Y"))->Fill(simY);
 715                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Phi"))->Fill(simPhi);
 716                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Rec_Pt"))->Fill(pt);
 717
 718                                                         ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_DCAxy_Pt"))->Fill(simPt,dcaxy);
 719
 720                                                         ((TH2D*)fHistoMap->Get(simparticle + "_Prim_Response_Matrix"))->Fill(pt, simPt);
 721
 722                                                         ((TH2D*)fHistoMap->Get(simparticle + "_Prim_DiffPt_RecPt"))->Fill(pt,simPt-pt);
 723
 724                                                         if( this->TOFmatch(iTrack) )
 725                                                         {
 726                                                                 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_P"))->Fill(pTOF, m2);
 727                                                                 ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_Pt"))->Fill(pt, m2);
 728                                                         }
 729                                                 }
 730                                                 else if(this->IsFromWeakDecay(iParticle))
 731                                                 {
 732                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_Pt"))->Fill(simPt);
 733
 734                                                         ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_DCAxy_Pt"))->Fill(simPt,dcaxy);
 735                                                 }
 736                                                 else // if(this->IsFromMaterial(iParticle)
 737                                                 {
 738                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Mat_Pt"))->Fill(simPt);
 739
 740                                                         ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Mat_DCAxy_Pt"))->Fill(simPt,dcaxy);
 741                                                 }
 742                                         }
 743                                         else // fake tracks
 744                                         {
 745                                                 ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Pt"))->Fill(simPt);
 746                                                 if(this->IsPhysicalPrimary(iParticle))
 747                                                 {
 748                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Prim_Pt"))->Fill(simPt);
 749                                                 }
 750                                                 else if(this->IsFromWeakDecay(iParticle))
 751                                                 {
 752                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Fdwn_Pt"))->Fill(simPt);
 753                                                 }
 754                                                 else
 755                                                 {
 756                                                         ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Mat_Pt"))->Fill(simPt);
 757                                                 }
 758                                         }
 759                                 }
 760                         }
 761                 }
 762
 763                 // get the pid
 764                 Int_t pid = fLnID->GetPID( fPartCode, pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta, fMaxNSigmaITS, fMaxNSigmaTPC, fMaxNSigmaTOF);
 765
 766                 // pid table for prior probabilities (only Bayes)
 767
 768                 Int_t offset = AliPID::kDeuteron - 5;
 769
 770                 if( fSimulation )
 771                 {
 772                         Int_t sim = simpid > AliPID::kProton ? simpid - offset : simpid;
 773                         Int_t rec = pid > AliPID::kProton ? pid - offset : pid;
 774
 775                         if((sim > -1) && (rec > -1)) // pid performance
 776                         {
 777                                 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(sim, rec);
 778                                 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(9, rec);
 779                         }
 780                         if(sim > -1)
 781                         {
 782                                 ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(sim, 9);
 783                         }
 784                 }
 785
 786                 // for bayes iteration
 787                 if(pid != -1)
 788                 {
 789                         Int_t iBin = pid > AliPID::kProton ? pid - offset : pid;
 790                         ((TH1D*)fHistoMap->Get(fSpecies + "_Stats_PID"))->Fill(iBin);
 791                 }
 792
 793                 if(pid != fPartCode) continue;
 794
 795                 Bool_t goodPid = 0;
 796
 797                 if(fSimulation)
 798                 {
 799                         goodPid = ( simpid == pid );
 800                 }
 801
 802                 ((TH1D*)fHistoMap->Get(particle + "_PID_Ntrk_pTPC"))->Fill(pTPC,fNtrk);
 803                 ((TH1D*)fHistoMap->Get(particle + "_PID_Zmult_pTPC"))->Fill(pTPC,fKNOmult);
 804
 805                 // pid performance
 806                 ((TH2D*)fHistoMap->Get(particle + "_PID_ITSdEdx_P"))->Fill(pITS, dEdxITS);
 807                 ((TH2D*)fHistoMap->Get(particle + "_PID_TPCdEdx_P"))->Fill(pTPC, dEdxTPC);
 808
 809                 if(this->TOFmatch(iTrack))
 810                 {
 811                         ((TH2D*)fHistoMap->Get(particle + "_PID_Beta_P"))->Fill(pTOF, beta);
 812                         ((TH2D*)fHistoMap->Get(particle + "_PID_Mass_P"))->Fill(pTOF, mass);
 813                         if(fSimulation && goodPid)
 814                         {
 815                                 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Mass"))->Fill(mass);
 816                         }
 817                 }
 818
 819                 // ---------------------------------
 820                 //  results in |y| < fMaxY
 821                 // ---------------------------------
 822
 823                 ((TH1D*)fHistoMap->Get(particle + "_PID_Y"))->Fill(y);
 824                 ((TH2D*)fHistoMap->Get(particle + "_PID_Pt_Y"))->Fill(y, pt);
 825
 826                 if(TMath::Abs(y) >= fMaxY) continue;
 827
 828                 ((TH1D*)fHistoMap->Get(particle + "_PID_Pt"))->Fill(pt);
 829                 ((TH1D*)fHistoMap->Get(particle + "_PID_Phi"))->Fill(phi);
 830
 831                 if(iTrack->IsOn(AliESDtrack::kTRDin))
 832                 {
 833                         ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_Pt"))->Fill(pt);
 834
 835                         if(iTrack->IsOn(AliESDtrack::kTRDout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_TRDout_Pt"))->Fill(pt);
 836                         if(iTrack->IsOn(AliESDtrack::kTOFout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TRDin_TOFout_Pt"))->Fill(pt);
 837                 }
 838
 839                 if(iTrack->IsOn(AliESDtrack::kTOFin))
 840                 {
 841                         ((TH1D*)fHistoMap->Get(particle + "_PID_TOFin_Pt"))->Fill(pt);
 842
 843                         if(iTrack->IsOn(AliESDtrack::kTOFout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TOFin_TOFout_Pt"))->Fill(pt);
 844                 }
 845
 846                 if( this->TOFmatch(iTrack) )
 847                 {
 848                         ((TH2D*)fHistoMap->Get(particle + "_PID_M2_Pt"))->Fill(pt, m2);
 849                         ((TH1D*)fHistoMap->Get(particle + "_PID_TOFmatch_Pt"))->Fill(pt);
 850                 }
 851
 852                 // secondaries
 853
 854                 Bool_t m2match = kTRUE;
 855
 856                 if( fTOFmatch && (fLnID->GetPidProcedure() > AliLnID::kMaxLikelihood))
 857                 {
 858                         if((m2 < fMinM2) || (m2 >= fMaxM2)) m2match = kFALSE;
 859                 }
 860
 861                 if(m2match)
 862                 {
 863                         ((TH2D*)fHistoMap->Get(particle + "_PID_DCAxy_Pt"))->Fill(pt, dcaxy);
 864                         ((TH2D*)fHistoMap->Get(particle + "_PID_DCAz_Pt"))->Fill(pt, dcaz);
 865                         ((TH2D*)fHistoMap->Get(particle + "_PID_NSigma_Pt"))->Fill(pt, nSigmaVtx);
 866                 }
 867
 868                 if(fSimulation && goodPid)
 869                 {
 870                         // for unfolding and pid contamination
 871                         if(!this->IsFakeTrack(iTrack))
 872                         {
 873                                 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Pt"))->Fill(simPt);
 874
 875                                 if( this->TOFmatch(iTrack) )
 876                                 {
 877                                         ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_M2_Pt"))->Fill(pt, m2);
 878                                 }
 879
 880                                 if(this->IsPhysicalPrimary(iParticle))
 881                                 {
 882                                         ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Prim_Pt"))->Fill(simPt);
 883                                 }
 884
 885                                 if(m2match)
 886                                 {
 887                                         if(this->IsPhysicalPrimary(iParticle))
 888                                         {
 889                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_DCAxy_Pt"))->Fill(pt, dcaxy);
 890                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_DCAz_Pt"))->Fill(pt, dcaz);
 891                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_NSigma_Pt"))->Fill(pt, nSigmaVtx);
 892                                         }
 893                                         else if(this->IsFromWeakDecay(iParticle))
 894                                         {
 895                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_DCAxy_Pt"))->Fill(pt, dcaxy);
 896                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_DCAz_Pt"))->Fill(pt, dcaz);
 897                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_NSigma_Pt"))->Fill(pt, nSigmaVtx);
 898                                         }
 899                                         else // from materials
 900                                         {
 901                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_DCAxy_Pt"))->Fill(pt, dcaxy);
 902                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_DCAz_Pt"))->Fill(pt, dcaz);
 903                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_NSigma_Pt"))->Fill(pt, nSigmaVtx);
 904                                         }
 905                                 }
 906                         }
 907                         else // fake tracks
 908                         {
 909                                 ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Pt"))->Fill(simPt);
 910
 911                                 if(m2match)
 912                                 {
 913                                         if(this->IsPhysicalPrimary(iParticle))
 914                                         {
 915                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Prim_DCAxy_Pt"))->Fill(pt, dcaxy);
 916                                         }
 917                                         else if(this->IsFromWeakDecay(iParticle))
 918                                         {
 919                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt"))->Fill(pt, dcaxy);
 920                                         }
 921                                         else // from materials
 922                                         {
 923                                                 ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Mat_DCAxy_Pt"))->Fill(pt, dcaxy);
 924                                         }
 925                                 }
 926                         }
 927                 }
 928         }
 929
 930         return nTracks;
 931 }
 932
 933 void AliAnalysisTaskB2::Terminate(Option_t* )
 934 {
 935   // The Terminate() function is the last function to be called during
 936   // a query. It always runs on the client, it can be used to present
 937   // the results graphically or save the results to file.
 938
 939 }
 940
 941 Bool_t AliAnalysisTaskB2::IsV0AND() const
 942 {
 943 //
 944 // signals in both V0A and V0C
 945 //
 946         return ( fTrigAna->IsOfflineTriggerFired(fESDevent, AliTriggerAnalysis::kV0A) &&
 947                          fTrigAna->IsOfflineTriggerFired(fESDevent, AliTriggerAnalysis::kV0C) );
 948 }
 949
 950 Bool_t AliAnalysisTaskB2::IsFastOnly(UInt_t triggerBits) const
 951 {
 952 //
 953 // kFastOnly trigger
 954 //
 955         return ( (triggerBits&AliVEvent::kFastOnly) == AliVEvent::kFastOnly );
 956 }
 957
 958 Bool_t AliAnalysisTaskB2::IsMB(UInt_t triggerBits) const
 959 {
 960 //
 961 // MB event
 962 //
 963         return ( (triggerBits&AliVEvent::kMB) == AliVEvent::kMB );
 964 }
 965
 966 Bool_t AliAnalysisTaskB2::TOFmatch(const AliESDtrack* trk) const
 967 {
 968 //
 969 // Check for TOF match signal
 970 //
 971         return ( trk->IsOn(AliESDtrack::kTOFout) && trk->IsOn(AliESDtrack::kTIME) );
 972 }
 973
 974 Bool_t AliAnalysisTaskB2::AcceptTOFtrack(const AliESDtrack* trk) const
 975 {
 976 //
 977 // Additional checks for TOF match signal
 978 //
 979         if( !this->TOFmatch(trk) ) return kFALSE;
 980         if( trk->GetIntegratedLength() < 350) return kFALSE;
 981         if( trk->GetTOFsignal() < 1e-6) return kFALSE;
 982
 983         return kTRUE;
 984 }
 985
 986 TParticle* AliAnalysisTaskB2::GetParticle(const AliESDtrack* trk) const
 987 {
 988 //
 989 // Particle that left the track
 990 //
 991         AliStack* stack = fMCevent->Stack();
 992
 993         Int_t label = TMath::Abs(trk->GetLabel()); // if negative then it shares points from other tracks
 994         if( label >= fMCevent->GetNumberOfTracks() ) return 0;
 995
 996         return stack->Particle(label);
 997 }
 998
 999 Bool_t AliAnalysisTaskB2::IsFakeTrack(const AliESDtrack* trk) const
1000 {
1001 //
1002 // Check if the track shares some clusters with different particles
1003 // (definition changed to label=0? )
1004 //
1005         return ( trk->GetLabel() < 0 );
1006 }
1007
1008 Bool_t AliAnalysisTaskB2::IsPhysicalPrimary(const TParticle* prt) const
1009 {
1010 //
1011 // Check if the particle is physical primary
1012 //
1013         AliStack* stack = fMCevent->Stack();
1014         Int_t index = stack->Particles()->IndexOf(prt);
1015
1016         return stack->IsPhysicalPrimary(index);
1017 }
1018
1019 Bool_t AliAnalysisTaskB2::IsFromMaterial(const TParticle* prt) const
1020 {
1021 //
1022 // Check if the particle is originated at the materials
1023 //
1024         AliStack* stack = fMCevent->Stack();
1025         Int_t index = stack->Particles()->IndexOf(prt);
1026
1027         return stack->IsSecondaryFromMaterial(index);
1028 }
1029
1030 Bool_t AliAnalysisTaskB2::IsFromWeakDecay(const TParticle* prt) const
1031 {
1032 //
1033 // Check if the particle comes from a weak decay
1034 //
1035         AliStack* stack = fMCevent->Stack();
1036         Int_t index = stack->Particles()->IndexOf(prt);
1037
1038         return stack->IsSecondaryFromWeakDecay(index);
1039 }
1040
1041 Double_t AliAnalysisTaskB2::GetSign(TParticle* prt) const
1042 {
1043 //
1044 // Sign of the particle
1045 //
1046         TParticlePDG* pdg = prt->GetPDG();
1047
1048         if(pdg != 0) return pdg->Charge();
1049
1050         return 0;
1051 }
1052
1053 Double_t AliAnalysisTaskB2::GetPhi(const AliESDtrack* trk) const
1054 {
1055 //
1056 // Azimuthal angle [0,2pi) using the pt at the DCA point
1057 //
1058         Double_t px = trk->Px();
1059         Double_t py = trk->Py();
1060
1061         return TMath::Pi()+TMath::ATan2(-py, -px);
1062 }
1063
1064 Double_t AliAnalysisTaskB2::GetTheta(const AliESDtrack* trk) const
1065 {
1066 //
1067 // Polar angle using the pt at the DCA point
1068 //
1069         Double_t p  = trk->GetP();
1070         Double_t pz = trk->Pz();
1071
1072         return (pz == 0) ? TMath::PiOver2() : TMath::ACos(pz/p);
1073 }
1074
1075 Double_t AliAnalysisTaskB2::GetRapidity(const AliESDtrack* trk, Int_t pid) const
1076 {
1077 //
1078 // Rapidity assuming the given particle hypothesis
1079 // and using the momentum at the DCA
1080 //
1081         Double_t m  = AliPID::ParticleMass(pid);
1082         Double_t p  = (pid>AliPID::kTriton) ? 2.*trk->GetP() : trk->GetP();
1083         Double_t e  = TMath::Sqrt(p*p + m*m);
1084         Double_t pz = (pid>AliPID::kTriton) ? 2.*trk->Pz() : trk->Pz();
1085
1086         if(e <= pz) return 1.e+16;
1087
1088         return 0.5*TMath::Log( (e+pz)/(e-pz) );
1089 }
1090
1091 Double_t AliAnalysisTaskB2::GetITSmomentum(const AliESDtrack* trk) const
1092 {
1093 //
1094 // Momentum for ITS pid
1095 //
1096         Double_t pDCA = trk->GetP();
1097         Double_t pTPC = trk->GetTPCmomentum();
1098
1099         return (pDCA+pTPC)/2.;
1100 }
1101
1102 Double_t AliAnalysisTaskB2::GetTOFmomentum(const AliESDtrack* trk) const
1103 {
1104 //
1105 // Momentum for TOF pid
1106 //
1107         Double_t pDCA = trk->GetP();
1108
1109         const AliExternalTrackParam* param = trk->GetOuterParam();
1110         Double_t pOut = param ? param->GetP() : trk->GetP();
1111
1112         return (pDCA+pOut)/2.;
1113 }
1114
1115 Double_t AliAnalysisTaskB2::GetBeta(const AliESDtrack* trk) const
1116 {
1117 //
1118 // Velocity
1119 //
1120         Double_t t = this->GetTimeOfFlight(trk); // ps
1121         Double_t l = trk->GetIntegratedLength(); // cm
1122
1123         if(t <= 0) return 1.e6; // 1M times the speed of light ;)
1124
1125         return (l/t)/2.99792458e-2;
1126 }
1127
1128 Double_t AliAnalysisTaskB2::GetMassSquare(const AliESDtrack* trk) const
1129 {
1130 //
1131 // Square mass
1132 //
1133         Double_t p = (fPartCode>AliPID::kTriton) ? 2.*this->GetTOFmomentum(trk) : this->GetTOFmomentum(trk);
1134         Double_t beta = this->GetBeta(trk);
1135
1136         return p*p*(1./(beta*beta) - 1.);
1137 }
1138
1139 Int_t AliAnalysisTaskB2::GetChargedMultiplicity(Double_t etaMax) const
1140 {
1141 //
1142 // Charged particle multiplicity using ALICE physical primary definition
1143 //
1144         AliStack* stack = fMCevent->Stack();
1145
1146         Int_t nch = 0;
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
1166 Double_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
1190 Int_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
1206 Double_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
1219 Int_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
1235 Int_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 }