PWGLF/SPECTRA/IdentifiedHighPt/grid/AliAnalysisTaskHighPtDeDx.cxx

   1 #include "AliAnalysisTaskHighPtDeDx.h"
   2
   3 // ROOT includes
   4 #include <TList.h>
   5 #include <TTree.h>
   6 #include <TMath.h>
   7 #include <TH1.h>
   8 #include <TParticle.h>
   9 #include <TFile.h>
  10
  11 // AliRoot includes
  12 #include <AliAnalysisManager.h>
  13 #include <AliAnalysisFilter.h>
  14 #include <AliESDInputHandler.h>
  15 #include <AliESDEvent.h>
  16 #include <AliESDVertex.h>
  17 #include <AliLog.h>
  18 #include <AliExternalTrackParam.h>
  19 #include <AliESDtrackCuts.h>
  20 #include <AliESDVZERO.h>
  21 #include <AliAODVZERO.h>
  22
  23 #include <AliMCEventHandler.h>
  24 #include <AliMCEvent.h>
  25 #include <AliStack.h>
  26
  27 #include <TTreeStream.h>
  28
  29 #include <AliHeader.h>
  30 #include <AliGenPythiaEventHeader.h>
  31 #include <AliGenDPMjetEventHeader.h>
  32
  33 #include "AliCentrality.h"
  34
  35 #include <AliAODTrack.h>
  36 #include <AliAODPid.h>
  37 #include <AliAODMCHeader.h>
  38
  39
  40 // STL includes
  41 #include <iostream>
  42 using namespace std;
  43
  44
  45 //
  46 // Responsible:
  47 // Alexandru Dobrin (Wayne State)
  48 // Peter Christiansen (Lund)
  49 //
  50
  51 /*
  52 To do:
  53
  54 Separate the code into two
  55
  56 */
  57
  58
  59
  60
  61 ClassImp(AliAnalysisTaskHighPtDeDx)
  62
  63 const Double_t AliAnalysisTaskHighPtDeDx::fgkClight = 2.99792458e-2;
  64
  65 //_____________________________________________________________________________
  66 AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx():
  67   AliAnalysisTaskSE(),
  68   fESD(0x0),
  69   fAOD(0x0),
  70   fMC(0x0),
  71   fMCStack(0x0),
  72   fMCArray(0x0),
  73   fTrackFilter(0x0),
  74   fTrackFilterGolden(0x0),
  75   fTrackFilterTPC(0x0),
  76   fAnalysisType("ESD"),
  77   fAnalysisMC(kFALSE),
  78   fAnalysisPbPb(kFALSE),
  79   fTPCBranch(kFALSE),
  80   ftrigBit1(0x0),
  81   ftrigBit2(0x0),
  82   fRandom(0x0),
  83   fEvent(0x0),
  84   fTrackArrayGlobalPar(0x0),
  85   fTrackArrayTPCPar(0x0),
  86   fTrackArrayMC(0x0),
  87   fVZEROArray(0x0),
  88   fVtxCut(10.0),
  89   fEtaCut(0.9),
  90   fMinPt(0.1),
  91   fMinCent(0.0),
  92   fMaxCent(100.0),
  93   fLowPtFraction(0.01),
  94   fTreeOption(0),
  95   fMcProcessType(-999),
  96   fTriggeredEventMB(-999),
  97   fVtxStatus(-999),
  98   fZvtx(-999),
  99   fZvtxMC(-999),
 100   fRun(-999),
 101   fEventId(-999),
 102   fListOfObjects(0),
 103   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
 104   fTree(0x0),
 105   fn1(0),
 106   fn2(0),
 107   fcent(0)
 108
 109
 110 {
 111   // Default constructor (should not be used)
 112
 113   //  fRandom = new TRandom(0); // 0 means random seed
 114 }
 115
 116 //______________________________________________________________________________
 117 AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx(const char *name):
 118   AliAnalysisTaskSE(name),
 119   fESD(0x0),
 120   fAOD(0x0),
 121   fMC(0x0),
 122   fMCStack(0x0),
 123   fMCArray(0x0),
 124   fTrackFilter(0x0),
 125   fTrackFilterGolden(0x0),
 126   fTrackFilterTPC(0x0),
 127   fAnalysisType("ESD"),
 128   fAnalysisMC(kFALSE),
 129   fAnalysisPbPb(kFALSE),
 130   fTPCBranch(kFALSE),
 131   ftrigBit1(0x0),
 132   ftrigBit2(0x0),
 133   fRandom(0x0),
 134   fEvent(0x0),
 135   fTrackArrayGlobalPar(0x0),
 136   fTrackArrayMC(0x0),
 137   fVtxCut(10.0),
 138   fEtaCut(0.9),
 139   fMinPt(0.1),
 140   fLowPtFraction(0.01),
 141   fTreeOption(0),
 142   fMcProcessType(-999),
 143   fTriggeredEventMB(-999),
 144   fVtxStatus(-999),
 145   fZvtx(-999),
 146   fZvtxMC(-999),
 147   fRun(-999),
 148   fEventId(-999),
 149   fListOfObjects(0),
 150   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
 151   fTree(0x0),
 152   fn1(0),
 153   fn2(0),
 154   fcent(0)
 155
 156 {
 157
 158   if(fTree)fTree=0;
 159   // Output slot #1 writes into a TList
 160   DefineOutput(1, TList::Class());
 161 }
 162
 163 //_____________________________________________________________________________
 164 AliAnalysisTaskHighPtDeDx::~AliAnalysisTaskHighPtDeDx()
 165 {
 166   // Destructor
 167   // histograms are in the output list and deleted when the output
 168   // list is deleted by the TSelector dtor
 169   if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 170     delete fListOfObjects;
 171     fListOfObjects = 0;
 172   }
 173   if (fRandom) delete fRandom;
 174   fRandom=0;
 175
 176   // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse
 177   // if (fListOfObjects  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
 178
 179
 180
 181 }
 182
 183 //______________________________________________________________________________
 184 void AliAnalysisTaskHighPtDeDx::UserCreateOutputObjects()
 185 {
 186   // This method is called once per worker node
 187   // Here we define the output: histograms and debug tree if requested
 188   // We also create the random generator here so it might get different seeds...
 189   fRandom = new TRandom(0); // 0 means random seed
 190
 191   OpenFile(1);
 192   fListOfObjects = new TList();
 193   fListOfObjects->SetOwner();
 194
 195   //
 196   // Histograms
 197   //
 198   fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
 199   fListOfObjects->Add(fEvents);
 200
 201   fn1=new TH1F("fn1","fn1",5001,-1,5000);
 202   fListOfObjects->Add(fn1);
 203
 204   fn2=new TH1F("fn2","fn2",5001,-1,5000);
 205   fListOfObjects->Add(fn2);
 206
 207   fcent=new TH1F("fcent","fcent",104,-2,102);
 208   fListOfObjects->Add(fcent);
 209
 210   fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
 211   fListOfObjects->Add(fVtx);
 212
 213   fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
 214   fListOfObjects->Add(fVtxBeforeCuts);
 215
 216   fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
 217   fListOfObjects->Add(fVtxAfterCuts);
 218
 219   if (fAnalysisMC) {
 220     fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
 221     fListOfObjects->Add(fVtxMC);
 222   }
 223
 224   if (fTreeOption) {
 225
 226     fTree = new TTree("tree","Event data");
 227     fEvent = new DeDxEvent();
 228     fTree->Branch("event", &fEvent);
 229
 230     fTrackArrayGlobalPar = new TClonesArray("DeDxTrack", 1000);
 231     fTree->Bronch("trackGlobalPar", "TClonesArray", &fTrackArrayGlobalPar);
 232     if(fTPCBranch){
 233       fTrackArrayTPCPar = new TClonesArray("DeDxTrack", 1000);
 234       fTree->Bronch("trackTPCPar", "TClonesArray", &fTrackArrayTPCPar);
 235     }
 236
 237     fVZEROArray = new TClonesArray("VZEROCell", 1000);
 238     fTree->Bronch("cellVZERO", "TClonesArray", &fVZEROArray);
 239
 240     if (fAnalysisMC) {
 241       fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
 242       fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
 243     }
 244
 245     fTree->SetDirectory(0);
 246
 247     fListOfObjects->Add(fTree);
 248
 249   }
 250
 251   // Post output data.
 252   PostData(1, fListOfObjects);
 253 }
 254
 255 //______________________________________________________________________________
 256 void AliAnalysisTaskHighPtDeDx::UserExec(Option_t *)
 257 {
 258   // Main loop
 259
 260   //
 261   // First we make sure that we have valid input(s)!
 262   //
 263   AliVEvent *event = InputEvent();
 264   if (!event) {
 265      Error("UserExec", "Could not retrieve event");
 266      return;
 267   }
 268
 269
 270   if (fAnalysisType == "ESD"){
 271     fESD = dynamic_cast<AliESDEvent*>(event);
 272     if(!fESD){
 273       Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 274       this->Dump();
 275       return;
 276     }
 277   } else {
 278     fAOD = dynamic_cast<AliAODEvent*>(event);
 279     if(!fAOD){
 280       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 281       this->Dump();
 282       return;
 283     }
 284   }
 285
 286
 287
 288   if (fAnalysisMC) {
 289
 290     if (fAnalysisType == "ESD"){
 291       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 292       if(!fMC){
 293         Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 294         this->Dump();
 295         return;
 296       }
 297
 298       fMCStack = fMC->Stack();
 299
 300       if(!fMCStack){
 301         Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
 302         this->Dump();
 303         return;
 304       }
 305     } else { // AOD
 306
 307       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 308       if(fMC)
 309         fMC->Dump();
 310
 311       fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
 312       if(!fMCArray){
 313         Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
 314         this->Dump();
 315         return;
 316       }
 317     }
 318   }
 319
 320
 321   // Get trigger decision
 322   fTriggeredEventMB = 0; //init
 323
 324
 325   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 326      ->IsEventSelected() & ftrigBit1 ){
 327     fn1->Fill(1);
 328     fTriggeredEventMB = 1;  //event triggered as minimum bias
 329   }
 330   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 331      ->IsEventSelected() & ftrigBit2 ){
 332     // From AliVEvent:
 333     //    kINT7         = BIT(1), // V0AND trigger, offline V0 selection
 334     fTriggeredEventMB += 2;
 335     fn2->Fill(1);
 336   }
 337
 338
 339   // Get process type for MC
 340   fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
 341
 342   // real data that are not triggered we skip
 343   if(!fAnalysisMC && !fTriggeredEventMB)
 344     return;
 345
 346   if (fAnalysisMC) {
 347
 348     if (fAnalysisType == "ESD"){
 349
 350       AliHeader* headerMC = fMC->Header();
 351       if (headerMC) {
 352
 353         AliGenEventHeader* genHeader = headerMC->GenEventHeader();
 354         TArrayF vtxMC(3); // primary vertex  MC
 355         vtxMC[0]=9999; vtxMC[1]=9999;  vtxMC[2]=9999; //initialize with dummy
 356         if (genHeader) {
 357           genHeader->PrimaryVertex(vtxMC);
 358         }
 359         fZvtxMC = vtxMC[2];
 360
 361         // PYTHIA:
 362         AliGenPythiaEventHeader* pythiaGenHeader =
 363           dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
 364         if (pythiaGenHeader) {  //works only for pythia
 365           fMcProcessType =  GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
 366         }
 367         // PHOJET:
 368         AliGenDPMjetEventHeader* dpmJetGenHeader =
 369           dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
 370         if (dpmJetGenHeader) {
 371           fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
 372         }
 373       }
 374     } else { // AOD
 375
 376       AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
 377       if(mcHeader) {
 378         fZvtxMC = mcHeader->GetVtxZ();
 379
 380         if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
 381           fMcProcessType =  GetPythiaEventProcessType(mcHeader->GetEventType());
 382         } else {
 383           fMcProcessType =  GetDPMjetEventProcessType(mcHeader->GetEventType());
 384         }
 385       }
 386     }
 387   }
 388
 389   // There are 3 cases
 390   // Vertex: NO  - status -1
 391   // Vertex: YES : outside cut - status 0
 392   //             : inside cut  - status 1
 393   // We have to be careful how we normalize because we probably want to
 394   // normalize to:
 395   // Nevents=(No vertex + outside + inside)/(out + in)*in
 396
 397
 398   if (fAnalysisType == "ESD")
 399     fZvtx = GetVertex(fESD);
 400   else // AOD
 401     fZvtx = GetVertex(fAOD);
 402
 403   fVtxStatus = -999;
 404
 405   if(fZvtx<-990) {
 406
 407     fVtxStatus = -1;
 408     if(fTriggeredEventMB)
 409       fVtx->Fill(0);
 410     if(fAnalysisMC)
 411       fVtxMC->Fill(0);
 412   } else {
 413
 414     if(fTriggeredEventMB)
 415       fVtx->Fill(1);
 416     if(fAnalysisMC)
 417       fVtxMC->Fill(1);
 418     fVtxBeforeCuts->Fill(fZvtx);
 419     fVtxStatus = 0;
 420     if (TMath::Abs(fZvtx) < fVtxCut) {
 421       fVtxAfterCuts->Fill(fZvtx);
 422       fVtxStatus = 1;
 423     }
 424   }
 425
 426   // store MC event data no matter what
 427   if(fAnalysisMC) {
 428
 429     if (fAnalysisType == "ESD") {
 430       ProcessMCTruthESD();
 431     } else { // AOD
 432       ProcessMCTruthAOD();
 433     }
 434   }
 435
 436
 437
 438   Float_t centrality = -10;
 439   // only analyze triggered events
 440   if(fTriggeredEventMB) {
 441
 442     if (fAnalysisType == "ESD"){
 443       if(fAnalysisPbPb){
 444         AliCentrality *centObject = fESD->GetCentrality();
 445         centrality = centObject->GetCentralityPercentile("V0M");
 446         if((centrality>fMaxCent)||(centrality<fMinCent))return;
 447       }
 448       fcent->Fill(centrality);
 449       AnalyzeESD(fESD);
 450     } else { // AOD
 451       if(fAnalysisPbPb){
 452         AliCentrality *centObject = fAOD->GetCentrality();
 453         if(centObject)
 454           centrality = centObject->GetCentralityPercentile("V0M");
 455         if((centrality>fMaxCent)||(centrality<fMinCent))return;
 456       }
 457       fcent->Fill(centrality);
 458       AnalyzeAOD(fAOD);
 459     }
 460   }
 461
 462   if( fTreeOption) {
 463
 464     fEvent->process = fMcProcessType;
 465     fEvent->trig    = fTriggeredEventMB;
 466     fEvent->zvtxMC  = fZvtxMC;
 467     fEvent->cent      = centrality;
 468
 469     fTree->Fill();
 470     fTrackArrayGlobalPar->Clear();
 471     if(fTPCBranch)
 472       fTrackArrayTPCPar->Clear();
 473     fVZEROArray->Clear();
 474
 475     if (fAnalysisMC)
 476       fTrackArrayMC->Clear();
 477   }
 478
 479   // Post output data.
 480   PostData(1, fListOfObjects);
 481 }
 482
 483 //________________________________________________________________________
 484 void AliAnalysisTaskHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
 485 {
 486   fRun  = esdEvent->GetRunNumber();
 487   fEventId = 0;
 488   if(esdEvent->GetHeader())
 489     fEventId = GetEventIdAsLong(esdEvent->GetHeader());
 490
 491   Short_t isPileup = esdEvent->IsPileupFromSPD();
 492
 493   //  Int_t     event     = esdEvent->GetEventNumberInFile();
 494   UInt_t    time      = esdEvent->GetTimeStamp();
 495   //  ULong64_t trigger   = esdEvent->GetTriggerMask();
 496   Float_t   magf      = esdEvent->GetMagneticField();
 497
 498
 499
 500
 501
 502   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
 503
 504     // accepted event
 505     fEvents->Fill(0);
 506
 507
 508     if(fVtxStatus!=1) return; // accepted vertex
 509     Int_t nESDTracks = esdEvent->GetNumberOfTracks();
 510
 511     if(nESDTracks<1)return;
 512     cout<<"Multiplicity="<<nESDTracks<<endl;
 513
 514     ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
 515     if(fTPCBranch)
 516       ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
 517
 518     fEvents->Fill(1);
 519     AliESDVZERO *esdV0 = esdEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
 520
 521     for (Int_t iCh=0; iCh<64; ++iCh) {
 522       Float_t multv=esdV0->GetMultiplicity(iCh);
 523       Int_t intexv=iCh;
 524       VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
 525       cellv0->cellindex=intexv;
 526       cellv0->cellmult= multv;
 527     }
 528
 529
 530
 531   } // end if triggered
 532
 533   if(fTreeOption) {
 534
 535     fEvent->run       = fRun;
 536     fEvent->eventid   = fEventId;
 537     fEvent->time      = time;
 538     //fEvent->cent      = centrality;
 539     fEvent->mag       = magf;
 540     fEvent->zvtx      = fZvtx;
 541     fEvent->vtxstatus = fVtxStatus;
 542     fEvent->pileup    = isPileup;
 543
 544   }
 545
 546
 547
 548
 549 }
 550
 551 //________________________________________________________________________
 552 void AliAnalysisTaskHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
 553 {
 554   fRun  = aodEvent->GetRunNumber();
 555   fEventId = 0;
 556   if(aodEvent->GetHeader())
 557     fEventId = GetEventIdAsLong(aodEvent->GetHeader());
 558
 559   UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
 560   Float_t   magf      = aodEvent->GetMagneticField();
 561
 562   //Int_t     trackmult = 0; // no pt cuts
 563   //Int_t     nadded    = 0;
 564
 565   Short_t isPileup = aodEvent->IsPileupFromSPD();
 566
 567
 568
 569
 570   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
 571
 572     // accepted event
 573     fEvents->Fill(0);
 574
 575     if(fVtxStatus!=1) return; // accepted vertex
 576     Int_t nAODTracks = aodEvent->GetNumberOfTracks();
 577     if(nAODTracks<1)return;
 578
 579     ProduceArrayTrksAOD( aodEvent, kGlobalTrk );
 580     if(fTPCBranch)
 581       ProduceArrayTrksAOD( aodEvent, kTPCTrk );
 582     fEvents->Fill(1);
 583
 584     AliAODVZERO *esdV0 = aodEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
 585
 586     for (Int_t iCh=0; iCh<64; ++iCh) {
 587       Float_t multv=esdV0->GetMultiplicity(iCh);
 588       Int_t intexv=iCh;
 589       VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
 590       cellv0->cellindex=intexv;
 591       cellv0->cellmult= multv;
 592     }
 593
 594
 595
 596
 597   } // end if triggered
 598
 599   if(fTreeOption) {
 600
 601     //Sort(fTrackArrayGlobalPar, kFALSE);
 602
 603     fEvent->run       = fRun;
 604     fEvent->eventid   = fEventId;
 605     fEvent->time      = time;
 606     //fEvent->cent      = centrality;
 607     fEvent->mag       = magf;
 608     fEvent->zvtx      = fZvtx;
 609     fEvent->vtxstatus = fVtxStatus;
 610     //fEvent->trackmult = trackmult;
 611     //fEvent->n         = nadded;
 612     fEvent->pileup    = isPileup;
 613   }
 614 }
 615
 616 //_____________________________________________________________________________
 617 Float_t AliAnalysisTaskHighPtDeDx::GetVertex(const AliVEvent* event) const
 618 {
 619   Float_t zvtx = -999;
 620
 621   const AliVVertex* primaryVertex = event->GetPrimaryVertex();
 622
 623   if(primaryVertex->GetNContributors()>0)
 624     zvtx = primaryVertex->GetZ();
 625
 626   return zvtx;
 627 }
 628
 629 //_____________________________________________________________________________
 630 Short_t AliAnalysisTaskHighPtDeDx::GetPidCode(Int_t pdgCode) const
 631 {
 632   // return our internal code for pions, kaons, and protons
 633
 634   Short_t pidCode = 6;
 635
 636   switch (TMath::Abs(pdgCode)) {
 637   case 211:
 638     pidCode = 1; // pion
 639     break;
 640   case 321:
 641     pidCode = 2; // kaon
 642     break;
 643   case 2212:
 644     pidCode = 3; // proton
 645     break;
 646   case 11:
 647     pidCode = 4; // electron
 648     break;
 649   case 13:
 650     pidCode = 5; // muon
 651     break;
 652   default:
 653     pidCode = 6;  // something else?
 654   };
 655
 656   return pidCode;
 657 }
 658
 659
 660 //_____________________________________________________________________________
 661 void AliAnalysisTaskHighPtDeDx::ProcessMCTruthESD()
 662 {
 663   // Fill the special MC histogram with the MC truth info
 664
 665   Short_t trackmult = 0;
 666   Short_t nadded    = 0;
 667   const Int_t nTracksMC = fMCStack->GetNtrack();
 668
 669   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
 670
 671     //Cuts
 672     if(!(fMCStack->IsPhysicalPrimary(iTracks)))
 673       continue;
 674
 675     TParticle* trackMC = fMCStack->Particle(iTracks);
 676
 677     Double_t chargeMC = trackMC->GetPDG()->Charge();
 678     if (chargeMC == 0)
 679       continue;
 680
 681     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
 682       continue;
 683
 684     trackmult++;
 685
 686     //   "charge:pt:p:eta:phi:pidCode"
 687     Float_t ptMC      = trackMC->Pt();
 688     Float_t pMC       = trackMC->P();
 689     Float_t etaMC     = trackMC->Eta();
 690     Float_t phiMC     = trackMC->Phi();
 691
 692     Int_t pdgCode = trackMC->GetPdgCode();
 693     Short_t pidCodeMC = 0;
 694     pidCodeMC = GetPidCode(pdgCode);
 695
 696     // Here we want to add some of the MC histograms!
 697
 698     // And therefore we first cut here!
 699     if (trackMC->Pt() < fMinPt) {
 700
 701       // Keep small fraction of low pT tracks
 702       if(fRandom->Rndm() > fLowPtFraction)
 703         continue;
 704     } // else {
 705     // Here we want to add the high pt part of the MC histograms!
 706     //    }
 707
 708     if(fTreeOption) {
 709
 710       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
 711       nadded++;
 712
 713       track->pMC   = pMC;
 714       track->ptMC  = ptMC;
 715       track->etaMC = etaMC;
 716       track->phiMC = phiMC;
 717       track->qMC   = Short_t(chargeMC);
 718       track->pidMC = pidCodeMC;
 719       track->pdgMC = pdgCode;
 720     }
 721
 722   }//MC track loop
 723
 724   if(fTreeOption) {
 725
 726     Sort(fTrackArrayMC, kTRUE);
 727
 728     fEvent->trackmultMC = trackmult;
 729     fEvent->nMC         = nadded;
 730   }
 731
 732 }
 733
 734 //_____________________________________________________________________________
 735 void AliAnalysisTaskHighPtDeDx::ProcessMCTruthAOD()
 736 {
 737   // Fill the special MC histogram with the MC truth info
 738
 739   Short_t trackmult = 0;
 740   Short_t nadded    = 0;
 741   const Int_t nTracksMC = fMCArray->GetEntriesFast();
 742
 743   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
 744
 745     AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
 746
 747     //Cuts
 748     if(!(trackMC->IsPhysicalPrimary()))
 749       continue;
 750
 751
 752     Double_t chargeMC = trackMC->Charge();
 753     if (chargeMC == 0)
 754       continue;
 755
 756     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
 757       continue;
 758
 759     trackmult++;
 760
 761     //   "charge:pt:p:eta:phi:pidCode"
 762     Float_t ptMC      = trackMC->Pt();
 763     Float_t pMC       = trackMC->P();
 764     Float_t etaMC     = trackMC->Eta();
 765     Float_t phiMC     = trackMC->Phi();
 766
 767     Int_t pdgCode = trackMC->PdgCode();
 768     Short_t pidCodeMC = 0;
 769     pidCodeMC = GetPidCode(pdgCode);
 770
 771     // Here we want to add some of the MC histograms!
 772
 773     // And therefore we first cut here!
 774     if (trackMC->Pt() < fMinPt) {
 775
 776       // Keep small fraction of low pT tracks
 777       if(fRandom->Rndm() > fLowPtFraction)
 778         continue;
 779     } // else {
 780     // Here we want to add the high pt part of the MC histograms!
 781     //    }
 782
 783     if(fTreeOption) {
 784
 785       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
 786       nadded++;
 787
 788       track->pMC   = pMC;
 789       track->ptMC  = ptMC;
 790       track->etaMC = etaMC;
 791       track->phiMC = phiMC;
 792       track->qMC   = Short_t(chargeMC);
 793       track->pidMC = pidCodeMC;
 794       track->pdgMC = pdgCode;
 795     }
 796
 797   }//MC track loop
 798
 799   if(fTreeOption) {
 800
 801     Sort(fTrackArrayMC, kTRUE);
 802
 803     fEvent->trackmultMC = trackmult;
 804     fEvent->nMC         = nadded;
 805   }
 806
 807 }
 808
 809 //_____________________________________________________________________________
 810 Short_t AliAnalysisTaskHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
 811   //
 812   // Get the process type of the event.  PYTHIA
 813   //
 814   // source PWG0   dNdpt
 815
 816   Short_t globalType = -1; //init
 817
 818   if(pythiaType==92||pythiaType==93){
 819     globalType = 2; //single diffractive
 820   }
 821   else if(pythiaType==94){
 822     globalType = 3; //double diffractive
 823   }
 824   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
 825   else {
 826     globalType = 1;  //non diffractive
 827   }
 828   return globalType;
 829 }
 830
 831 //_____________________________________________________________________________
 832 Short_t AliAnalysisTaskHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
 833   //
 834   // get the process type of the event.  PHOJET
 835   //
 836   //source PWG0   dNdpt
 837   // can only read pythia headers, either directly or from cocktalil header
 838   Short_t globalType = -1;
 839
 840   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
 841     globalType = 1;
 842   }
 843   else if (dpmJetType==5 || dpmJetType==6) {
 844     globalType = 2;
 845   }
 846   else if (dpmJetType==7) {
 847     globalType = 3;
 848   }
 849   return globalType;
 850 }
 851
 852 //_____________________________________________________________________________
 853 ULong64_t AliAnalysisTaskHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
 854 {
 855   // To have a unique id for each event in a run!
 856   // Modified from AliRawReader.h
 857   return ((ULong64_t)header->GetBunchCrossNumber()+
 858           (ULong64_t)header->GetOrbitNumber()*3564+
 859           (ULong64_t)header->GetPeriodNumber()*16777215*3564);
 860 }
 861
 862
 863 //____________________________________________________________________
 864 TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
 865 {
 866   //
 867   // Finds the first mother among the primary particles of the particle identified by <label>,
 868   // i.e. the primary that "caused" this particle
 869   //
 870   // Taken from AliPWG0Helper class
 871   //
 872
 873   Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
 874   if (motherLabel < 0)
 875     return 0;
 876
 877   return stack->Particle(motherLabel);
 878 }
 879
 880 //____________________________________________________________________
 881 Int_t AliAnalysisTaskHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
 882 {
 883   //
 884   // Finds the first mother among the primary particles of the particle identified by <label>,
 885   // i.e. the primary that "caused" this particle
 886   //
 887   // returns its label
 888   //
 889   // Taken from AliPWG0Helper class
 890   //
 891   const Int_t nPrim  = stack->GetNprimary();
 892
 893   while (label >= nPrim) {
 894
 895     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
 896
 897     TParticle* particle = stack->Particle(label);
 898     if (!particle) {
 899
 900       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
 901       return -1;
 902     }
 903
 904     // find mother
 905     if (particle->GetMother(0) < 0) {
 906
 907       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
 908       return -1;
 909     }
 910
 911     label = particle->GetMother(0);
 912   }
 913
 914   return label;
 915 }
 916
 917 //____________________________________________________________________
 918 AliAODMCParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
 919 {
 920   //
 921   // Finds the first mother among the primary particles of the particle identified by <label>,
 922   // i.e. the primary that "caused" this particle
 923   //
 924   // Taken from AliPWG0Helper class
 925   //
 926
 927   AliAODMCParticle* mcPart = startParticle;
 928
 929   while (mcPart)
 930     {
 931
 932       if(mcPart->IsPrimary())
 933         return mcPart;
 934
 935       Int_t mother = mcPart->GetMother();
 936
 937       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
 938     }
 939
 940   return 0;
 941 }
 942
 943 //_____________________________________________________________________________
 944 void AliAnalysisTaskHighPtDeDx::Sort(TClonesArray* array, Bool_t isMC)
 945 {
 946   const Int_t n = array->GetEntriesFast();
 947   if(n==0) {
 948
 949     if(isMC)
 950       fEvent->ptmaxMC = 0;
 951     else
 952       fEvent->ptmax   = 0;
 953
 954     return;
 955   }
 956
 957   Float_t ptArray[n];
 958   Int_t   indexArray[n];
 959
 960   for(Int_t i = 0; i < n; i++) {
 961
 962     Float_t pt = 0;
 963     if(isMC) {
 964       DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
 965       pt = track->ptMC;
 966     } else {
 967       DeDxTrack* track = (DeDxTrack*)array->At(i);
 968       pt = track->pt;
 969     }
 970     ptArray[i]    = pt;
 971     indexArray[i] = i;
 972   }
 973
 974   TMath::Sort(n, ptArray, indexArray);
 975
 976   // set max pt
 977   if(isMC) {
 978     DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
 979     fEvent->ptmaxMC = track->ptMC;
 980   } else {
 981     DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
 982     fEvent->ptmax   = track->pt;
 983   }
 984
 985   // set order of each track
 986   for(Int_t i = 0; i < n; i++) {
 987
 988     if(isMC) {
 989       DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
 990       track->orderMC = i;
 991     } else {
 992       DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
 993       track->order = i;
 994     }
 995   }
 996 }
 997 //__________________________________________________________________
 998 void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
 999
