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