PWGLF/SPECTRA/IdentifiedHighPt/grid/AliAnalysisTaskHighPtDeDxV0.cxx

   1 #include "AliAnalysisTaskHighPtDeDxV0.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
  20 #include <AliMCEventHandler.h>
  21 #include <AliMCEvent.h>
  22 #include <AliStack.h>
  23
  24 #include <AliHeader.h>
  25 #include <AliGenPythiaEventHeader.h>
  26 #include <AliGenDPMjetEventHeader.h>
  27
  28 #include <AliAODVertex.h>
  29 #include <AliESDv0.h>
  30 #include <AliKFVertex.h>
  31
  32 #include "AliCentrality.h"
  33 #include "AliESDtrackCuts.h"
  34
  35 #include <AliAODTrack.h>
  36 #include <AliAODPid.h>
  37 #include <AliAODMCHeader.h>
  38
  39 // STL includes
  40 #include <iostream>
  41 using namespace std;
  42
  43 //_____________________________________________________________________________
  44 //
  45 // Responsible:
  46 // Peter Christiansen (Lund)
  47 //
  48 //_____________________________________________________________________________
  49
  50
  51 ClassImp(AliAnalysisTaskHighPtDeDxV0)
  52
  53 const Double_t AliAnalysisTaskHighPtDeDxV0::fgkClight = 2.99792458e-2;
  54
  55 //_____________________________________________________________________________
  56 AliAnalysisTaskHighPtDeDxV0::AliAnalysisTaskHighPtDeDxV0():
  57   AliAnalysisTaskSE(),
  58   fESD(0x0),
  59   fAOD(0x0),
  60   fMC(0x0),
  61   fMCStack(0x0),
  62   fMCArray(0x0),
  63   fTrackFilter(0x0),
  64   fTrackFilterGolden(0x0),
  65   fTrackFilterTPC(0x0),
  66   fAnalysisType("ESD"),
  67   fAnalysisMC(kFALSE),
  68   fAnalysisPbPb(kFALSE),
  69   ftrigBit1(0x0),
  70   ftrigBit2(0x0),
  71   fEvent(0x0),
  72   fV0ArrayGlobalPar(0x0),
  73   fV0ArrayTPCPar(0x0),
  74   fTrackArrayMC(0x0),
  75   fVtxCut(10.0),
  76   fEtaCut(0.9),
  77   fMinPt(0.1),
  78   fMinCent(0.0),
  79   fMaxCent(100.0),
  80   fMassCut(0.1),
  81   fTreeOption(0),
  82   fRequireRecV0(kTRUE),
  83   fStoreMcIn(kFALSE),
  84   fMcProcessType(-999),
  85   fTriggeredEventMB(-999),
  86   fVtxStatus(-999),
  87   fZvtx(-999),
  88   fZvtxMC(-999),
  89   fRun(-999),
  90   fEventId(-999),
  91   fListOfObjects(0),
  92   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
  93   fTree(0x0),
  94   fnv0(0)
  95 {
  96   // Default constructor (should not be used)
  97
  98 }
  99
 100 //______________________________________________________________________________
 101 AliAnalysisTaskHighPtDeDxV0::AliAnalysisTaskHighPtDeDxV0(const char *name):
 102   AliAnalysisTaskSE(name),
 103   fESD(0x0),
 104   fAOD(0x0),
 105   fMC(0x0),
 106   fMCStack(0x0),
 107   fMCArray(0x0),
 108   fTrackFilter(0x0),
 109   fTrackFilterGolden(0x0),
 110   fTrackFilterTPC(0x0),
 111   fAnalysisType("ESD"),
 112   fAnalysisMC(kFALSE),
 113   fAnalysisPbPb(kFALSE),
 114   ftrigBit1(0x0),
 115   ftrigBit2(0x0),
 116   fEvent(0x0),
 117   fV0ArrayGlobalPar(0x0),
 118   fV0ArrayTPCPar(0x0),
 119   fTrackArrayMC(0x0),
 120   fVtxCut(10.0),
 121   fEtaCut(0.9),
 122   fMinPt(0.1),
 123   fMinCent(0.0),
 124   fMaxCent(100.0),
 125   fMassCut(0.1),
 126   fTreeOption(0),
 127   fRequireRecV0(kTRUE),
 128   fStoreMcIn(kFALSE),
 129   fMcProcessType(-999),
 130   fTriggeredEventMB(-999),
 131   fVtxStatus(-999),
 132   fZvtx(-999),
 133   fZvtxMC(-999),
 134   fRun(-999),
 135   fEventId(-999),
 136   fListOfObjects(0),
 137   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
 138   fTree(0x0),
 139   fnv0(0)
 140 {
 141   // Output slot #1 writes into a TList
 142   DefineOutput(1, TList::Class());
 143 }
 144
 145 //_____________________________________________________________________________
 146 AliAnalysisTaskHighPtDeDxV0::~AliAnalysisTaskHighPtDeDxV0()
 147 {
 148   // Destructor
 149   // histograms are in the output list and deleted when the output
 150   // list is deleted by the TSelector dtor
 151   if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 152     delete fListOfObjects;
 153     fListOfObjects = 0;
 154   }
 155
 156   // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse
 157   // if (fListOfObjects  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
 158
 159
 160
 161 }
 162
 163 //______________________________________________________________________________
 164 void AliAnalysisTaskHighPtDeDxV0::UserCreateOutputObjects()
 165 {
 166   // This method is called once per worker node
 167   // Here we define the output: histograms and debug tree if requested
 168
 169   OpenFile(1);
 170   fListOfObjects = new TList();
 171   fListOfObjects->SetOwner();
 172
 173   //
 174   // Histograms
 175   //
 176   fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 1, 0, 1);
 177   fListOfObjects->Add(fEvents);
 178
 179   fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
 180   fListOfObjects->Add(fVtx);
 181
 182   fnv0 = new TH1F("fnv0","# V0's", 10000,0,10000);
 183   fListOfObjects->Add(fnv0);
 184
 185   fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
 186   fListOfObjects->Add(fVtxBeforeCuts);
 187
 188   fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
 189   fListOfObjects->Add(fVtxAfterCuts);
 190
 191   if (fAnalysisMC) {
 192     fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
 193     fListOfObjects->Add(fVtxMC);
 194
 195   }
 196
 197   if (fTreeOption) {
 198
 199     fTree = new TTree("tree","Event data");
 200     fEvent = new DeDxEvent();
 201     fTree->Branch("event", &fEvent);
 202
 203     fV0ArrayGlobalPar = new TClonesArray("DeDxV0", 1000);
 204     fTree->Bronch("v0GlobalPar", "TClonesArray", &fV0ArrayGlobalPar);
 205
 206     fV0ArrayTPCPar = new TClonesArray("DeDxV0", 1000);
 207     fTree->Bronch("v0TPCPar", "TClonesArray", &fV0ArrayTPCPar);
 208
 209     if (fAnalysisMC && fStoreMcIn) {
 210       fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
 211       fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
 212     }
 213
 214     fTree->SetDirectory(0);
 215
 216     fListOfObjects->Add(fTree);
 217
 218   }
 219
 220   // Post output data.
 221   PostData(1, fListOfObjects);
 222 }
 223
 224 //______________________________________________________________________________
 225 void AliAnalysisTaskHighPtDeDxV0::UserExec(Option_t *)
 226 {
 227   // Main loop
 228
 229   //
 230   // Comment: This method matches completely the same method for the high pT
 231   // tracks
 232   //
 233
 234
 235   //
 236   // First we make sure that we have valid input(s)!
 237   //
 238   if (fAnalysisType == "ESD"){
 239     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 240     if(!fESD){
 241       Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 242       this->Dump();
 243       return;
 244     }
 245   } else {
 246     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 247     if(!fAOD){
 248       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 249       this->Dump();
 250       return;
 251     }
 252   }
 253
 254   if (fAnalysisMC) {
 255
 256     if (fAnalysisType == "ESD"){
 257       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 258       if(!fMC){
 259         Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 260         this->Dump();
 261         return;
 262       }
 263
 264       fMCStack = fMC->Stack();
 265
 266       if(!fMCStack){
 267         Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
 268         this->Dump();
 269         return;
 270       }
 271     } else { // AOD
 272
 273       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 274       if(fMC)
 275         fMC->Dump();
 276
 277       fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
 278       if(!fMCArray){
 279         Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
 280         this->Dump();
 281         return;
 282       }
 283     }
 284   }
 285
 286     // Get trigger decision
 287   fTriggeredEventMB = 0; //init
 288
 289
 290
 291   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 292      ->IsEventSelected() & ftrigBit1 ){
 293     fTriggeredEventMB = 1;  //event triggered as minimum bias
 294   }
 295   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 296      ->IsEventSelected() & ftrigBit2 ){
 297     // From AliVEvent:
 298     //    kINT7         = BIT(1), // V0AND trigger, offline V0 selection
 299     fTriggeredEventMB += 2;
 300   }
 301
 302
 303   // Get process type for MC
 304   fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
 305
 306   // real data that are not triggered we skip
 307   if(!fAnalysisMC && !fTriggeredEventMB)
 308     return;
 309
 310   if (fAnalysisMC) {
 311
 312     if (fAnalysisType == "ESD"){
 313
 314       AliHeader* headerMC = fMC->Header();
 315       if (headerMC) {
 316
 317         AliGenEventHeader* genHeader = headerMC->GenEventHeader();
 318         TArrayF vtxMC(3); // primary vertex  MC
 319         vtxMC[0]=9999; vtxMC[1]=9999;  vtxMC[2]=9999; //initialize with dummy
 320         if (genHeader) {
 321           genHeader->PrimaryVertex(vtxMC);
 322         }
 323         fZvtxMC = vtxMC[2];
 324
 325         // PYTHIA:
 326         AliGenPythiaEventHeader* pythiaGenHeader =
 327           dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
 328         if (pythiaGenHeader) {  //works only for pythia
 329           fMcProcessType =  GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
 330         }
 331         // PHOJET:
 332         AliGenDPMjetEventHeader* dpmJetGenHeader =
 333           dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
 334         if (dpmJetGenHeader) {
 335           fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
 336         }
 337       }
 338     } else { // AOD
 339
 340       AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
 341       if(mcHeader) {
 342         fZvtxMC = mcHeader->GetVtxZ();
 343
 344         if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
 345           fMcProcessType =  GetPythiaEventProcessType(mcHeader->GetEventType());
 346         } else {
 347           fMcProcessType =  GetDPMjetEventProcessType(mcHeader->GetEventType());
 348         }
 349       }
 350     }
 351   }
 352
 353   // There are 3 cases
 354   // Vertex: NO  - status -1
 355   // Vertex: YES : outside cut - status 0
 356   //             : inside cut  - status 1
 357   // We have to be careful how we normalize because we probably want to
 358   // normalize to:
 359   // Nevents=(No vertex + outside + inside)/(out + in)*in
 360
 361   if (fAnalysisType == "ESD")
 362     fZvtx = GetVertex(fESD);
 363   else // AOD
 364     fZvtx = GetVertex(fAOD);
 365
 366   fVtxStatus = -999;
 367
 368   if(fZvtx<-990) {
 369
 370     fVtxStatus = -1;
 371     if(fTriggeredEventMB)
 372       fVtx->Fill(0);
 373     if(fAnalysisMC)
 374       fVtxMC->Fill(0);
 375   } else {
 376
 377     if(fTriggeredEventMB)
 378       fVtx->Fill(1);
 379     if(fAnalysisMC)
 380       fVtxMC->Fill(1);
 381     fVtxBeforeCuts->Fill(fZvtx);
 382     fVtxStatus = 0;
 383     if (TMath::Abs(fZvtx) < fVtxCut) {
 384       fVtxAfterCuts->Fill(fZvtx);
 385       fVtxStatus = 1;
 386     }
 387   }
 388
 389   // store MC event data no matter what
 390   if(fAnalysisMC && fStoreMcIn) {
 391
 392     if (fAnalysisType == "ESD") {
 393       ProcessMCTruthESD();
 394     } else { // AOD
 395       ProcessMCTruthAOD();
 396     }
 397   }
 398   Float_t centrality = -10;
 399   // only analyze triggered events
 400   if(fTriggeredEventMB) {
 401
 402     if (fAnalysisType == "ESD"){
 403       if(fAnalysisPbPb){
 404         AliCentrality *centObject = fESD->GetCentrality();
 405         centrality = centObject->GetCentralityPercentile("V0M");
 406         if((centrality>fMaxCent)||(centrality<fMinCent))return;
 407       }
 408       Int_t nv0test=fESD->GetNumberOfV0s();
 409       cout<<"&&&&&&&&&&&&&&&&   hello world"<<endl;
 410       fnv0->Fill(nv0test);
 411       AnalyzeESD(fESD);
 412
 413     } else { // AOD
 414       if(fAnalysisPbPb){
 415         AliCentrality *centObject = fAOD->GetCentrality();
 416         centrality = centObject->GetCentralityPercentile("V0M");
 417         if((centrality>fMaxCent)||(centrality<fMinCent))return;
 418       }
 419       Int_t nv0test=fAOD->GetNumberOfV0s();
 420       fnv0->Fill(nv0test);
 421       AnalyzeAOD(fAOD);
 422
 423     }
 424   }
 425
 426
 427
 428   if( fTreeOption) {
 429
 430     fEvent->process = fMcProcessType;
 431     fEvent->trig    = fTriggeredEventMB;
 432     fEvent->zvtxMC  = fZvtxMC;
 433     fEvent->cent      = centrality;
 434
 435     if(!fRequireRecV0 || fEvent->n > 0) // only fill if there are accepted V0s
 436       fTree->Fill();
 437     fV0ArrayGlobalPar->Clear();
 438     fV0ArrayTPCPar->Clear();
 439     if (fAnalysisMC && fStoreMcIn)
 440       fTrackArrayMC->Clear();
 441   }
 442
 443   // Post output data.
 444   PostData(1, fListOfObjects);
 445 }
 446
 447 //________________________________________________________________________
 448 void AliAnalysisTaskHighPtDeDxV0::AnalyzeESD(AliESDEvent* esdEvent)
 449 {
 450   fRun  = esdEvent->GetRunNumber();
 451   fEventId = 0;
 452   if(esdEvent->GetHeader())
 453     fEventId = GetEventIdAsLong(esdEvent->GetHeader());
 454
 455   //  Int_t     event     = esdEvent->GetEventNumberInFile();
 456   UInt_t    time      = esdEvent->GetTimeStamp();
 457   //  ULong64_t trigger   = esdEvent->GetTriggerMask();
 458   Float_t   magf      = esdEvent->GetMagneticField();
 459
 460
 461
 462   Short_t isPileup = esdEvent->IsPileupFromSPD();
 463
 464   // centrality
 465   Float_t centrality = 120;
 466   AliCentrality *centObject = esdEvent->GetCentrality();
 467   if(centObject)
 468     centrality = centObject->GetCentralityPercentile("V0M");
 469
 470   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
 471
 472     if(fVtxStatus==1) { // accepted vertex
 473
 474       // accepted event
 475       fEvents->Fill(0);
 476
 477       ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
 478       ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
 479
 480
 481     } // end if vtx status
 482   } // end if triggered
 483
 484   if(fTreeOption) {
 485
 486     fEvent->run       = fRun;
 487     fEvent->eventid   = fEventId;
 488     fEvent->time      = time;
 489     //fEvent->cent      = centrality;
 490     fEvent->mag       = magf;
 491     fEvent->zvtx      = fZvtx;
 492     fEvent->vtxstatus = fVtxStatus;
 493     //fEvent->trackmult = trackmult;
 494     //fEvent->n         = nadded;
 495     fEvent->pileup    = isPileup;
 496   }
 497 }
 498
 499 //________________________________________________________________________
 500 void AliAnalysisTaskHighPtDeDxV0::AnalyzeAOD(AliAODEvent* aodEvent)
 501 {
 502   fRun  = aodEvent->GetRunNumber();
 503   fEventId = 0;
 504   if(aodEvent->GetHeader())
 505     fEventId = GetEventIdAsLong(aodEvent->GetHeader());
 506
 507   //  Int_t     event     = aodEvent->GetEventNumberInFile();
 508   UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
 509   //  ULong64_t trigger   = aodEvent->GetTriggerMask();
 510   Float_t   magf      = aodEvent->GetMagneticField();
 511
 512   Int_t     trackmult = 0; // no pt cuts
 513   Int_t     nadded    = 0;
 514
 515   Short_t isPileup = aodEvent->IsPileupFromSPD();
 516
 517   // centrality
 518   Float_t centrality = 120;
 519   AliCentrality *centObject = aodEvent->GetCentrality();
 520   if(centObject)
 521     centrality = centObject->GetCentralityPercentile("V0M");
 522
 523   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
 524
 525     if(fVtxStatus==1) { // accepted vertex
 526
 527       // accepted event
 528       fEvents->Fill(0);
 529
 530       ProduceArrayTrksAOD( aodEvent, kGlobalTrk );//produce array with global track parameters
 531       ProduceArrayTrksAOD( aodEvent, kTPCTrk );//array with tpc track parametes
 532
 533
 534     } // end if vtx status
 535   } // end if triggered
 536
 537   if(fTreeOption) {
 538
 539     fEvent->run       = fRun;
 540     fEvent->eventid   = fEventId;
 541     fEvent->time      = time;
 542     //fEvent->cent      = centrality;
 543     fEvent->mag       = magf;
 544     fEvent->zvtx      = fZvtx;
 545     fEvent->vtxstatus = fVtxStatus;
 546     //fEvent->trackmult = trackmult;
 547     //fEvent->n         = nadded;
 548     fEvent->pileup    = isPileup;
 549   }
 550
 551
 552
 553
 554
 555
 556
 557 }
 558
 559 //_____________________________________________________________________________
 560 Float_t AliAnalysisTaskHighPtDeDxV0::GetVertex(const AliVEvent* event) const
 561 {
 562   Float_t zvtx = -999;
 563
 564   const AliVVertex* primaryVertex = event->GetPrimaryVertex();
 565
 566   if(primaryVertex->GetNContributors()>0)
 567     zvtx = primaryVertex->GetZ();
 568
 569   return zvtx;
 570 }
 571
 572 //_____________________________________________________________________________
 573 Short_t AliAnalysisTaskHighPtDeDxV0::GetPidCode(Int_t pdgCode) const
 574 {
 575   // return our internal code for pions, kaons, and protons
 576
 577   Short_t pidCode = 6;
 578
 579   switch (TMath::Abs(pdgCode)) {
 580   case 211:
 581     pidCode = 1; // pion
 582     break;
 583   case 321:
 584     pidCode = 2; // kaon
 585     break;
 586   case 2212:
 587     pidCode = 3; // proton
 588     break;
 589   case 11:
 590     pidCode = 4; // electron
 591     break;
 592   case 13:
 593     pidCode = 5; // proton
 594     break;
 595   default:
 596     pidCode = 6;  // something else?
 597   };
 598
 599   return pidCode;
 600 }
 601
 602
 603 //_____________________________________________________________________________
 604 void AliAnalysisTaskHighPtDeDxV0::ProcessMCTruthESD()
 605 {
 606   // Fill the special MC histogram with the MC truth info
 607
 608   Short_t trackmult = 0;
 609   Short_t nadded    = 0;
 610   const Int_t nTracksMC = fMCStack->GetNtrack();
 611
 612   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
 613
 614     //Cuts
 615     if(!(fMCStack->IsPhysicalPrimary(iTracks)))
 616       continue;
 617
 618     TParticle* trackMC = fMCStack->Particle(iTracks);
 619
 620     Double_t chargeMC = trackMC->GetPDG()->Charge();
 621     if (chargeMC != 0) // select charge 0 particles!
 622       continue;
 623
 624     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
 625       continue;
 626
 627     trackmult++;
 628
 629     //   "charge:pt:p:eta:phi:pidCode"
 630     Float_t ptMC      = trackMC->Pt();
 631     Float_t pMC       = trackMC->P();
 632     Float_t etaMC     = trackMC->Eta();
 633     Float_t phiMC     = trackMC->Phi();
 634
 635     Int_t pdgCode = trackMC->GetPdgCode();
 636     if(!(pdgCode==310 || pdgCode == 3122 || pdgCode == -3122))
 637       continue;
 638     Short_t pidCodeMC = 0;
 639     pidCodeMC = GetPidCode(pdgCode);
 640
 641     // Warning: In the data code we cut on the daughters.....
 642     if (trackMC->Pt() < fMinPt)
 643       continue;
 644
 645     if(fTreeOption) {
 646
 647       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
 648       nadded++;
 649
 650       track->pMC   = pMC;
 651       track->ptMC  = ptMC;
 652       track->etaMC = etaMC;
 653       track->phiMC = phiMC;
 654       track->qMC   = Short_t(chargeMC);
 655       track->pidMC = pidCodeMC;
 656       track->pdgMC = pdgCode;
 657     }
 658
 659   }//MC track loop
 660
 661   if(fTreeOption) {
 662
 663     Sort(fTrackArrayMC, kTRUE);
 664
 665     fEvent->trackmultMC = trackmult;
 666     fEvent->nMC         = nadded;
 667   }
 668
 669 }
 670
 671 //_____________________________________________________________________________
 672 void AliAnalysisTaskHighPtDeDxV0::ProcessMCTruthAOD()
 673 {
 674   // Fill the special MC histogram with the MC truth info
 675
 676   Short_t trackmult = 0;
 677   Short_t nadded    = 0;
 678   const Int_t nTracksMC = fMCArray->GetEntriesFast();
 679
 680   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
 681
 682     AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
 683
 684     //Cuts
 685     if(!(trackMC->IsPhysicalPrimary()))
 686       continue;
 687
 688
 689     Double_t chargeMC = trackMC->Charge();
 690     if (chargeMC != 0) // Keep the neutral particles
 691       continue;
 692
 693     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
 694       continue;
 695
 696     trackmult++;
 697
 698     //   "charge:pt:p:eta:phi:pidCode"
 699     Float_t ptMC      = trackMC->Pt();
 700     Float_t pMC       = trackMC->P();
 701     Float_t etaMC     = trackMC->Eta();
 702     Float_t phiMC     = trackMC->Phi();
 703
 704     Int_t pdgCode = trackMC->PdgCode();
 705     if(!(pdgCode==310 || pdgCode == 3122 || pdgCode == -3122))
 706       continue;
 707     Short_t pidCodeMC = 0;
 708     pidCodeMC = GetPidCode(pdgCode);
 709
 710     if (trackMC->Pt() < fMinPt)
 711       continue;
 712
 713     if(fTreeOption) {
 714
 715       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
 716       nadded++;
 717
 718       track->pMC   = pMC;
 719       track->ptMC  = ptMC;
 720       track->etaMC = etaMC;
 721       track->phiMC = phiMC;
 722       track->qMC   = Short_t(chargeMC);
 723       track->pidMC = pidCodeMC;
 724       track->pdgMC = pdgCode;
 725     }
 726
 727   }//MC track loop
 728
 729   if(fTreeOption) {
 730
 731     Sort(fTrackArrayMC, kTRUE);
 732
 733     fEvent->trackmultMC = trackmult;
 734     fEvent->nMC         = nadded;
 735   }
 736
 737 }
 738
 739 //_____________________________________________________________________________
 740 Short_t AliAnalysisTaskHighPtDeDxV0::GetPythiaEventProcessType(Int_t pythiaType) {
 741   //
 742   // Get the process type of the event.  PYTHIA
 743   //
 744   // source PWG0   dNdpt
 745
 746   Short_t globalType = -1; //init
 747
 748   if(pythiaType==92||pythiaType==93){
 749     globalType = 2; //single diffractive
 750   }
 751   else if(pythiaType==94){
 752     globalType = 3; //double diffractive
 753   }
 754   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
 755   else {
 756     globalType = 1;  //non diffractive
 757   }
 758   return globalType;
 759 }
 760
 761 //_____________________________________________________________________________
 762 Short_t AliAnalysisTaskHighPtDeDxV0::GetDPMjetEventProcessType(Int_t dpmJetType) {
 763   //
 764   // get the process type of the event.  PHOJET
 765   //
 766   //source PWG0   dNdpt
 767   // can only read pythia headers, either directly or from cocktalil header
 768   Short_t globalType = -1;
 769
 770   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
 771     globalType = 1;
 772   }
 773   else if (dpmJetType==5 || dpmJetType==6) {
 774     globalType = 2;
 775   }
 776   else if (dpmJetType==7) {
 777     globalType = 3;
 778   }
 779   return globalType;
 780 }
 781
 782 //_____________________________________________________________________________
 783 ULong64_t AliAnalysisTaskHighPtDeDxV0::GetEventIdAsLong(AliVHeader* header) const
 784 {
 785   // To have a unique id for each event in a run!
 786   // Modified from AliRawReader.h
 787   return ((ULong64_t)header->GetBunchCrossNumber()+
 788           (ULong64_t)header->GetOrbitNumber()*3564+
 789           (ULong64_t)header->GetPeriodNumber()*16777215*3564);
 790 }
 791
 792
 793 //____________________________________________________________________
 794 TParticle* AliAnalysisTaskHighPtDeDxV0::FindPrimaryMother(AliStack* stack, Int_t label)
 795 {
 796   //
 797   // Finds the first mother among the primary particles of the particle identified by <label>,
 798   // i.e. the primary that "caused" this particle
 799   //
 800   // Taken from AliPWG0Helper class
 801   //
 802
 803   Int_t nSteps = 0;
 804
 805   Int_t motherLabel = FindPrimaryMotherLabel(stack, label, nSteps);
 806   if (motherLabel < 0)
 807     return 0;
 808
 809   return stack->Particle(motherLabel);
 810 }
 811
 812 //____________________________________________________________________
 813 Int_t AliAnalysisTaskHighPtDeDxV0::FindPrimaryMotherLabel(AliStack* stack, Int_t label, Int_t& nSteps)
 814 {
 815   //
 816   // Finds the first mother among the primary particles of the particle identified by <label>,
 817   // i.e. the primary that "caused" this particle
 818   //
 819   // returns its label
 820   //
 821   // Taken from AliPWG0Helper class
 822   //
 823   nSteps = 0;
 824   const Int_t nPrim  = stack->GetNprimary();
 825
 826   while (label >= nPrim) {
 827
 828     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
 829
 830     nSteps++; // 1 level down
 831
 832     TParticle* particle = stack->Particle(label);
 833     if (!particle) {
 834
 835       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
 836       return -1;
 837     }
 838
 839     // find mother
 840     if (particle->GetMother(0) < 0) {
 841
 842       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
 843       return -1;
 844     }
 845
 846     label = particle->GetMother(0);
 847   }
 848
 849   return label;
 850 }
 851
 852 //____________________________________________________________________
 853 AliAODMCParticle* AliAnalysisTaskHighPtDeDxV0::FindPrimaryMotherAOD(AliAODMCParticle* startParticle, Int_t& nSteps)
 854 {
 855   //
 856   // Finds the first mother among the primary particles of the particle identified by <label>,
 857   // i.e. the primary that "caused" this particle
 858   //
 859   // Taken from AliPWG0Helper class
 860   //
 861
 862   nSteps = 0;
 863
 864   AliAODMCParticle* mcPart = startParticle;
 865
 866   while (mcPart)
 867     {
 868
 869       if(mcPart->IsPrimary())
 870         return mcPart;
 871
 872       Int_t mother = mcPart->GetMother();
 873
 874       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
 875       nSteps++; // 1 level down
 876     }
 877
 878   return 0;
 879 }
 880
 881 //_____________________________________________________________________________
 882 void AliAnalysisTaskHighPtDeDxV0::Sort(TClonesArray* array, Bool_t isMC)
 883 {
 884   const Int_t n = array->GetEntriesFast();
 885   if(n==0) {
 886
 887     if(isMC)
 888       fEvent->ptmaxMC = 0;
 889     else
 890       fEvent->ptmax   = 0;
 891
 892     return;
 893   }
 894
 895   Float_t ptArray[n];
 896   Int_t   indexArray[n];
 897
 898   for(Int_t i = 0; i < n; i++) {
 899
 900     Float_t pt = 0;
 901     if(isMC) {
 902       DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
 903       pt = track->ptMC;
 904     } else {
 905       DeDxTrack* track = (DeDxTrack*)array->At(i);
 906       pt = track->pt;
 907     }
 908     ptArray[i]    = pt;
 909     indexArray[i] = i;
 910   }
 911
 912   TMath::Sort(n, ptArray, indexArray);
 913
 914   // set max pt
 915   if(isMC) {
 916     DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
 917     fEvent->ptmaxMC = track->ptMC;
 918   } else {
 919     DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
 920     fEvent->ptmax   = track->pt;
 921   }
 922
 923   // set order of each track
 924   for(Int_t i = 0; i < n; i++) {
 925
 926     if(isMC) {
 927       DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
 928       track->orderMC = i;
 929     } else {
 930       DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
 931       track->order = i;
 932     }
 933   }
 934 }
 935 //__________________________________________________________________
 936 void AliAnalysisTaskHighPtDeDxV0::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
 937   Int_t nv0s = ESDevent->GetNumberOfV0s();
 938   if(nv0s<1)return;
 939   Int_t     trackmult = 0; // no pt cuts
 940   Int_t     nadded    = 0;
 941   const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
 942   if (!myBestPrimaryVertex) return;
 943   if (!(myBestPrimaryVertex->GetStatus())) return;
 944   Double_t  lPrimaryVtxPosition[3];
 945   myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
 946
 947   Double_t  lPrimaryVtxCov[6];
 948   myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
 949   Double_t  lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
 950
 951   AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
 952
 953
 954   if( analysisMode == kGlobalTrk ){
 955     if(fV0ArrayGlobalPar)
 956       fV0ArrayGlobalPar->Clear();
 957   } else if( analysisMode == kTPCTrk ){
 958     if(fV0ArrayTPCPar)
 959       fV0ArrayTPCPar->Clear();
 960   }
 961
 962   if( analysisMode==kGlobalTrk ){
 963
 964
 965
 966
 967     // ################################
 968     // #### BEGINNING OF V0 CODE ######
 969     // ################################
 970
 971
 972     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
 973
 974       // This is the begining of the V0 loop
 975       AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
 976       if (!esdV0) continue;
 977
 978       // AliESDTrack (V0 Daughters)
 979       UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
 980       UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
 981
 982       AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
 983       AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
 984       if (!pTrack || !nTrack) {
 985         Printf("ERROR: Could not retreive one of the daughter track");
 986         continue;
 987       }
 988
 989       // Remove like-sign
 990       if (pTrack->GetSign() == nTrack->GetSign()) {
 991         //cout<< "like sign, continue"<< endl;
 992         continue;
 993       }
 994
 995       // Eta cut on decay products
 996       if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
 997         continue;
 998
 999       // Pt cut on decay products
1000       if (esdV0->Pt() < fMinPt)
1001         //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1002         continue;
1003
1004       //filter for positive track
1005       UShort_t filterFlag_p = 0;
1006       Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1007       Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1008       Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1009
1010       UInt_t selectDebug_p = 0;
1011       if (fTrackFilterGolden) {
1012         selectDebug_p = fTrackFilterGolden->IsSelected(pTrack);
1013         if (selectDebug_p) {
1014           filterFlag_p +=1;
1015           filterCut_Set3_p=kTRUE;
1016         }
1017       }
1018
1019       if (fTrackFilterTPC) {
1020
1021         selectDebug_p = fTrackFilterTPC->IsSelected(pTrack);
1022         if (selectDebug_p){//only tracks which pass the TPC-only track cuts
1023           filterFlag_p +=2;
1024           filterCut_Set1_p=kTRUE;
1025
1026         }
1027
1028       }
1029
1030       if (fTrackFilter) {
1031         selectDebug_p = fTrackFilter->IsSelected(pTrack);
1032         if (selectDebug_p) {
1033           filterCut_Set2_p=kTRUE;
1034         }
1035       }
1036
1037       if(filterFlag_p ==0 )
1038         continue;
1039
1040
1041       //filter for negative track
1042       UShort_t filterFlag_n = 0;
1043       Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1044       Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1045       Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1046
1047       UInt_t selectDebug_n = 0;
1048       if (fTrackFilterGolden) {
1049         selectDebug_n = fTrackFilterGolden->IsSelected(nTrack);
1050         if (selectDebug_n) {
1051           filterFlag_n +=1;
1052           filterCut_Set3_n=kTRUE;
1053         }
1054       }
1055
1056       if (fTrackFilterTPC) {
1057
1058         selectDebug_n = fTrackFilterTPC->IsSelected(nTrack);
1059         if (selectDebug_n){//only tracks which pass the TPC-only track cuts
1060           filterFlag_n +=2;
1061           filterCut_Set1_n=kTRUE;
1062
1063         }
1064
1065       }
1066
1067       if (fTrackFilter) {
1068         selectDebug_n = fTrackFilter->IsSelected(nTrack);
1069         if (selectDebug_n) {
1070           filterCut_Set2_n=kTRUE;
1071         }
1072       }
1073
1074
1075       // Check if switch does anything!
1076       Bool_t isSwitched = kFALSE;
1077       if (pTrack->GetSign() < 0) { // switch
1078
1079         isSwitched = kTRUE;
1080         AliESDtrack* helpTrack = nTrack;
1081         nTrack = pTrack;
1082         pTrack = helpTrack;
1083       }
1084
1085       Float_t alpha = esdV0->AlphaV0();
1086       Float_t ptarm = esdV0->PtArmV0();
1087       // Double_t pVtxPos= v0->PrimaryVtxPosition();
1088
1089       Double_t  lV0Position[3];
1090       esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1091
1092       Double_t lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
1093       Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1094                                             TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
1095                                             TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
1096       AliKFVertex primaryVtxKF( *myPrimaryVertex );
1097       AliKFParticle::SetField(ESDevent->GetMagneticField());
1098
1099       // Also implement switch here!!!!!!
1100       AliKFParticle* negEKF  = 0; // e-
1101       AliKFParticle* posEKF  = 0; // e+
1102       AliKFParticle* negPiKF = 0; // pi -
1103       AliKFParticle* posPiKF = 0; // pi +
1104       AliKFParticle* posPKF  = 0; // p
1105       AliKFParticle* negAPKF = 0; // p-bar
1106
1107       if(!isSwitched) {
1108         negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1109         posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1110         negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1111         posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1112         posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1113         negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1114       } else { // switch + and -
1115         negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1116         posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1117         negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1118         posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1119         posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1120         negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1121       }
1122
1123       AliKFParticle v0GKF;  // Gamma e.g. from pi0
1124       v0GKF+=(*negEKF);
1125       v0GKF+=(*posEKF);
1126       v0GKF.SetProductionVertex(primaryVtxKF);
1127
1128       AliKFParticle v0K0sKF; // K0 short
1129       v0K0sKF+=(*negPiKF);
1130       v0K0sKF+=(*posPiKF);
1131       v0K0sKF.SetProductionVertex(primaryVtxKF);
1132
1133       AliKFParticle v0LambdaKF; // Lambda
1134       v0LambdaKF+=(*negPiKF);
1135       v0LambdaKF+=(*posPKF);
1136       v0LambdaKF.SetProductionVertex(primaryVtxKF);
1137
1138       AliKFParticle v0AntiLambdaKF; // Lambda-bar
1139       v0AntiLambdaKF+=(*posPiKF);
1140       v0AntiLambdaKF+=(*negAPKF);
1141       v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1142
1143       Double_t deltaInvMassG     = v0GKF.GetMass();
1144       Double_t deltaInvMassK0s   = v0K0sKF.GetMass()-0.498;
1145       Double_t deltaInvMassL     = v0LambdaKF.GetMass()-1.116;
1146       Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
1147
1148       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1149          TMath::Abs(deltaInvMassL) > fMassCut &&
1150          TMath::Abs(deltaInvMassAntiL) > fMassCut)
1151         continue;
1152
1153       // Extract track information
1154
1155       Short_t pcharge  = pTrack->Charge();
1156       Float_t ppt      = pTrack->Pt();
1157       Float_t pp       = pTrack->P();
1158       Float_t peta     = pTrack->Eta();
1159       Float_t pphi     = pTrack->Phi();
1160       Short_t pncl     = pTrack->GetTPCsignalN();
1161       Short_t pneff    = Short_t(pTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1162       Float_t pdedx    = pTrack->GetTPCsignal();
1163
1164       Float_t ptpcchi  = 0;
1165       if(pTrack->GetTPCNcls() > 0)
1166         ptpcchi = pTrack->GetTPCchi2()/Float_t(pTrack->GetTPCNcls());
1167
1168       Short_t ncharge  = nTrack->Charge();
1169       Float_t npt      = nTrack->Pt();
1170       Float_t np       = nTrack->P();
1171       Float_t neta     = nTrack->Eta();
1172       Float_t nphi     = nTrack->Phi();
1173       Short_t nncl     = nTrack->GetTPCsignalN();
1174       Short_t nneff    = Short_t(nTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1175       Float_t ndedx    = nTrack->GetTPCsignal();
1176       Float_t ntpcchi  = 0;
1177       if(nTrack->GetTPCNcls() > 0)
1178         ntpcchi = nTrack->GetTPCchi2()/Float_t(nTrack->GetTPCNcls());
1179
1180       Float_t b[2];
1181       Float_t bCov[3];
1182       pTrack->GetImpactParameters(b,bCov);
1183       Float_t pdcaxy   = b[0];
1184       Float_t pdcaz    = b[1];
1185       nTrack->GetImpactParameters(b,bCov);
1186       Float_t ndcaxy   = b[0];
1187       Float_t ndcaz    = b[1];
1188
1189       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
1190       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
1191       Float_t p_ptMC        = 0;
1192       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
1193       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track
1194       Int_t   p_pdgMother   = 0;
1195       Float_t n_ptMC        = 0;
1196       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
1197       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track
1198       Int_t   n_pdgMother   = 0;
1199       if(fAnalysisMC) {
1200
1201         Int_t p_mother_label = 0;
1202         Int_t p_mother_steps = 0;
1203         Int_t n_mother_label = 0;
1204         Int_t n_mother_steps = 0;
1205
1206         // positive track
1207         const Int_t p_label = TMath::Abs(pTrack->GetLabel());
1208         TParticle* p_mcTrack = fMCStack->Particle(p_label);
1209         if (p_mcTrack){
1210
1211           if(fMCStack->IsPhysicalPrimary(p_label))
1212             p_primaryFlag = 1;
1213
1214           Int_t p_pdgCode = p_mcTrack->GetPdgCode();
1215           p_pidCode = GetPidCode(p_pdgCode);
1216
1217           p_ptMC      = p_mcTrack->Pt();
1218
1219           p_mother_label = FindPrimaryMotherLabel(fMCStack, p_label,
1220                                                   p_mother_steps);
1221           if(p_mother_label>0) {
1222             TParticle* p_mother = fMCStack->Particle(p_mother_label);
1223             p_pdgMother = p_mother->GetPdgCode();
1224           }
1225         }
1226
1227         // negative track
1228         const Int_t n_label = TMath::Abs(nTrack->GetLabel());
1229         TParticle* n_mcTrack = fMCStack->Particle(n_label);
1230         if (n_mcTrack){
1231
1232           if(fMCStack->IsPhysicalPrimary(n_label))
1233             n_primaryFlag = 1;
1234
1235           Int_t n_pdgCode = n_mcTrack->GetPdgCode();
1236           n_pidCode = GetPidCode(n_pdgCode);
1237
1238           n_ptMC      = n_mcTrack->Pt();
1239
1240           n_mother_label = FindPrimaryMotherLabel(fMCStack, n_label,
1241                                                   n_mother_steps);
1242           if(n_mother_label>0) {
1243             TParticle* n_mother = fMCStack->Particle(n_mother_label);
1244             n_pdgMother = n_mother->GetPdgCode();
1245           }
1246         }
1247
1248         // Check if V0 is primary = first and the same mother of both partciles
1249         if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
1250           pdgV0 = p_pdgMother;
1251           if(p_mother_steps == 1 && n_mother_steps == 1) {
1252             primaryV0 = 1;
1253           }
1254         }
1255       }
1256
1257
1258
1259       if(fTreeOption) {
1260
1261
1262
1263         DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
1264         nadded++;
1265
1266         // v0 data
1267         v0data->p       = esdV0->P();
1268         v0data->pt      = esdV0->Pt();
1269         v0data->eta     = esdV0->Eta();
1270         v0data->phi     = esdV0->Phi();
1271         v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
1272         v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
1273         v0data->dmassG  = deltaInvMassG;
1274         v0data->dmassK0 = deltaInvMassK0s;
1275         v0data->dmassL  = deltaInvMassL;
1276         v0data->dmassAL = deltaInvMassAntiL;
1277         v0data->alpha   = alpha;
1278         v0data->ptarm   = ptarm;
1279         v0data->decayr  = lV0Radius;
1280         v0data->decayl  = lV0DecayLength;
1281
1282         // New parameters
1283         v0data->status  = esdV0->GetOnFlyStatus();
1284         v0data->chi2    = esdV0->GetChi2V0();
1285         v0data->cospt   = esdV0->GetV0CosineOfPointingAngle();
1286         // cospt: as I understand this means that the pointing to the vertex
1287         // is fine so I remove the dcaxy and dcaz for the V= class
1288         v0data->dcadaughters = esdV0->GetDcaV0Daughters();
1289         v0data->primary = primaryV0;
1290         v0data->pdg     = pdgV0;
1291
1292         // positive track
1293         v0data->ptrack.p       = pp;
1294         v0data->ptrack.pt      = ppt;
1295         //        v0data->ptrack.ptcon   = ppt_con;
1296         //        v0data->ptrack.tpcchi  = ptpcchi;
1297         v0data->ptrack.eta     = peta;
1298         v0data->ptrack.phi     = pphi;
1299         v0data->ptrack.q       = pcharge;
1300         v0data->ptrack.ncl     = pncl;
1301         v0data->ptrack.neff    = pneff;
1302         v0data->ptrack.dedx    = pdedx;
1303         v0data->ptrack.dcaxy   = pdcaxy;
1304         v0data->ptrack.dcaz    = pdcaz;
1305         v0data->ptrack.pid     = p_pidCode;
1306         v0data->ptrack.primary = p_primaryFlag;
1307         v0data->ptrack.pttrue  = p_ptMC;
1308         v0data->ptrack.mother  = p_pdgMother;
1309         v0data->ptrack.filter  = filterFlag_p;
1310         v0data->ptrack.filterset1 = filterCut_Set1_p;
1311         v0data->ptrack.filterset2 = filterCut_Set2_p;
1312         v0data->ptrack.filterset3 = filterCut_Set3_p;
1313
1314         // negative track
1315         v0data->ntrack.p       = np;
1316         v0data->ntrack.pt      = npt;
1317         //        v0data->ntrack.ptcon   = npt_con;
1318         //        v0data->ntrack.tpcchi  = ntpcchi;
1319         v0data->ntrack.eta     = neta;
1320         v0data->ntrack.phi     = nphi;
1321         v0data->ntrack.q       = ncharge;
1322         v0data->ntrack.ncl     = nncl;
1323         v0data->ntrack.neff    = nneff;
1324         v0data->ntrack.dedx    = ndedx;
1325         v0data->ntrack.dcaxy   = ndcaxy;
1326         v0data->ntrack.dcaz    = ndcaz;
1327         v0data->ntrack.pid     = n_pidCode;
1328         v0data->ntrack.primary = n_primaryFlag;
1329         v0data->ntrack.pttrue  = n_ptMC;
1330         v0data->ntrack.mother  = n_pdgMother;
1331         v0data->ntrack.filter  = filterFlag_n;
1332         v0data->ntrack.filterset1 = filterCut_Set1_n;
1333         v0data->ntrack.filterset2 = filterCut_Set2_n;
1334         v0data->ntrack.filterset3 = filterCut_Set3_n;
1335
1336
1337       }
1338
1339       // clean up loop over v0
1340
1341       delete negPiKF;
1342       delete posPiKF;
1343       delete posPKF;
1344       delete negAPKF;
1345
1346
1347
1348     }
1349
1350     // clean up event
1351     //delete myPrimaryVertex;
1352   }
1353   else if( analysisMode==kTPCTrk ){
1354     cout<<"&&&&&&&&&&&&&&&&&&&&&                Hello world"<<endl;
1355     const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
1356     if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
1357
1358
1359     // ################################
1360     // #### BEGINNING OF V0 CODE ######
1361     // ################################
1362
1363
1364     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1365
1366       // This is the begining of the V0 loop
1367       AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
1368       if (!esdV0) continue;
1369
1370       // AliESDTrack (V0 Daughters)
1371       UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
1372       UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
1373
1374       AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
1375       AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
1376       if (!pTrack || !nTrack) {
1377         Printf("ERROR: Could not retreive one of the daughter track");
1378         continue;
1379       }
1380
1381       // Remove like-sign
1382       if (pTrack->GetSign() == nTrack->GetSign()) {
1383         //cout<< "like sign, continue"<< endl;
1384         continue;
1385       }
1386
1387       AliESDtrack *pTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),pTrack->GetID());
1388       AliESDtrack *nTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),nTrack->GetID());
1389
1390       if (!pTrackTPC || !nTrackTPC) {
1391         Printf("ERROR: Could not retreive one of the daughter TPC track");
1392         continue;
1393       }
1394
1395       //filter for positive track
1396       UInt_t selectDebug_p = 0;
1397       UShort_t filterFlag_p = 0;
1398       selectDebug_p = fTrackFilterTPC->IsSelected(pTrackTPC);
1399       //if(selectDebug_p==0) continue;
1400
1401       //filter for negative track
1402       UInt_t selectDebug_n = 0;
1403       UShort_t filterFlag_n = 0;
1404       selectDebug_n = fTrackFilterTPC->IsSelected(nTrackTPC);
1405       if(selectDebug_n==0 || selectDebug_p==0) continue;
1406
1407       if(selectDebug_p)
1408         filterFlag_p += 1;
1409
1410       if(pTrackTPC->Pt()>0.){
1411         // only constrain tracks above threshold
1412         AliExternalTrackParam exParamp;
1413         // take the B-field from the ESD, no 3D fieldMap available at this point
1414         Bool_t relate_p = false;
1415         relate_p = pTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1416                                                 kVeryBig,&exParamp);
1417         if(!relate_p){
1418           delete pTrackTPC;
1419           continue;
1420         }
1421         pTrackTPC->Set(exParamp.GetX(),exParamp.GetAlpha(),exParamp.GetParameter(),
1422                        exParamp.GetCovariance());
1423       }
1424       else continue;
1425
1426
1427       //filter for negative track
1428
1429       if(selectDebug_n)
1430         filterFlag_n += 1;
1431
1432
1433       if(nTrackTPC->Pt()>0.){
1434         // only constrain tracks above threshold
1435         AliExternalTrackParam exParamn;
1436         // take the B-field from the ESD, no 3D fieldMap available at this point
1437         Bool_t relate_n = false;
1438         relate_n = nTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1439                                                 kVeryBig,&exParamn);
1440         if(!relate_n){
1441           delete nTrackTPC;
1442           continue;
1443         }
1444         nTrackTPC->Set(exParamn.GetX(),exParamn.GetAlpha(),exParamn.GetParameter(),
1445                        exParamn.GetCovariance());
1446       }
1447       else continue;
1448
1449
1450       // Eta cut on decay products
1451       if(TMath::Abs(pTrackTPC->Eta()) > fEtaCut || TMath::Abs(nTrackTPC->Eta()) > fEtaCut)
1452         continue;
1453
1454       // Pt cut on decay products
1455       if (esdV0->Pt() < fMinPt)
1456         //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1457         continue;
1458
1459
1460       // Check if switch does anything!
1461       Bool_t isSwitched = kFALSE;
1462       if (pTrackTPC->GetSign() < 0) { // switch
1463
1464         isSwitched = kTRUE;
1465         AliESDtrack* helpTrack = nTrack;
1466         nTrackTPC = pTrackTPC;
1467         pTrackTPC = helpTrack;
1468       }
1469
1470
1471       // Extract track information
1472
1473       Short_t pcharge  = pTrackTPC->Charge();
1474       Float_t ppt      = pTrackTPC->Pt();
1475       Float_t pp       = pTrackTPC->P();
1476       Float_t peta     = pTrackTPC->Eta();
1477       Float_t pphi     = pTrackTPC->Phi();
1478       Short_t pncl     = pTrackTPC->GetTPCsignalN();
1479       Short_t pneff    = Short_t(pTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1480       Float_t pdedx    = pTrackTPC->GetTPCsignal();
1481
1482       Float_t ptpcchi  = 0;
1483       if(pTrackTPC->GetTPCNcls() > 0)
1484         ptpcchi = pTrackTPC->GetTPCchi2()/Float_t(pTrackTPC->GetTPCNcls());
1485
1486       Short_t ncharge  = nTrackTPC->Charge();
1487       Float_t npt      = nTrackTPC->Pt();
1488       Float_t np       = nTrackTPC->P();
1489       Float_t neta     = nTrackTPC->Eta();
1490       Float_t nphi     = nTrackTPC->Phi();
1491       Short_t nncl     = nTrackTPC->GetTPCsignalN();
1492       Short_t nneff    = Short_t(nTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1493       Float_t ndedx    = nTrackTPC->GetTPCsignal();
1494       Float_t ntpcchi  = 0;
1495       if(nTrackTPC->GetTPCNcls() > 0)
1496         ntpcchi = nTrackTPC->GetTPCchi2()/Float_t(nTrackTPC->GetTPCNcls());
1497
1498       Float_t bp[2]={0,0};
1499       Float_t bCovp[3]={0,0,0};
1500       pTrackTPC->GetImpactParameters(bp,bCovp);
1501       Float_t pdcaxy   = bp[0];
1502       Float_t pdcaz    = bp[1];
1503
1504       Float_t bn[2]={0,0};
1505       Float_t bCovn[3]={0,0,0};
1506       nTrackTPC->GetImpactParameters(bn,bCovn);
1507       Float_t ndcaxy   = bn[0];
1508       Float_t ndcaz    = bn[1];
1509
1510
1511       Float_t alpha = esdV0->AlphaV0();
1512       Float_t ptarm = esdV0->PtArmV0();
1513       // Double_t pVtxPos= v0->PrimaryVtxPosition();
1514
1515       Double_t  lV0Position[3];
1516       esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1517
1518       Double_t lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
1519       Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1520                                             TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
1521                                             TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
1522       AliKFVertex primaryVtxKF( *myPrimaryVertex );
1523       AliKFParticle::SetField(ESDevent->GetMagneticField());
1524
1525       // Also implement switch here!!!!!!
1526       AliKFParticle* negEKF  = 0; // e-
1527       AliKFParticle* posEKF  = 0; // e+
1528       AliKFParticle* negPiKF = 0; // pi -
1529       AliKFParticle* posPiKF = 0; // pi +
1530       AliKFParticle* posPKF  = 0; // p
1531       AliKFParticle* negAPKF = 0; // p-bar
1532
1533
1534       if(!isSwitched) {
1535         /*
1536                 negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1537                 posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1538                 negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1539                 posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1540                 posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1541                 negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1542         */
1543
1544         negEKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) , 11);
1545         posEKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) ,-11);
1546         negPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-211);
1547         posPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 211);
1548         posPKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 2212);
1549         negAPKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-2212);
1550
1551
1552
1553       } else { // switch + and -
1554         /*
1555         negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1556         posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1557         negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1558         posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1559         posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1560         negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1561         */
1562
1563         negEKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)), 11);
1564         posEKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)),-11);
1565         negPiKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)),-211);
1566         posPiKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)), 211);
1567         posPKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)), 2212);
1568         negAPKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)),-2212);
1569
1570
1571       }
1572
1573
1574
1575
1576
1577       AliKFParticle v0GKF;  // Gamma e.g. from pi0
1578       v0GKF+=(*negEKF);
1579       v0GKF+=(*posEKF);
1580       v0GKF.SetProductionVertex(primaryVtxKF);
1581
1582       AliKFParticle v0K0sKF; // K0 short
1583       v0K0sKF+=(*negPiKF);
1584       v0K0sKF+=(*posPiKF);
1585       v0K0sKF.SetProductionVertex(primaryVtxKF);
1586
1587       AliKFParticle v0LambdaKF; // Lambda
1588       v0LambdaKF+=(*negPiKF);
1589       v0LambdaKF+=(*posPKF);
1590       v0LambdaKF.SetProductionVertex(primaryVtxKF);
1591
1592       AliKFParticle v0AntiLambdaKF; // Lambda-bar
1593       v0AntiLambdaKF+=(*posPiKF);
1594       v0AntiLambdaKF+=(*negAPKF);
1595       v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1596
1597       Double_t deltaInvMassG     = v0GKF.GetMass();
1598       Double_t deltaInvMassK0s   = v0K0sKF.GetMass()-0.498;
1599       Double_t deltaInvMassL     = v0LambdaKF.GetMass()-1.116;
1600       Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
1601
1602       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1603          TMath::Abs(deltaInvMassL) > fMassCut &&
1604          TMath::Abs(deltaInvMassAntiL) > fMassCut)
1605         continue;
1606
1607       // TODO: Whe should these be different? Different mass hypothesis = energy loss
1608       // This is not important for us as we focus on the decay products!
1609       // Double_t ptK0s        = v0K0sKF.GetPt();
1610       // Double_t ptL          = v0LambdaKF.GetPt();
1611       // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();
1612
1613       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
1614       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
1615       Float_t p_ptMC        = 0;
1616       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
1617       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track
1618       Int_t   p_pdgMother   = 0;
1619       Float_t n_ptMC        = 0;
1620       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
1621       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track
1622       Int_t   n_pdgMother   = 0;
1623       if(fAnalysisMC) {
1624
1625         Int_t p_mother_label = 0;
1626         Int_t p_mother_steps = 0;
1627         Int_t n_mother_label = 0;
1628         Int_t n_mother_steps = 0;
1629
1630         // positive track
1631         const Int_t p_label = TMath::Abs(pTrackTPC->GetLabel());
1632         TParticle* p_mcTrack = fMCStack->Particle(p_label);
1633         if (p_mcTrack){
1634
1635           if(fMCStack->IsPhysicalPrimary(p_label))
1636             p_primaryFlag = 1;
1637
1638           Int_t p_pdgCode = p_mcTrack->GetPdgCode();
1639           p_pidCode = GetPidCode(p_pdgCode);
1640
1641           p_ptMC      = p_mcTrack->Pt();
1642
1643           p_mother_label = FindPrimaryMotherLabel(fMCStack, p_label,
1644                                                   p_mother_steps);
1645           if(p_mother_label>0) {
1646             TParticle* p_mother = fMCStack->Particle(p_mother_label);
1647             p_pdgMother = p_mother->GetPdgCode();
1648           }
1649         }
1650
1651         // negative track
1652         const Int_t n_label = TMath::Abs(nTrackTPC->GetLabel());
1653         TParticle* n_mcTrack = fMCStack->Particle(n_label);
1654         if (n_mcTrack){
1655
1656           if(fMCStack->IsPhysicalPrimary(n_label))
1657             n_primaryFlag = 1;
1658
1659           Int_t n_pdgCode = n_mcTrack->GetPdgCode();
1660           n_pidCode = GetPidCode(n_pdgCode);
1661
1662           n_ptMC      = n_mcTrack->Pt();
1663
1664           n_mother_label = FindPrimaryMotherLabel(fMCStack, n_label,
1665                                                   n_mother_steps);
1666           if(n_mother_label>0) {
1667             TParticle* n_mother = fMCStack->Particle(n_mother_label);
1668             n_pdgMother = n_mother->GetPdgCode();
1669           }
1670         }
1671
1672         // Check if V0 is primary = first and the same mother of both partciles
1673         if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
1674           pdgV0 = p_pdgMother;
1675           if(p_mother_steps == 1 && n_mother_steps == 1) {
1676             primaryV0 = 1;
1677           }
1678         }
1679       }
1680
1681
1682       if(fTreeOption) {
1683
1684         DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
1685         nadded++;
1686
1687         // v0 data
1688         v0datatpc->p       = esdV0->P();
1689         v0datatpc->pt      = esdV0->Pt();
1690         v0datatpc->eta     = esdV0->Eta();
1691         v0datatpc->phi     = esdV0->Phi();
1692         v0datatpc->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
1693         v0datatpc->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
1694         v0datatpc->dmassG  = deltaInvMassG;
1695         v0datatpc->dmassK0 = deltaInvMassK0s;
1696         v0datatpc->dmassL  = deltaInvMassL;
1697         v0datatpc->dmassAL = deltaInvMassAntiL;
1698         v0datatpc->alpha   = alpha;
1699         v0datatpc->ptarm   = ptarm;
1700         v0datatpc->decayr  = lV0Radius;
1701         v0datatpc->decayl  = lV0DecayLength;
1702
1703         // New parameters
1704         v0datatpc->status  = esdV0->GetOnFlyStatus();
1705         v0datatpc->chi2    = esdV0->GetChi2V0();
1706         v0datatpc->cospt   = esdV0->GetV0CosineOfPointingAngle();
1707         // cospt: as I understand this means that the pointing to the vertex
1708         // is fine so I remove the dcaxy and dcaz for the V= class
1709         v0datatpc->dcadaughters = esdV0->GetDcaV0Daughters();
1710         v0datatpc->primary = primaryV0;
1711         v0datatpc->pdg     = pdgV0;
1712
1713         // positive track
1714         v0datatpc->ptrack.p       = pp;
1715         v0datatpc->ptrack.pt      = ppt;
1716         //        v0data->ptrack.ptcon   = ppt_con;
1717         //        v0data->ptrack.tpcchi  = ptpcchi;
1718         v0datatpc->ptrack.eta     = peta;
1719         v0datatpc->ptrack.phi     = pphi;
1720         v0datatpc->ptrack.q       = pcharge;
1721         v0datatpc->ptrack.ncl     = pncl;
1722         v0datatpc->ptrack.neff    = pneff;
1723         v0datatpc->ptrack.dedx    = pdedx;
1724         v0datatpc->ptrack.dcaxy   = pdcaxy;
1725         v0datatpc->ptrack.dcaz    = pdcaz;
1726         v0datatpc->ptrack.pid     = p_pidCode;
1727         v0datatpc->ptrack.primary = p_primaryFlag;
1728         v0datatpc->ptrack.pttrue  = p_ptMC;
1729         v0datatpc->ptrack.mother  = p_pdgMother;
1730         v0datatpc->ptrack.filter  = filterFlag_p;
1731         v0datatpc->ptrack.filterset1 = 0;
1732         v0datatpc->ptrack.filterset2 = 0;
1733         v0datatpc->ptrack.filterset3 = 0;
1734
1735         // negative track
1736         v0datatpc->ntrack.p       = np;
1737         v0datatpc->ntrack.pt      = npt;
1738         //        v0data->ntrack.ptcon   = npt_con;
1739         //        v0data->ntrack.tpcchi  = ntpcchi;
1740         v0datatpc->ntrack.eta     = neta;
1741         v0datatpc->ntrack.phi     = nphi;
1742         v0datatpc->ntrack.q       = ncharge;
1743         v0datatpc->ntrack.ncl     = nncl;
1744         v0datatpc->ntrack.neff    = nneff;
1745         v0datatpc->ntrack.dedx    = ndedx;
1746         v0datatpc->ntrack.dcaxy   = ndcaxy;
1747         v0datatpc->ntrack.dcaz    = ndcaz;
1748         v0datatpc->ntrack.pid     = n_pidCode;
1749         v0datatpc->ntrack.primary = n_primaryFlag;
1750         v0datatpc->ntrack.pttrue  = n_ptMC;
1751         v0datatpc->ntrack.mother  = n_pdgMother;
1752         v0datatpc->ntrack.filter  = filterFlag_n;
1753         v0datatpc->ntrack.filterset1 = 0;
1754         v0datatpc->ntrack.filterset2 = 0;
1755         v0datatpc->ntrack.filterset3 = 0;
1756       }
1757
1758       // clean up loop over v0
1759
1760       delete negPiKF;
1761       delete posPiKF;
1762       delete posPKF;
1763       delete negAPKF;
1764
1765
1766
1767     }
1768
1769
1770
1771
1772
1773   }
1774   delete myPrimaryVertex;
1775
1776   if(fTreeOption) {
1777
1778     if( analysisMode==kGlobalTrk ){
1779
1780
1781       fEvent->trackmult = trackmult;
1782       fEvent->n         = nadded;
1783
1784
1785     }
1786
1787
1788   }
1789
1790
1791
1792 }
1793 //_______________________________________________________________________________________________________________
1794 void AliAnalysisTaskHighPtDeDxV0::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
1795   Int_t nv0s = AODevent->GetNumberOfV0s();
1796   if(nv0s<1)return;
1797   Int_t     trackmult = 0; // no pt cuts
1798   Int_t     nadded    = 0;
1799
1800   AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
1801   if (!myBestPrimaryVertex) return;
1802
1803
1804   if( analysisMode == kGlobalTrk ){
1805     if(fV0ArrayGlobalPar)
1806       fV0ArrayGlobalPar->Clear();
1807   } else if( analysisMode == kTPCTrk ){
1808     if(fV0ArrayTPCPar)
1809       fV0ArrayTPCPar->Clear();
1810   }
1811
1812   if( analysisMode==kGlobalTrk ){
1813
1814
1815     // ################################
1816     // #### BEGINNING OF V0 CODE ######
1817     // ################################
1818     // This is the begining of the V0 loop
1819     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1820       AliAODv0 *aodV0 = AODevent->GetV0(iV0);
1821       if (!aodV0) continue;
1822
1823       // common part
1824
1825       // AliAODTrack (V0 Daughters)
1826       AliAODVertex* vertex = aodV0->GetSecondaryVtx();
1827       if (!vertex) {
1828         Printf("ERROR: Could not retrieve vertex");
1829         continue;
1830       }
1831
1832       AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
1833       AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
1834       if (!pTrack || !nTrack) {
1835         Printf("ERROR: Could not retrieve one of the daughter track");
1836         continue;
1837       }
1838
1839       // Remove like-sign
1840       if (pTrack->Charge() == nTrack->Charge()) {
1841         //cout<< "like sign, continue"<< endl;
1842         continue;
1843       }
1844
1845       // Make sure charge ordering is ok
1846       if (pTrack->Charge() < 0) {
1847         AliAODTrack* helpTrack = pTrack;
1848         pTrack = nTrack;
1849         nTrack = helpTrack;
1850       }
1851
1852       // Eta cut on decay products
1853       if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1854         continue;
1855
1856       // Pt cut on decay products
1857       if (aodV0->Pt() < fMinPt)
1858         //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1859         continue;
1860
1861       //check positive tracks
1862       UShort_t filterFlag_p = 0;
1863       Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1864       Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1865       Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1866
1867       if (fTrackFilterGolden) {
1868         // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1869
1870         if(pTrack->TestFilterBit(32)) {
1871           filterFlag_p +=1;
1872           filterCut_Set3_p = kTRUE;
1873         }
1874       }
1875
1876
1877       if (fTrackFilterTPC) {
1878         // TPC only cuts is bit 1 according to above macro
1879         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
1880         if(pTrack->TestFilterBit(1)){
1881           filterFlag_p +=2;
1882           filterCut_Set1_p = kTRUE;
1883
1884         }
1885       }
1886
1887       if(filterFlag_p==0)
1888         continue;
1889
1890       //check negative tracks
1891       UShort_t filterFlag_n = 0;
1892       Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1893       Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1894       Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1895
1896
1897       if (fTrackFilterGolden) {
1898
1899         // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1900         if(nTrack->TestFilterBit(32)) {
1901           filterFlag_n +=1;
1902           filterCut_Set3_n = kTRUE;
1903         }
1904       }
1905
1906
1907       if (fTrackFilterTPC) {
1908         // TPC only cuts is bit 1 according to above macro
1909         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
1910         if(nTrack->TestFilterBit(1)){
1911           filterFlag_n +=2;
1912           filterCut_Set1_n = kTRUE;
1913
1914         }
1915       }
1916
1917       if(filterFlag_n==0)
1918         continue;
1919
1920
1921       Float_t alpha = aodV0->AlphaV0();
1922       Float_t ptarm = aodV0->PtArmV0();
1923       // Double_t pVtxPos= v0->PrimaryVtxPosition();
1924
1925       Double_t lV0Radius      = aodV0->RadiusV0();
1926       Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
1927
1928       Double_t deltaInvMassG     = aodV0->InvMass2Prongs(0,1,11,11);
1929       Double_t deltaInvMassK0s   = aodV0->MassK0Short()-0.498;
1930       Double_t deltaInvMassL     = aodV0->MassLambda()-1.116;
1931       Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
1932
1933       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1934          TMath::Abs(deltaInvMassL) > fMassCut &&
1935          TMath::Abs(deltaInvMassAntiL) > fMassCut)
1936         continue;
1937
1938       // TODO: Why should these be different? Different mass hypothesis = energy loss
1939       // This is not important for us as we focus on the decay products!
1940       // Double_t ptK0s        = v0K0sKF.GetPt();
1941       // Double_t ptL          = v0LambdaKF.GetPt();
1942       // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();
1943
1944       // Extract track information
1945
1946       Double_t b[2], cov[3];
1947       if(!pTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1948         filterFlag_p += 32; // propagation failed!!!!!
1949
1950       Float_t pdcaxy   = b[0];
1951       Float_t pdcaz    = b[1];
1952       if(!nTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1953         filterFlag_n += 32; // propagation failed!!!!!
1954       Float_t ndcaxy   = b[0];
1955       Float_t ndcaz    = b[1];
1956
1957       Short_t pcharge  = pTrack->Charge();
1958       Float_t ppt      = pTrack->Pt();
1959       Float_t pp       = pTrack->P();
1960       Float_t peta     = pTrack->Eta();
1961       Float_t pphi     = pTrack->Phi();
1962       //        Float_t ptpcchi  = pTrack->Chi2perNDF();
1963
1964       AliAODPid* pPid = pTrack->GetDetPid();
1965       Short_t pncl     = -10;
1966       Short_t pneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1967       Float_t pdedx    = -10;
1968       Float_t pbeta = -99;
1969       if(pPid) {
1970         pncl     = pPid->GetTPCsignalN();
1971         pdedx    = pPid->GetTPCsignal();
1972         //TOF
1973         if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
1974           Double_t tof[5];
1975           pPid->GetIntegratedTimes(tof);
1976           pbeta = tof[0]/pPid->GetTOFsignal();
1977         }
1978       }
1979
1980       Short_t ncharge  = nTrack->Charge();
1981       Float_t npt      = nTrack->Pt();
1982       Float_t np       = nTrack->P();
1983       Float_t neta     = nTrack->Eta();
1984       Float_t nphi     = nTrack->Phi();
1985       //        Float_t ntpcchi  = nTrack->Chi2perNDF();
1986
1987       AliAODPid* nPid = nTrack->GetDetPid();
1988       Short_t nncl     = -10;
1989       Short_t nneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1990       Float_t ndedx    = -10;
1991       Float_t nbeta = -99;
1992       if(pPid) {
1993         nncl     = nPid->GetTPCsignalN();
1994         ndedx    = nPid->GetTPCsignal();
1995         //TOF
1996         if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
1997           Double_t tof[5];
1998           nPid->GetIntegratedTimes(tof);
1999           nbeta = tof[0]/nPid->GetTOFsignal();
2000         }
2001       }
2002
2003       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
2004       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
2005       Float_t p_ptMC        = 0;
2006       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
2007       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track
2008       Int_t   p_pdgMother   = 0;
2009       Float_t n_ptMC        = 0;
2010       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
2011       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track
2012       Int_t   n_pdgMother   = 0;
2013       if(fAnalysisMC) {
2014
2015         AliAODMCParticle* p_mother = 0;
2016         Int_t p_mother_steps = 0;
2017         AliAODMCParticle* n_mother = 0;
2018         Int_t n_mother_steps = 0;
2019
2020         // positive track
2021         const Int_t p_label = TMath::Abs(pTrack->GetLabel());
2022
2023         AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
2024         if (p_mcTrack){
2025
2026           if(p_mcTrack->IsPhysicalPrimary())
2027             p_primaryFlag = 1;
2028
2029           Int_t p_pdgCode = p_mcTrack->GetPdgCode();
2030           p_pidCode = GetPidCode(p_pdgCode);
2031
2032           p_ptMC      = p_mcTrack->Pt();
2033
2034           p_mother = FindPrimaryMotherAOD(p_mcTrack, p_mother_steps);
2035           if(p_mother)
2036             p_pdgMother = p_mother->GetPdgCode();
2037         }
2038
2039         // negative track
2040         const Int_t n_label = TMath::Abs(pTrack->GetLabel());
2041
2042         AliAODMCParticle* n_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(n_label));
2043         if (n_mcTrack){
2044
2045           if(n_mcTrack->IsPhysicalPrimary())
2046             n_primaryFlag = 1;
2047
2048           Int_t n_pdgCode = n_mcTrack->GetPdgCode();
2049           n_pidCode = GetPidCode(n_pdgCode);
2050
2051           n_ptMC      = n_mcTrack->Pt();
2052
2053           n_mother = FindPrimaryMotherAOD(n_mcTrack, n_mother_steps);
2054           if(n_mother)
2055             n_pdgMother = n_mother->GetPdgCode();
2056         }
2057
2058         // Check if V0 is primary = first and the same mother of both partciles
2059         if(p_mother && n_mother && p_mother == n_mother) {
2060           pdgV0 = p_pdgMother;
2061           if(p_mother_steps == 1 && n_mother_steps == 1) {
2062             primaryV0 = 1;
2063           }
2064         }
2065       }
2066
2067       if(fTreeOption) {
2068
2069         DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
2070         nadded++;
2071
2072         // v0 data
2073         v0data->p       = aodV0->P();
2074         v0data->pt      = aodV0->Pt();
2075         v0data->eta     = aodV0->Eta();
2076         v0data->phi     = aodV0->Phi();
2077         v0data->pdca    = aodV0->DcaPosToPrimVertex();
2078         v0data->ndca    = aodV0->DcaNegToPrimVertex();
2079         v0data->dmassG  = deltaInvMassG;
2080         v0data->dmassK0 = deltaInvMassK0s;
2081         v0data->dmassL  = deltaInvMassL;
2082         v0data->dmassAL = deltaInvMassAntiL;
2083         v0data->alpha   = alpha;
2084         v0data->ptarm   = ptarm;
2085         v0data->decayr  = lV0Radius;
2086         v0data->decayl  = lV0DecayLength;
2087         // v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
2088         // v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
2089
2090         // New parameters
2091         v0data->status  = aodV0->GetOnFlyStatus();
2092         v0data->chi2    = aodV0->Chi2V0();
2093         v0data->cospt   = aodV0->CosPointingAngle(myBestPrimaryVertex);
2094         // cospt: as I understand this means that the pointing to the vertex
2095         // is fine so I remove the dcaxy and dcaz for the V= class
2096         v0data->dcav0   = aodV0->DcaV0ToPrimVertex();
2097         v0data->dcadaughters = aodV0->DcaV0Daughters();
2098         v0data->primary = primaryV0;
2099         v0data->pdg     = pdgV0;
2100
2101         // positive track
2102         v0data->ptrack.p       = pp;
2103         v0data->ptrack.pt      = ppt;
2104         //        v0data->ptrack.ptcon   = ppt_con;
2105         //        v0data->ptrack.tpcchi  = ptpcchi;
2106         v0data->ptrack.eta     = peta;
2107         v0data->ptrack.phi     = pphi;
2108         v0data->ptrack.q       = pcharge;
2109         v0data->ptrack.ncl     = pncl;
2110         v0data->ptrack.neff    = pneff;
2111         v0data->ptrack.dedx    = pdedx;
2112         v0data->ptrack.dcaxy   = pdcaxy;
2113         v0data->ptrack.dcaz    = pdcaz;
2114         v0data->ptrack.pid     = p_pidCode;
2115         v0data->ptrack.primary = p_primaryFlag;
2116         v0data->ptrack.pttrue  = p_ptMC;
2117         v0data->ptrack.mother  = p_pdgMother;
2118         v0data->ptrack.filter  = filterFlag_p;
2119         v0data->ptrack.filterset1 = filterCut_Set1_p;
2120         v0data->ptrack.filterset2 = filterCut_Set2_p;
2121         v0data->ptrack.filterset3 = filterCut_Set3_p;
2122
2123
2124         // negative track
2125         v0data->ntrack.p       = np;
2126         v0data->ntrack.pt      = npt;
2127         //        v0data->ntrack.ptcon   = npt_con;
2128         //        v0data->ntrack.tpcchi  = ntpcchi;
2129         v0data->ntrack.eta     = neta;
2130         v0data->ntrack.phi     = nphi;
2131         v0data->ntrack.q       = ncharge;
2132         v0data->ntrack.ncl     = nncl;
2133         v0data->ntrack.neff    = nneff;
2134         v0data->ntrack.dedx    = ndedx;
2135         v0data->ntrack.dcaxy   = ndcaxy;
2136         v0data->ntrack.dcaz    = ndcaz;
2137         v0data->ntrack.pid     = n_pidCode;
2138         v0data->ntrack.primary = n_primaryFlag;
2139         v0data->ntrack.pttrue  = n_ptMC;
2140         v0data->ntrack.mother  = n_pdgMother;
2141         v0data->ntrack.filter  = filterFlag_n;
2142         v0data->ntrack.filterset1 = filterCut_Set1_n;
2143         v0data->ntrack.filterset2 = filterCut_Set2_n;
2144         v0data->ntrack.filterset3 = filterCut_Set3_n;
2145
2146
2147       }
2148     }//end loop over v0's
2149
2150
2151   }else if( analysisMode==kTPCTrk ){
2152
2153     const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
2154     cout<<"&&&&&&&&&&&&&&&&&&&&&                Hello world  0"<<endl;
2155     if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
2156
2157     cout<<"&&&&&&&&&&&&&&&&&&&&&                Hello world"<<endl;
2158     // ################################
2159     // #### BEGINNING OF V0 CODE ######
2160     // ################################
2161     // This is the begining of the V0 loop
2162
2163
2164
2165      for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
2166        AliAODv0 *aodV0 = AODevent->GetV0(iV0);
2167        if (!aodV0) continue;
2168
2169        // common part
2170
2171        // AliAODTrack (V0 Daughters)
2172        AliAODVertex* vertex = aodV0->GetSecondaryVtx();
2173        if (!vertex) {
2174          Printf("ERROR: Could not retrieve vertex");
2175          continue;
2176        }
2177
2178        AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
2179        AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
2180        if (!pTrack || !nTrack) {
2181          Printf("ERROR: Could not retrieve one of the daughter track");
2182          continue;
2183        }
2184
2185
2186        // Remove like-sign
2187        if (pTrack->Charge() == nTrack->Charge()) {
2188          //cout<< "like sign, continue"<< endl;
2189          continue;
2190        }
2191
2192        // Make sure charge ordering is ok
2193        if (pTrack->Charge() < 0) {
2194          AliAODTrack* helpTrack = pTrack;
2195          pTrack = nTrack;
2196          nTrack = helpTrack;
2197        }
2198
2199        // Eta cut on decay products
2200        if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
2201          continue;
2202
2203        cout<<"Eta positive track:"<<pTrack->Eta()<<endl;
2204
2205
2206        // Pt cut on decay products
2207        if (aodV0->Pt() < fMinPt)
2208          //     if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
2209          continue;
2210
2211        //check positive tracks
2212        UShort_t filterFlag_p = 0;
2213        Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
2214        Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
2215        Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
2216
2217        // TPC only cuts is bit 1 according to above macro
2218        // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
2219        if(pTrack->TestFilterBit(128)) {
2220          cout<<"este track paso el corte bit 128"<<endl;
2221          filterFlag_p +=1;
2222        }
2223        cout<<"filterFlag_p="<<filterFlag_p<<endl;
2224
2225
2226
2227
2228        if(filterFlag_p==0)
2229          continue;
2230
2231
2232
2233        //check negative tracks
2234        UShort_t filterFlag_n = 0;
2235        Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
2236        Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
2237        Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
2238
2239        // TPC only cuts is bit 1 according to above macro
2240        // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
2241        if(nTrack->TestFilterBit(128)) {
2242          filterFlag_n +=1;
2243        }
2244
2245        if(filterFlag_n==0)
2246          continue;
2247
2248
2249        Float_t alpha = aodV0->AlphaV0();
2250        Float_t ptarm = aodV0->PtArmV0();
2251        // Double_t pVtxPos= v0->PrimaryVtxPosition();
2252
2253        Double_t lV0Radius      = aodV0->RadiusV0();
2254        Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
2255
2256        Double_t deltaInvMassG     = aodV0->InvMass2Prongs(0,1,11,11);
2257        Double_t deltaInvMassK0s   = aodV0->MassK0Short()-0.498;
2258        Double_t deltaInvMassL     = aodV0->MassLambda()-1.116;
2259        Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
2260
2261        if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
2262           TMath::Abs(deltaInvMassL) > fMassCut &&
2263           TMath::Abs(deltaInvMassAntiL) > fMassCut)
2264          continue;
2265
2266        // TODO: Why should these be different? Different mass hypothesis = energy loss
2267        // This is not important for us as we focus on the decay products!
2268        // Double_t ptK0s        = v0K0sKF.GetPt();
2269        // Double_t ptL          = v0LambdaKF.GetPt();
2270        // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();
2271
2272        // Extract track information
2273
2274        Double_t b[2], cov[3];
2275        if(!pTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov))
2276          filterFlag_p += 32; // propagation failed!!!!!
2277
2278        Float_t pdcaxy   = b[0];
2279        Float_t pdcaz    = b[1];
2280        if(!nTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov))
2281          filterFlag_n += 32; // propagation failed!!!!!
2282        Float_t ndcaxy   = b[0];
2283        Float_t ndcaz    = b[1];
2284
2285        Short_t pcharge  = pTrack->Charge();
2286        Float_t ppt      = pTrack->Pt();
2287        Float_t pp       = pTrack->P();
2288        Float_t peta     = pTrack->Eta();
2289        Float_t pphi     = pTrack->Phi();
2290        //       Float_t ptpcchi  = pTrack->Chi2perNDF();
2291
2292        AliAODPid* pPid = pTrack->GetDetPid();
2293        Short_t pncl     = -10;
2294        Short_t pneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2295        Float_t pdedx    = -10;
2296        Float_t pbeta = -99;
2297        if(pPid) {
2298          pncl     = pPid->GetTPCsignalN();
2299          pdedx    = pPid->GetTPCsignal();
2300          //TOF
2301          if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
2302            Double_t tof[5];
2303            pPid->GetIntegratedTimes(tof);
2304            pbeta = tof[0]/pPid->GetTOFsignal();
2305          }
2306        }
2307
2308        Short_t ncharge  = nTrack->Charge();
2309        Float_t npt      = nTrack->Pt();
2310        Float_t np       = nTrack->P();
2311        Float_t neta     = nTrack->Eta();
2312        Float_t nphi     = nTrack->Phi();
2313        //       Float_t ntpcchi  = nTrack->Chi2perNDF();
2314
2315        AliAODPid* nPid = nTrack->GetDetPid();
2316        Short_t nncl     = -10;
2317        Short_t nneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2318        Float_t ndedx    = -10;
2319        Float_t nbeta = -99;
2320        if(pPid) {
2321          nncl     = nPid->GetTPCsignalN();
2322          ndedx    = nPid->GetTPCsignal();
2323          //TOF
2324          if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
2325            Double_t tof[5];
2326            nPid->GetIntegratedTimes(tof);
2327            nbeta = tof[0]/nPid->GetTOFsignal();
2328          }
2329        }
2330
2331        Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
2332        Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
2333        Float_t p_ptMC        = 0;
2334        Short_t p_pidCode     = 0; // 0 = real data / no mc track!
2335        Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track
2336        Int_t   p_pdgMother   = 0;
2337        Float_t n_ptMC        = 0;
2338        Short_t n_pidCode     = 0; // 0 = real data / no mc track!
2339        Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track
2340        Int_t   n_pdgMother   = 0;
2341        if(fAnalysisMC) {
2342
2343          AliAODMCParticle* p_mother = 0;
2344          Int_t p_mother_steps = 0;
2345          AliAODMCParticle* n_mother = 0;
2346          Int_t n_mother_steps = 0;
2347
2348          // positive track
2349          const Int_t p_label = TMath::Abs(pTrack->GetLabel());
2350
2351          AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
2352          if (p_mcTrack){
2353
2354            if(p_mcTrack->IsPhysicalPrimary())
2355              p_primaryFlag = 1;
2356
2357            Int_t p_pdgCode = p_mcTrack->GetPdgCode();
2358            p_pidCode = GetPidCode(p_pdgCode);
2359
2360            p_ptMC      = p_mcTrack->Pt();
2361
2362            p_mother = FindPrimaryMotherAOD(p_mcTrack, p_mother_steps);
2363            if(p_mother)
2364              p_pdgMother = p_mother->GetPdgCode();
2365          }
2366
2367          // negative track
2368          const Int_t n_label = TMath::Abs(pTrack->GetLabel());
2369
2370          AliAODMCParticle* n_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(n_label));
2371          if (n_mcTrack){
2372
2373            if(n_mcTrack->IsPhysicalPrimary())
2374              n_primaryFlag = 1;
2375
2376            Int_t n_pdgCode = n_mcTrack->GetPdgCode();
2377            n_pidCode = GetPidCode(n_pdgCode);
2378
2379            n_ptMC      = n_mcTrack->Pt();
2380
2381            n_mother = FindPrimaryMotherAOD(n_mcTrack, n_mother_steps);
2382            if(n_mother)
2383              n_pdgMother = n_mother->GetPdgCode();
2384          }
2385
2386          // Check if V0 is primary = first and the same mother of both partciles
2387          if(p_mother && n_mother && p_mother == n_mother) {
2388            pdgV0 = p_pdgMother;
2389            if(p_mother_steps == 1 && n_mother_steps == 1) {
2390              primaryV0 = 1;
2391            }
2392          }
2393        }
2394
2395        if(fTreeOption) {
2396
2397          DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
2398          nadded++;
2399
2400          // v0 data
2401          v0datatpc->p       = aodV0->P();
2402          v0datatpc->pt      = aodV0->Pt();
2403          v0datatpc->eta     = aodV0->Eta();
2404          v0datatpc->phi     = aodV0->Phi();
2405          v0datatpc->pdca    = aodV0->DcaPosToPrimVertex();
2406          v0datatpc->ndca    = aodV0->DcaNegToPrimVertex();
2407          v0datatpc->dmassG  = deltaInvMassG;
2408          v0datatpc->dmassK0 = deltaInvMassK0s;
2409          v0datatpc->dmassL  = deltaInvMassL;
2410          v0datatpc->dmassAL = deltaInvMassAntiL;
2411          v0datatpc->alpha   = alpha;
2412          v0datatpc->ptarm   = ptarm;
2413          v0datatpc->decayr  = lV0Radius;
2414          v0datatpc->decayl  = lV0DecayLength;
2415          // v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
2416          // v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
2417
2418          // New parameters
2419          v0datatpc->status  = aodV0->GetOnFlyStatus();
2420          v0datatpc->chi2    = aodV0->Chi2V0();
2421          v0datatpc->cospt   = aodV0->CosPointingAngle(myBestPrimaryVertex);
2422          // cospt: as I understand this means that the pointing to the vertex
2423          // is fine so I remove the dcaxy and dcaz for the V= class
2424          v0datatpc->dcav0   = aodV0->DcaV0ToPrimVertex();
2425          v0datatpc->dcadaughters = aodV0->DcaV0Daughters();
2426          v0datatpc->primary = primaryV0;
2427          v0datatpc->pdg     = pdgV0;
2428
2429          // positive track
2430          v0datatpc->ptrack.p       = pp;
2431          v0datatpc->ptrack.pt      = ppt;
2432          //       v0data->ptrack.ptcon   = ppt_con;
2433          //       v0data->ptrack.tpcchi  = ptpcchi;
2434          v0datatpc->ptrack.eta     = peta;
2435          v0datatpc->ptrack.phi     = pphi;
2436          v0datatpc->ptrack.q       = pcharge;
2437          v0datatpc->ptrack.ncl     = pncl;
2438          v0datatpc->ptrack.neff    = pneff;
2439          v0datatpc->ptrack.dedx    = pdedx;
2440          v0datatpc->ptrack.dcaxy   = pdcaxy;
2441          v0datatpc->ptrack.dcaz    = pdcaz;
2442          v0datatpc->ptrack.pid     = p_pidCode;
2443          v0datatpc->ptrack.primary = p_primaryFlag;
2444          v0datatpc->ptrack.pttrue  = p_ptMC;
2445          v0datatpc->ptrack.mother  = p_pdgMother;
2446          v0datatpc->ptrack.filter  = filterFlag_p;
2447          v0datatpc->ptrack.filterset1 = 0;
2448          v0datatpc->ptrack.filterset2 = 0;
2449          v0datatpc->ptrack.filterset3 = 0;
2450
2451
2452          // negative track
2453          v0datatpc->ntrack.p       = np;
2454          v0datatpc->ntrack.pt      = npt;
2455          //       v0data->ntrack.ptcon   = npt_con;
2456          //       v0data->ntrack.tpcchi  = ntpcchi;
2457          v0datatpc->ntrack.eta     = neta;
2458          v0datatpc->ntrack.phi     = nphi;
2459          v0datatpc->ntrack.q       = ncharge;
2460          v0datatpc->ntrack.ncl     = nncl;
2461          v0datatpc->ntrack.neff    = nneff;
2462          v0datatpc->ntrack.dedx    = ndedx;
2463          v0datatpc->ntrack.dcaxy   = ndcaxy;
2464          v0datatpc->ntrack.dcaz    = ndcaz;
2465          v0datatpc->ntrack.pid     = n_pidCode;
2466          v0datatpc->ntrack.primary = n_primaryFlag;
2467          v0datatpc->ntrack.pttrue  = n_ptMC;
2468          v0datatpc->ntrack.mother  = n_pdgMother;
2469          v0datatpc->ntrack.filter  = filterFlag_n;
2470          v0datatpc->ntrack.filterset1 = 0;
2471          v0datatpc->ntrack.filterset2 = 0;
2472          v0datatpc->ntrack.filterset3 = 0;
2473
2474
2475        }
2476      }//end v0's loop
2477
2478   }
2479
2480   if(fTreeOption) {
2481
2482     if( analysisMode==kGlobalTrk ){
2483
2484
2485       fEvent->trackmult = trackmult;
2486       fEvent->n         = nadded;
2487
2488
2489     }
2490
2491
2492   }
2493
2494
2495
2496 }
2497
2498
2499