1000   const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
1001   Int_t trackmult=0;
1002   Int_t nadded=0;
1003   if( analysisMode == kGlobalTrk ){
1004     if(fTrackArrayGlobalPar)
1005       fTrackArrayGlobalPar->Clear();
1006   } else if( analysisMode == kTPCTrk ){
1007     if(fTrackArrayTPCPar)
1008       fTrackArrayTPCPar->Clear();
1009   }
1010
1011   if( analysisMode==kGlobalTrk ){
1012
1013     for(Int_t iT = 0; iT < nESDTracks; iT++) {
1014
1015       AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1016
1017
1018       if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1019         continue;
1020
1021       UShort_t filterFlag = 0;
1022       Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1023       Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1024       Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1025
1026       UInt_t selectDebug = 0;
1027       if (fTrackFilterGolden) {
1028         selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
1029         if (selectDebug) {
1030           filterFlag +=1;
1031           filterCut_Set3=kTRUE;
1032         }
1033       }
1034
1035       if (fTrackFilterTPC) {
1036
1037         selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
1038         if (selectDebug){//only tracks which pass the TPC-only track cuts
1039           trackmult++;
1040           filterFlag +=2;
1041           filterCut_Set1=kTRUE;
1042
1043         }
1044
1045       }
1046
1047       if (fTrackFilter) {
1048         selectDebug = fTrackFilter->IsSelected(esdTrack);
1049         if (selectDebug) {
1050           filterCut_Set2=kTRUE;
1051         }
1052       }
1053
1054
1055       if(filterFlag==0)
1056         continue;
1057
1058
1059       //
1060       // Track was accepted
1061       //
1062
1063       // Here we want to add histograms!
1064
1065       if (esdTrack->Pt() < fMinPt) {
1066
1067         // Keep small fraction of low pT tracks
1068         if(fRandom->Rndm() > fLowPtFraction)
1069           continue;
1070       } // else {
1071       // Here we want to add the high pt part of the histograms!
1072       //    }
1073
1074       Short_t charge  = esdTrack->Charge();
1075       Float_t pt      = esdTrack->Pt();
1076       Float_t p       = esdTrack->P();
1077       Float_t eta     = esdTrack->Eta();
1078       Float_t phi     = esdTrack->Phi();
1079       Short_t ncl     = esdTrack->GetTPCsignalN();
1080       Short_t neff    = Short_t(esdTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1081       //          Short_t nclf    = esdTrack->GetTPCNclsF();
1082       Float_t dedx    = esdTrack->GetTPCsignal();
1083       Float_t tpcchi  = 0;
1084       if(esdTrack->GetTPCNcls() > 0)
1085         tpcchi = esdTrack->GetTPCchi2()/Float_t(esdTrack->GetTPCNcls());
1086       Float_t b[2];
1087       Float_t bCov[3];
1088       esdTrack->GetImpactParameters(b,bCov);
1089       Float_t dcaxy   = b[0];
1090       Float_t dcaz    = b[1];
1091       Double_t p_con[3] = {0, 0, 0};
1092       esdTrack->GetConstrainedPxPyPz(p_con);
1093
1094
1095       //        Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1096       // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
1097       // Float_t pttpc   = tpcParam->Pt();
1098       // Float_t ptpc    = tpcParam->P();
1099
1100       Float_t ptMC        = 0;
1101       Short_t pidCode     = 0; // 0 = real data / no mc track!
1102       Short_t primaryFlag = 0; // 0 = real data / not primary mc track
1103       Int_t   pdgMother   = 0;
1104
1105
1106       //with Globaltrack parameters????
1107       if(fAnalysisMC) {
1108
1109         const Int_t label = TMath::Abs(esdTrack->GetLabel());
1110         TParticle* mcTrack = fMCStack->Particle(label);
1111         if (mcTrack){
1112
1113           if(fMCStack->IsPhysicalPrimary(label))
1114             primaryFlag = 1;
1115
1116           Int_t pdgCode = mcTrack->GetPdgCode();
1117           pidCode = GetPidCode(pdgCode);
1118
1119           ptMC      = mcTrack->Pt();
1120
1121           TParticle* mother = FindPrimaryMother(fMCStack, label);
1122           pdgMother = mother->GetPdgCode();
1123         }
1124       }
1125
1126
1127       //TOF
1128       Float_t beta = -99;
1129       if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
1130         if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
1131           beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
1132       }
1133
1134       if(fTreeOption) {
1135
1136         DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
1137         nadded++;
1138
1139         track->p          = p;
1140         track->pt         = pt;
1141         //        track->ptcon   = pt_con;
1142         //        track->tpcchi  = tpcchi;
1143         track->eta        = eta;
1144         track->phi        = phi;
1145         track->q          = charge;
1146         track->filter     = filterFlag;
1147         track->ncl        = ncl;
1148         track->neff       = neff;
1149         track->dedx       = dedx;
1150         track->beta       = beta;
1151         track->dcaxy      = dcaxy;
1152         track->dcaz       = dcaz;
1153         track->pid        = pidCode;
1154         track->primary    = primaryFlag;
1155         track->pttrue     = ptMC;
1156         track->mother     = pdgMother;
1157         track->filterset1 = filterCut_Set1;
1158         track->filterset2 = filterCut_Set2;
1159         track->filterset3 = filterCut_Set3;
1160
1161
1162       }
1163     }//end of track loop
1164
1165   }//end global
1166
1167   else if( analysisMode==kTPCTrk ){
1168     const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
1169     if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
1170
1171
1172     for(Int_t iT = 0; iT < nESDTracks; iT++) {
1173
1174       AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1175
1176       AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),esdTrack->GetID());
1177
1178       if(!trackTPC) continue;
1179
1180       UInt_t selectDebug = 0;
1181       selectDebug = fTrackFilterTPC->IsSelected(trackTPC);
1182       if(selectDebug==0) continue;
1183
1184       if(trackTPC->Pt()>0.){
1185         // only constrain tracks above threshold
1186         AliExternalTrackParam exParam;
1187         // take the B-field from the ESD, no 3D fieldMap available at this point
1188         Bool_t relate = false;
1189         relate = trackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1190                                            kVeryBig,&exParam);
1191         if(!relate){
1192           delete trackTPC;
1193           continue;
1194         }
1195         trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
1196                       exParam.GetCovariance());
1197       }
1198       else continue;
1199
1200
1201
1202
1203       if(TMath::Abs(trackTPC->Eta()) > fEtaCut)
1204         continue;
1205
1206       //
1207       // Track was accepted
1208       //
1209
1210       // Here we want to add histograms!
1211
1212       if (trackTPC->Pt() < fMinPt) {
1213
1214         // Keep small fraction of low pT tracks
1215         if(fRandom->Rndm() > fLowPtFraction)
1216           continue;
1217       } // else {
1218       // Here we want to add the high pt part of the histograms!
1219       //    }
1220
1221       Short_t charge  = trackTPC->Charge();
1222       Float_t pt      = trackTPC->Pt();
1223       Float_t p       = trackTPC->P();
1224       Float_t eta     = trackTPC->Eta();
1225       Float_t phi     = trackTPC->Phi();
1226       Short_t ncl     = trackTPC->GetTPCsignalN();
1227       Short_t neff    = Short_t(trackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1228       //          Short_t nclf    = esdTrack->GetTPCNclsF();
1229       Float_t dedx    = trackTPC->GetTPCsignal();
1230       Float_t tpcchi  = 0;
1231       if(trackTPC->GetTPCNcls() > 0)
1232         tpcchi = trackTPC->GetTPCchi2()/Float_t(trackTPC->GetTPCNcls());
1233       //Float_t b[2];
1234       //Float_t bCov[3];
1235       //trackTPC->GetImpactParameters(b,bCov);
1236       //Float_t dcaxy   = b[0];
1237       //Float_t dcaz    = b[1];
1238
1239       Float_t dcaxy = 0.;
1240       Float_t dcaz = 0.;
1241       trackTPC->GetImpactParameters(dcaxy,dcaz);
1242
1243
1244       Double_t p_con[3] = {0, 0, 0};
1245       trackTPC->GetConstrainedPxPyPz(p_con);
1246
1247
1248       // Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1249       // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
1250       // Float_t pttpc   = tpcParam->Pt();
1251       // Float_t ptpc    = tpcParam->P();
1252
1253       Float_t ptMC        = 0;
1254       Short_t pidCode     = 0; // 0 = real data / no mc track!
1255       Short_t primaryFlag = 0; // 0 = real data / not primary mc track
1256       Int_t   pdgMother   = 0;
1257
1258
1259       //with Globaltrack parameters????
1260       if(fAnalysisMC) {
1261
1262         const Int_t label = TMath::Abs(esdTrack->GetLabel());
1263         TParticle* mcTrack = fMCStack->Particle(label);
1264         if (mcTrack){
1265
1266           if(fMCStack->IsPhysicalPrimary(label))
1267             primaryFlag = 1;
1268
1269           Int_t pdgCode = mcTrack->GetPdgCode();
1270           pidCode = GetPidCode(pdgCode);
1271
1272           ptMC      = mcTrack->Pt();
1273
1274           TParticle* mother = FindPrimaryMother(fMCStack, label);
1275           pdgMother = mother->GetPdgCode();
1276         }
1277       }
1278
1279
1280       //TOF
1281       Float_t beta = -99;
1282       if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
1283         if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
1284           beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
1285       }
1286
1287       if(fTreeOption) {
1288
1289         DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
1290         nadded++;
1291
1292         tracktpc->p          = p;
1293         tracktpc->pt         = pt;
1294         //        track->ptcon   = pt_con;
1295         //        track->tpcchi  = tpcchi;
1296         tracktpc->eta        = eta;
1297         tracktpc->phi        = phi;
1298         tracktpc->q          = charge;
1299         tracktpc->filter     = 1;
1300         tracktpc->ncl        = ncl;
1301         tracktpc->neff       = neff;
1302         tracktpc->dedx       = dedx;
1303         tracktpc->beta       = beta;//computed with Global tracks
1304         tracktpc->dcaxy      = dcaxy;
1305         tracktpc->dcaz       = dcaz;
1306         tracktpc->pid        = pidCode;
1307         tracktpc->primary    = primaryFlag;
1308         tracktpc->pttrue     = ptMC;
1309         tracktpc->mother     = pdgMother;
1310         tracktpc->filterset1 = 0;
1311         tracktpc->filterset2 = 0;
1312         tracktpc->filterset3 = 0;
1313
1314       }
1315
1316
1317       if(trackTPC)
1318         delete trackTPC;
1319
1320
1321     }//end of track loop
1322
1323   }//end of: if isglobal==kFALSE
1324
1325
1326   if(fTreeOption) {
1327
1328     if( analysisMode==kGlobalTrk ){
1329       Sort(fTrackArrayGlobalPar, kFALSE);
1330
1331       fEvent->trackmult = trackmult;
1332       fEvent->n         = nadded;
1333     }
1334     else if( analysisMode==kTPCTrk ){
1335       Sort(fTrackArrayTPCPar, kFALSE);
1336     }
1337
1338   }
1339
1340
1341 }
1342 //__________________________________________________________________
1343 void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
1344
1345
1346   Int_t     trackmult = 0; // no pt cuts
1347   Int_t     nadded    = 0;
1348   Int_t nAODTracks = AODevent->GetNumberOfTracks();
1349   if( analysisMode == kGlobalTrk ){
1350
1351     if(fTrackArrayGlobalPar)
1352       fTrackArrayGlobalPar->Clear();
1353
1354   }
1355   if( analysisMode == kTPCTrk ){
1356     if(fTrackArrayTPCPar)
1357       fTrackArrayTPCPar->Clear();
1358   }
1359
1360
1361
1362
1363   //const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1364   //if( vertexSPD->GetNContributors() < 1 || TMath::Abs( vertexSPD->GetZ()) > 10.0 ) return;
1365
1366
1367
1368   if( analysisMode==kGlobalTrk ){
1369
1370
1371      const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1372     for(Int_t iT = 0; iT < nAODTracks; iT++) {
1373
1374       AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1375
1376
1377
1378       UShort_t filterFlag = 0;
1379       Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1380       Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1381       Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1382
1383
1384
1385       if (fTrackFilterGolden) {
1386
1387         // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1388         if(aodTrack->TestFilterBit(32)) {
1389           filterFlag +=1;
1390           filterCut_Set3 = kTRUE;
1391         }
1392       }
1393
1394
1395       if (fTrackFilterTPC) {
1396         // TPC only cuts is bit 1 according to above macro
1397         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
1398         if(aodTrack->TestFilterBit(1)){
1399           filterFlag +=2;
1400           filterCut_Set1 = kTRUE;
1401           trackmult++;
1402
1403         }
1404       }
1405
1406
1407
1408       if(filterFlag==0)
1409         continue;
1410
1411
1412       Double_t b[2], cov[3];
1413       if (!(aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov)))
1414         filterFlag = 32; // propagation failed!!!!!;
1415       Float_t dcaxy   = b[0];
1416       Float_t dcaz    = b[1];
1417
1418
1419
1420
1421
1422       // As I understand this routine recalculates the momentum so it should be called first!
1423       //Double_t b[2], cov[3];
1424
1425
1426       //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1427       //        filterFlag = 32; // propagation failed!!!!!
1428
1429       if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1430         continue;
1431       if (aodTrack->Pt() < fMinPt) {
1432
1433         // Keep small fraction of low pT tracks
1434         if(fRandom->Rndm() > fLowPtFraction)
1435           continue;
1436       } // else {
1437       // Here we want to add the high pt part of the histograms!
1438       //    }
1439
1440
1441
1442
1443       Short_t charge  = aodTrack->Charge();
1444       Float_t pt      = aodTrack->Pt();
1445       Float_t p       = aodTrack->P();
1446       Float_t eta     = aodTrack->Eta();
1447       Float_t phi     = aodTrack->Phi();
1448       AliAODPid* aodPid = aodTrack->GetDetPid();
1449       Short_t ncl     = -10;
1450       Short_t neff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1451       //          Short_t nclf    = aodTrack->GetTPCNclsF();
1452       Float_t dedx    = -10;
1453       Float_t beta = -99;
1454       if(aodPid) {
1455         ncl     = aodPid->GetTPCsignalN();
1456         dedx    = aodPid->GetTPCsignal();
1457         //TOF
1458         if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
1459           Double_t tof[5];
1460           aodPid->GetIntegratedTimes(tof);
1461           beta = tof[0]/aodPid->GetTOFsignal();
1462         }
1463       }
1464       //        Float_t tpcchi = aodTrack->Chi2perNDF();
1465
1466       // Double_t p_con[3] = {0, 0, 0};
1467       // aodTrack->GetConstrainedPxPyPz(p_con);
1468       //        Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1469       // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
1470       // Float_t pttpc   = tpcParam->Pt();
1471       // Float_t ptpc    = tpcParam->P();
1472
1473       Float_t ptMC        = 0;
1474       Short_t pidCode     = 0; // 0 = real data / no mc track!
1475       Short_t primaryFlag = 0; // 0 = real data / not primary mc track
1476       Int_t   pdgMother   = 0;
1477
1478       if(fAnalysisMC) {
1479
1480         const Int_t label = TMath::Abs(aodTrack->GetLabel());
1481
1482         AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1483         if (mcTrack){
1484
1485           if(mcTrack->IsPhysicalPrimary())
1486             primaryFlag = 1;
1487
1488           Int_t pdgCode = mcTrack->GetPdgCode();
1489           pidCode = GetPidCode(pdgCode);
1490
1491           ptMC      = mcTrack->Pt();
1492
1493           AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
1494           pdgMother = mother->GetPdgCode();
1495         }
1496       }
1497
1498
1499       if(fTreeOption) {
1500
1501         DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
1502         nadded++;
1503
1504         track->p          = p;
1505         track->pt         = pt;
1506         //        track->ptcon   = pt_con;
1507         //        track->tpcchi  = tpcchi;
1508         track->eta        = eta;
1509         track->phi        = phi;
1510         track->q          = charge;
1511         track->filter     = filterFlag;
1512         track->ncl        = ncl;
1513         track->neff       = neff;
1514         track->dedx       = dedx;
1515         track->beta       = beta;
1516         track->dcaxy      = dcaxy;
1517         track->dcaz       = dcaz;
1518         track->pid        = pidCode;
1519         track->primary    = primaryFlag;
1520         track->pttrue     = ptMC;
1521         track->mother     = pdgMother;
1522         track->filterset1 = filterCut_Set1;
1523         track->filterset2 = filterCut_Set2;
1524         track->filterset3 = filterCut_Set3;
1525       }
1526     }//end of track loop
1527
1528
1529
1530   }//end of global track analysis
1531
1532
1533
1534   else if( analysisMode==kTPCTrk ){
1535
1536     const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1537     if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
1538
1539
1540     //const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1541     for(Int_t iT = 0; iT < nAODTracks; iT++) {
1542
1543       AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1544
1545
1546
1547       UShort_t filterFlag = 0;
1548       Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1549       Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1550       Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1551
1552
1553
1554       // TPC only cuts is bit 1 according to above macro
1555       // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
1556       if(aodTrack->TestFilterBit(128)){
1557         filterFlag +=2;
1558         trackmult++;
1559       }
1560
1561
1562       if(filterFlag==0)
1563         continue;
1564
1565
1566       Double_t b[2], cov[3];
1567       //AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1568       //Double_t b[2] = {-99., -99.};
1569       //Double_t bCov[3] = {-99., -99., -99.};
1570       //aodTrack->PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov);
1571       if (!(aodTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov)))
1572         filterFlag = 32; // propagation failed!!!!!;
1573       Float_t dcaxy   = b[0];
1574       Float_t dcaz    = b[1];
1575
1576
1577
1578
1579       // As I understand this routine recalculates the momentum so it should be called first!
1580       //Double_t b[2], cov[3];
1581
1582
1583       //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1584       //        filterFlag = 32; // propagation failed!!!!!
1585
1586       if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1587         continue;
1588       if (aodTrack->Pt() < fMinPt) {
1589
1590         // Keep small fraction of low pT tracks
1591         if(fRandom->Rndm() > fLowPtFraction)
1592           continue;
1593       } // else {
1594       // Here we want to add the high pt part of the histograms!
1595       //    }
1596
1597
1598
1599
1600       Short_t charge  = aodTrack->Charge();
1601       Float_t pt      = aodTrack->Pt();
1602       Float_t p       = aodTrack->P();
1603       Float_t eta     = aodTrack->Eta();
1604       Float_t phi     = aodTrack->Phi();
1605       AliAODPid* aodPid = aodTrack->GetDetPid();
1606       Short_t ncl     = -10;
1607       Short_t neff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1608       //          Short_t nclf    = aodTrack->GetTPCNclsF();
1609       Float_t dedx    = -10;
1610       Float_t beta = -99;
1611       if(aodPid) {
1612         ncl     = aodPid->GetTPCsignalN();
1613         dedx    = aodPid->GetTPCsignal();
1614         //TOF
1615         if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
1616           Double_t tof[5];
1617           aodPid->GetIntegratedTimes(tof);
1618           beta = tof[0]/aodPid->GetTOFsignal();
1619         }
1620       }
1621       //        Float_t tpcchi = aodTrack->Chi2perNDF();
1622
1623       // Double_t p_con[3] = {0, 0, 0};
1624       // aodTrack->GetConstrainedPxPyPz(p_con);
1625       //        Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1626       // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
1627       // Float_t pttpc   = tpcParam->Pt();
1628       // Float_t ptpc    = tpcParam->P();
1629
1630       Float_t ptMC        = 0;
1631       Short_t pidCode     = 0; // 0 = real data / no mc track!
1632       Short_t primaryFlag = 0; // 0 = real data / not primary mc track
1633       Int_t   pdgMother   = 0;
1634
1635       if(fAnalysisMC) {
1636
1637         const Int_t label = TMath::Abs(aodTrack->GetLabel());
1638
1639         AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1640         if (mcTrack){
1641
1642           if(mcTrack->IsPhysicalPrimary())
1643             primaryFlag = 1;
1644
1645           Int_t pdgCode = mcTrack->GetPdgCode();
1646           pidCode = GetPidCode(pdgCode);
1647
1648           ptMC      = mcTrack->Pt();
1649
1650           AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
1651           pdgMother = mother->GetPdgCode();
1652         }
1653       }
1654
1655
1656       if(fTreeOption) {
1657
1658         DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
1659         nadded++;
1660
1661         tracktpc->p          = p;
1662         tracktpc->pt         = pt;
1663         //        track->ptcon   = pt_con;
1664         //        track->tpcchi  = tpcchi;
1665         tracktpc->eta        = eta;
1666         tracktpc->phi        = phi;
1667         tracktpc->q          = charge;
1668         tracktpc->filter     = filterFlag;
1669         tracktpc->ncl        = ncl;
1670         tracktpc->neff       = neff;
1671         tracktpc->dedx       = dedx;
1672         tracktpc->beta       = beta;
1673         tracktpc->dcaxy      = dcaxy;
1674         tracktpc->dcaz       = dcaz;
1675         tracktpc->pid        = pidCode;
1676         tracktpc->primary    = primaryFlag;
1677         tracktpc->pttrue     = ptMC;
1678         tracktpc->mother     = pdgMother;
1679         tracktpc->filterset1 = filterCut_Set1;
1680         tracktpc->filterset2 = filterCut_Set2;
1681         tracktpc->filterset3 = filterCut_Set3;
1682       }
1683     }//end of track loop
1684
1685
1686   }//end of global track analysis
1687
1688
1689   if(fTreeOption) {
1690
1691     if( analysisMode==kGlobalTrk ){
1692       Sort(fTrackArrayGlobalPar, kFALSE);
1693
1694       fEvent->trackmult = trackmult;
1695       fEvent->n         = nadded;
1696
1697
1698     }
1699     if( analysisMode==kTPCTrk ){
1700       Sort(fTrackArrayTPCPar, kFALSE);
1701     }
1702
1703   }
1704
1705 }