PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtTrackDumpTask.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 #include "iostream"
  17
  18 #include <TPDGCode.h>
  19 #include <TDatabasePDG.h>
  20
  21 #include "TChain.h"
  22 #include "TTreeStream.h"
  23 #include "TTree.h"
  24 #include "TH1F.h"
  25 #include "TCanvas.h"
  26 #include "TList.h"
  27 #include "TObjArray.h"
  28 #include "TFile.h"
  29 #include "TMatrixD.h"
  30 #include "TRandom3.h"
  31
  32 #include "AliHeader.h"
  33 #include "AliGenEventHeader.h"
  34 #include "AliInputEventHandler.h"
  35 #include "AliStack.h"
  36 #include "AliTrackReference.h"
  37
  38 #include "AliPhysicsSelection.h"
  39 #include "AliAnalysisTask.h"
  40 #include "AliAnalysisManager.h"
  41 #include "AliESDEvent.h"
  42 #include "AliESDfriend.h"
  43 #include "AliMCEvent.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliESDVertex.h"
  46 #include "AliTracker.h"
  47 #include "AliGeomManager.h"
  48
  49 #include "AliCentrality.h"
  50 #include "AliESDVZERO.h"
  51 #include "AliMultiplicity.h"
  52
  53 #include "AliESDtrackCuts.h"
  54 #include "AliMCEventHandler.h"
  55 #include "AlidNdPt.h"
  56 #include "AlidNdPtEventCuts.h"
  57 #include "AlidNdPtAcceptanceCuts.h"
  58
  59 #include "AlidNdPtTrackDumpTask.h"
  60 #include "AliKFParticle.h"
  61 #include "AliESDv0.h"
  62
  63 using namespace std;
  64
  65 ClassImp(AlidNdPtTrackDumpTask)
  66
  67 //_____________________________________________________________________________
  68 AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name)
  69   : AliAnalysisTaskSE(name)
  70   , fESD(0)
  71   , fMC(0)
  72   , fESDfriend(0)
  73   , fOutput(0)
  74   , fPitList(0)
  75   , fUseMCInfo(kFALSE)
  76   , fUseESDfriends(kFALSE)
  77   , fdNdPtEventCuts(0)
  78   , fdNdPtAcceptanceCuts(0)
  79   , fdNdPtRecAcceptanceCuts(0)
  80   , fEsdTrackCuts(0)
  81   , fTrigger(AliTriggerAnalysis::kMB1)
  82   , fAnalysisMode(AlidNdPtHelper::kTPC)
  83   , fTreeSRedirector(0)
  84   , fCentralityEstimator(0)
  85   , fLowPtTrackDownscaligF(0)
  86   , fLowPtV0DownscaligF(0)
  87   , fProcessAll(kFALSE)
  88   , fProcessCosmics(kFALSE)
  89   , fTree1(0)
  90   , fTree2(0)
  91   , fTree3(0)
  92   , fTree4(0)
  93   , fTree5(0)
  94   , fTree6(0)
  95 {
  96   // Constructor
  97
  98   // Define input and output slots here
  99   DefineOutput(1, TTree::Class());
 100   DefineOutput(2, TTree::Class());
 101   DefineOutput(3, TTree::Class());
 102   DefineOutput(4, TTree::Class());
 103   DefineOutput(5, TTree::Class());
 104   DefineOutput(6, TTree::Class());
 105 }
 106
 107 //_____________________________________________________________________________
 108 AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()
 109 {
 110   if(fOutput) delete fOutput;  fOutput =0;
 111   if(fTreeSRedirector) delete fTreeSRedirector;  fTreeSRedirector =0;
 112
 113   if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL;
 114   if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;
 115   if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL;
 116   if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;
 117
 118
 119 }
 120
 121 //____________________________________________________________________________
 122 Bool_t AlidNdPtTrackDumpTask::Notify()
 123 {
 124   static Int_t count = 0;
 125   count++;
 126   TTree *chain = (TChain*)GetInputData(0);
 127   if(chain)
 128   {
 129     Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());
 130   }
 131
 132 return kTRUE;
 133 }
 134
 135 //_____________________________________________________________________________
 136 void AlidNdPtTrackDumpTask::UserCreateOutputObjects()
 137 {
 138   // Create histograms
 139   // Called once
 140   fOutput = new TList;
 141   fOutput->SetOwner();
 142
 143   //
 144   // create temporary file for output tree
 145   fTreeSRedirector = new TTreeSRedirector("jotwinow_Temp_Trees.root");
 146
 147   fTree1 = new TTree;
 148   fTree2 = new TTree;
 149   fTree3 = new TTree;
 150   fTree4 = new TTree;
 151   fTree5 = new TTree;
 152   fTree6 = new TTree;
 153
 154   fOutput->Add(fTree1);
 155   fOutput->Add(fTree2);
 156   fOutput->Add(fTree3);
 157   fOutput->Add(fTree4);
 158   fOutput->Add(fTree5);
 159   fOutput->Add(fTree6);
 160
 161   PostData(1, fTree1);
 162   PostData(2, fTree2);
 163   PostData(3, fTree3);
 164   PostData(4, fTree4);
 165   PostData(5, fTree5);
 166   PostData(6, fTree6);
 167
 168 }
 169
 170 //_____________________________________________________________________________
 171 void AlidNdPtTrackDumpTask::UserExec(Option_t *)
 172 {
 173   //
 174   // Called for each event
 175   //
 176
 177   // ESD event
 178   fESD = (AliESDEvent*) (InputEvent());
 179   if (!fESD) {
 180     Printf("ERROR: ESD event not available");
 181     return;
 182   }
 183
 184   // MC event
 185   if(fUseMCInfo) {
 186     fMC = MCEvent();
 187     if (!fMC) {
 188       Printf("ERROR: MC event not available");
 189       return;
 190     }
 191   }
 192
 193   if(fUseESDfriends) {
 194     fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 195       if(!fESDfriend) {
 196         Printf("ERROR: ESD friends not available");
 197     }
 198   }
 199
 200   //
 201   if(fProcessAll) {
 202     ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC
 203   }
 204   else {
 205     Process(fESD,fMC,fESDfriend);    // only global and TPC tracks
 206   }
 207
 208   //
 209   ProcessV0(fESD,fMC,fESDfriend);
 210   ProcessLaser(fESD,fMC,fESDfriend);
 211   ProcessdEdx(fESD,fMC,fESDfriend);
 212
 213   if (fProcessCosmics) { ProcessCosmics(fESD); }
 214   if(IsUseMCInfo()) { ProcessMCEff(fESD,fMC,fESDfriend); }
 215 }
 216
 217 //_____________________________________________________________________________
 218 void AlidNdPtTrackDumpTask::ProcessCosmics(AliESDEvent *const event)
 219 {
 220   //
 221   // Select real events with high-pT tracks
 222   //
 223   if(!event) {
 224     AliDebug(AliLog::kError, "event not available");
 225     return;
 226   }
 227
 228   //
 229   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 230   if (!inputHandler)
 231   {
 232     Printf("ERROR: Could not receive input handler");
 233     return;
 234   }
 235
 236   // get file name
 237   TTree *chain = (TChain*)GetInputData(0);
 238   if(!chain) {
 239     Printf("ERROR: Could not receive input chain");
 240     return;
 241   }
 242   TObjString fileName(chain->GetCurrentFile()->GetName());
 243
 244
 245     // check for cosmic pairs
 246     //
 247     // find cosmic pairs trigger by random trigger
 248     //
 249     //
 250     AliESDVertex *vertexSPD =  (AliESDVertex *)event->GetPrimaryVertexSPD();
 251     AliESDVertex *vertexTPC =  (AliESDVertex *)event->GetPrimaryVertexTPC();
 252     const Double_t kMinPt=0.8;
 253     const Double_t kMinPtMax=0.8;
 254     const Double_t kMinNcl=50;
 255     const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
 256     Int_t ntracks=event->GetNumberOfTracks();
 257     //  Float_t dcaTPC[2]={0,0};
 258     // Float_t covTPC[3]={0,0,0};
 259
 260     UInt_t specie = event->GetEventSpecie();  // skip laser events
 261     if (specie==AliRecoParam::kCalib) return;
 262
 263
 264
 265     for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
 266       AliESDtrack *track0 = event->GetTrack(itrack0);
 267       if (!track0) continue;
 268       if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
 269
 270       if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue;
 271       if (track0->Pt() < kMinPt) continue;
 272       if (track0->GetTPCncls() < kMinNcl) continue;
 273       if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
 274       if (track0->GetKinkIndex(0)>0) continue;
 275       const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
 276       //rm primaries
 277       //
 278       //track0->GetImpactParametersTPC(dcaTPC,covTPC);
 279       //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
 280       //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
 281       //    const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
 282       for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
 283         AliESDtrack *track1 = event->GetTrack(itrack1);
 284         if (!track1) continue;
 285         if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
 286         if (track1->GetKinkIndex(0)>0) continue;
 287         if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue;
 288         if (track1->Pt() < kMinPt) continue;
 289         if (track1->GetTPCncls()<kMinNcl) continue;
 290         if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
 291         if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
 292         //track1->GetImpactParametersTPC(dcaTPC,covTPC);
 293         //      if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
 294         //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
 295         //
 296         const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
 297         //
 298         Bool_t isPair=kTRUE;
 299         for (Int_t ipar=0; ipar<5; ipar++){
 300           if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
 301           if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
 302         }
 303         if (!isPair) continue;
 304         if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
 305         //delta with correct sign
 306         /*
 307         TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
 308         TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
 309         TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
 310         */
 311         if  (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y   opposite sign
 312         if  (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
 313         if  (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
 314         if (!isPair) continue;
 315         TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
 316         Int_t eventNumber = event->GetEventNumberInFile();
 317         //Bool_t hasFriend = kFALSE;
 318         //Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);
 319         //printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS);
 320         //      const AliExternalTrackParam * trackIn1 = track1->GetInnerParam();
 321         //
 322         //
 323         Int_t ntracksSPD = vertexSPD->GetNContributors();
 324         Int_t ntracksTPC = vertexTPC->GetNContributors();
 325         Int_t runNumber     = event->GetRunNumber();
 326         Int_t timeStamp    = event->GetTimeStamp();
 327         ULong64_t triggerMask = event->GetTriggerMask();
 328         Float_t magField    = event->GetMagneticField();
 329         TObjString triggerClass = event->GetFiredTriggerClasses().Data();
 330
 331        //
 332       // Dump to the tree
 333       // vertex
 334       // TPC-ITS tracks
 335       //
 336       if(!fTreeSRedirector) return;
 337           (*fTreeSRedirector)<<"CosmicPairs"<<
 338             "fileName.="<<&fileName<<         // file name
 339             "runNumber="<<runNumber<<              //  run number
 340             "evtTimeStamp="<<timeStamp<<            //  time stamp of event
 341             "evtNumberInFile="<<eventNumber<<          //  event number
 342             "trigger="<<triggerMask<<      //  trigger
 343             "triggerClass="<<&triggerClass<<      //  trigger
 344             "Bz="<<magField<<             //  magnetic field
 345             //
 346             "multSPD="<<ntracksSPD<<
 347             "multTPC="<<ntracksTPC<<
 348             "vertSPD.="<<vertexSPD<<         //primary vertex -SPD
 349             "vertTPC.="<<vertexTPC<<         //primary vertex -TPC
 350             "t0.="<<track0<<              //track0
 351             "t1.="<<track1<<              //track1
 352             "\n";
 353         }
 354       }
 355 }
 356
 357
 358 //_____________________________________________________________________________
 359 void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
 360 {
 361   //
 362   // Select real events with high-pT tracks
 363   //
 364   if(!esdEvent) {
 365     AliDebug(AliLog::kError, "esdEvent not available");
 366     return;
 367   }
 368
 369   // get selection cuts
 370   AlidNdPtEventCuts *evtCuts = GetEventCuts();
 371   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
 372   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
 373
 374   if(!evtCuts || !accCuts  || !esdTrackCuts) {
 375     Printf("ERROR cuts not available");
 376     return;
 377   }
 378
 379   // trigger selection
 380   Bool_t isEventTriggered = kTRUE;
 381   AliPhysicsSelection *physicsSelection = NULL;
 382   AliTriggerAnalysis* triggerAnalysis = NULL;
 383
 384   //
 385   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 386   if (!inputHandler)
 387   {
 388     Printf("ERROR: Could not receive input handler");
 389     return;
 390   }
 391
 392   // get file name
 393   TTree *chain = (TChain*)GetInputData(0);
 394   if(!chain) {
 395     Printf("ERROR: Could not receive input chain");
 396     return;
 397   }
 398   TObjString fileName(chain->GetCurrentFile()->GetName());
 399
 400   // trigger
 401   if(evtCuts->IsTriggerRequired())
 402   {
 403     // always MB
 404     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
 405
 406     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 407     if(!physicsSelection) return;
 408     //SetPhysicsTriggerSelection(physicsSelection);
 409
 410     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
 411       // set trigger (V0AND)
 412       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 413       if(!triggerAnalysis) return;
 414       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
 415     }
 416   }
 417
 418   // centrality determination
 419   Float_t centralityF = -1;
 420   AliCentrality *esdCentrality = esdEvent->GetCentrality();
 421   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
 422
 423   // use MC information
 424   AliHeader* header = 0;
 425   AliGenEventHeader* genHeader = 0;
 426   AliStack* stack = 0;
 427   TArrayF vtxMC(3);
 428
 429   Int_t multMCTrueTracks = 0;
 430   if(IsUseMCInfo())
 431   {
 432     //
 433     if(!mcEvent) {
 434       AliDebug(AliLog::kError, "mcEvent not available");
 435       return;
 436     }
 437     // get MC event header
 438     header = mcEvent->Header();
 439     if (!header) {
 440       AliDebug(AliLog::kError, "Header not available");
 441       return;
 442     }
 443     // MC particle stack
 444     stack = mcEvent->Stack();
 445     if (!stack) {
 446       AliDebug(AliLog::kError, "Stack not available");
 447       return;
 448     }
 449
 450     // get MC vertex
 451     genHeader = header->GenEventHeader();
 452     if (!genHeader) {
 453       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 454       return;
 455     }
 456     genHeader->PrimaryVertex(vtxMC);
 457
 458     // multipliticy of all MC primary tracks
 459     // in Zv, pt and eta ranges)
 460     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
 461
 462   } // end bUseMC
 463
 464   // get reconstructed vertex
 465   //const AliESDVertex* vtxESD = 0;
 466   AliESDVertex* vtxESD = 0;
 467   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
 468         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
 469   }
 470   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
 471      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
 472   }
 473   else {
 474         return;
 475   }
 476
 477   if(!vtxESD) return;
 478
 479   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
 480   //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZv());
 481   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
 482
 483
 484   // check event cuts
 485   if(isEventOK && isEventTriggered)
 486   {
 487
 488     //
 489     // get IR information
 490     //
 491     AliESDHeader *esdHeader = 0;
 492     esdHeader = esdEvent->GetHeader();
 493     if(!esdHeader) return;
 494     //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
 495     //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
 496
 497     // Use when Peter commit the changes in the header
 498     Int_t ir1 = 0;
 499     Int_t ir2 = 0;
 500
 501     //
 502     Double_t vert[3] = {0};
 503     vert[0] = vtxESD->GetXv();
 504     vert[1] = vtxESD->GetYv();
 505     vert[2] = vtxESD->GetZv();
 506     Int_t mult = vtxESD->GetNContributors();
 507     Double_t bz = esdEvent->GetMagneticField();
 508     Double_t runNumber = esdEvent->GetRunNumber();
 509     Double_t evtTimeStamp = esdEvent->GetTimeStamp();
 510     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
 511
 512     // high pT tracks
 513     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 514     {
 515       AliESDtrack *track = esdEvent->GetTrack(iTrack);
 516       if(!track) continue;
 517       if(track->Charge()==0) continue;
 518       if(!esdTrackCuts->AcceptTrack(track)) continue;
 519       if(!accCuts->AcceptTrack(track)) continue;
 520
 521       // downscale low-pT tracks
 522       Double_t scalempt= TMath::Min(track->Pt(),10.);
 523       Double_t downscaleF = gRandom->Rndm();
 524       downscaleF *= fLowPtTrackDownscaligF;
 525       if(TMath::Exp(2*scalempt)<downscaleF) continue;
 526       //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
 527
 528       AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
 529       if (!tpcInner) continue;
 530       // transform to the track reference frame
 531       Bool_t isOK = kFALSE;
 532       isOK = tpcInner->Rotate(track->GetAlpha());
 533       isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
 534       if(!isOK) continue;
 535
 536       // Dump to the tree
 537       // vertex
 538       // TPC-ITS tracks
 539       //
 540       if(!fTreeSRedirector) return;
 541       (*fTreeSRedirector)<<"highPt"<<
 542         "fileName.="<<&fileName<<
 543         "runNumber="<<runNumber<<
 544         "evtTimeStamp="<<evtTimeStamp<<
 545         "evtNumberInFile="<<evtNumberInFile<<
 546         "Bz="<<bz<<
 547         "vtxESD.="<<vtxESD<<
 548         "IRtot="<<ir1<<
 549         "IRint2="<<ir2<<
 550         "mult="<<mult<<
 551         "esdTrack.="<<track<<
 552         "centralityF="<<centralityF<<
 553         "\n";
 554     }
 555   }
 556
 557 }
 558
 559
 560 //_____________________________________________________________________________
 561 void AlidNdPtTrackDumpTask::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const /*esdFriend*/)
 562 {
 563   //
 564   // Process laser events
 565   //
 566   if(!esdEvent) {
 567     AliDebug(AliLog::kError, "esdEvent not available");
 568     return;
 569   }
 570
 571   // get file name
 572   TTree *chain = (TChain*)GetInputData(0);
 573   if(!chain) {
 574     Printf("ERROR: Could not receive input chain");
 575     return;
 576   }
 577   TObjString fileName(chain->GetCurrentFile()->GetName());
 578
 579   // laser events
 580   const AliESDHeader* esdHeader = esdEvent->GetHeader();
 581   if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib)
 582   {
 583     Int_t countLaserTracks = 0;
 584     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 585     {
 586       AliESDtrack *track = esdEvent->GetTrack(iTrack);
 587       if(!track) continue;
 588
 589       if(track->GetTPCInnerParam()) countLaserTracks++;
 590     }
 591
 592     if(countLaserTracks > 100)
 593     {
 594       Double_t runNumber = esdEvent->GetRunNumber();
 595       Double_t evtTimeStamp = esdEvent->GetTimeStamp();
 596       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
 597       Double_t bz = esdEvent->GetMagneticField();
 598
 599       if(!fTreeSRedirector) return;
 600       (*fTreeSRedirector)<<"Laser"<<
 601         "fileName.="<<&fileName<<
 602         "runNumber="<<runNumber<<
 603         "evtTimeStamp="<<evtTimeStamp<<
 604         "evtNumberInFile="<<evtNumberInFile<<
 605         "Bz="<<bz<<
 606         "multTPCtracks="<<countLaserTracks<<
 607         "\n";
 608     }
 609   }
 610 }
 611
 612 //_____________________________________________________________________________
 613 void AlidNdPtTrackDumpTask::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
 614 {
 615   //
 616   // Select real events with high-pT tracks
 617   // Calculate and stor in the output tree:
 618   //  TPC constrained tracks
 619   //  InnerParams constrained tracks
 620   //  TPC-ITS tracks
 621   //  ITSout-InnerParams tracks
 622   //  chi2 distance between TPC constrained and TPC-ITS tracks
 623   //  chi2 distance between TPC refitted constrained and TPC-ITS tracks
 624   //  chi2 distance between ITSout and InnerParams tracks
 625   //  MC information:
 626   //   track references at ITSin, TPCin; InnerParam at first TPC track reference,
 627   //   particle ID, mother ID, production mechanism ...
 628   //
 629
 630   if(!esdEvent) {
 631     AliDebug(AliLog::kError, "esdEvent not available");
 632     return;
 633   }
 634
 635   // get selection cuts
 636   AlidNdPtEventCuts *evtCuts = GetEventCuts();
 637   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
 638   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
 639
 640   if(!evtCuts || !accCuts  || !esdTrackCuts) {
 641     AliDebug(AliLog::kError, "cuts not available");
 642     return;
 643   }
 644
 645   // trigger selection
 646   Bool_t isEventTriggered = kTRUE;
 647   AliPhysicsSelection *physicsSelection = NULL;
 648   AliTriggerAnalysis* triggerAnalysis = NULL;
 649
 650   //
 651   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 652   if (!inputHandler)
 653   {
 654     Printf("ERROR: Could not receive input handler");
 655     return;
 656   }
 657
 658   // get file name
 659   TTree *chain = (TChain*)GetInputData(0);
 660   if(!chain) {
 661     Printf("ERROR: Could not receive input chain");
 662     return;
 663   }
 664   TObjString fileName(chain->GetCurrentFile()->GetName());
 665
 666   // trigger
 667   if(evtCuts->IsTriggerRequired())
 668   {
 669     // always MB
 670     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
 671
 672     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 673     if(!physicsSelection) return;
 674     //SetPhysicsTriggerSelection(physicsSelection);
 675
 676     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
 677       // set trigger (V0AND)
 678       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 679       if(!triggerAnalysis) return;
 680       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
 681     }
 682   }
 683
 684   // centrality determination
 685   Float_t centralityF = -1;
 686   AliCentrality *esdCentrality = esdEvent->GetCentrality();
 687   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
 688
 689   // use MC information
 690   AliHeader* header = 0;
 691   AliGenEventHeader* genHeader = 0;
 692   AliStack* stack = 0;
 693   TArrayF vtxMC(3);
 694
 695   Int_t multMCTrueTracks = 0;
 696   if(IsUseMCInfo())
 697   {
 698     //
 699     if(!mcEvent) {
 700       AliDebug(AliLog::kError, "mcEvent not available");
 701       return;
 702     }
 703     // get MC event header
 704     header = mcEvent->Header();
 705     if (!header) {
 706       AliDebug(AliLog::kError, "Header not available");
 707       return;
 708     }
 709     // MC particle stack
 710     stack = mcEvent->Stack();
 711     if (!stack) {
 712       AliDebug(AliLog::kError, "Stack not available");
 713       return;
 714     }
 715
 716     // get MC vertex
 717     genHeader = header->GenEventHeader();
 718     if (!genHeader) {
 719       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 720       return;
 721     }
 722     genHeader->PrimaryVertex(vtxMC);
 723
 724     // multipliticy of all MC primary tracks
 725     // in Zv, pt and eta ranges)
 726     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
 727
 728   } // end bUseMC
 729
 730   // get reconstructed vertex
 731   //const AliESDVertex* vtxESD = 0;
 732   AliESDVertex* vtxESD = 0;
 733   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
 734         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
 735   }
 736   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
 737      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
 738   }
 739   else {
 740         return;
 741   }
 742
 743   if(!vtxESD) return;
 744
 745   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
 746   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
 747   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
 748
 749
 750   // check event cuts
 751   if(isEventOK && isEventTriggered)
 752   {
 753     //
 754     // get IR information
 755     //
 756     AliESDHeader *esdHeader = 0;
 757     esdHeader = esdEvent->GetHeader();
 758     if(!esdHeader) return;
 759     //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
 760     //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
 761     //
 762     Int_t ir1 = 0;
 763     Int_t ir2 = 0;
 764
 765     //
 766     Double_t vert[3] = {0};
 767     vert[0] = vtxESD->GetXv();
 768     vert[1] = vtxESD->GetYv();
 769     vert[2] = vtxESD->GetZv();
 770     Int_t mult = vtxESD->GetNContributors();
 771     Double_t bz = esdEvent->GetMagneticField();
 772     Double_t runNumber = esdEvent->GetRunNumber();
 773     Double_t evtTimeStamp = esdEvent->GetTimeStamp();
 774     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
 775
 776     // high pT tracks
 777     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 778     {
 779       AliESDtrack *track = esdEvent->GetTrack(iTrack);
 780       if(!track) continue;
 781       if(track->Charge()==0) continue;
 782       if(!esdTrackCuts->AcceptTrack(track)) continue;
 783       if(!accCuts->AcceptTrack(track)) continue;
 784
 785       // downscale low-pT tracks
 786       Double_t scalempt= TMath::Min(track->Pt(),10.);
 787       Double_t downscaleF = gRandom->Rndm();
 788       downscaleF *= fLowPtTrackDownscaligF;
 789       if(TMath::Exp(2*scalempt)<downscaleF) continue;
 790       //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
 791
 792       // Dump to the tree
 793       // vertex
 794       // TPC constrained tracks
 795       // InnerParams constrained tracks
 796       // TPC-ITS tracks
 797       // ITSout-InnerParams tracks
 798       // chi2 distance between TPC constrained and TPC-ITS tracks
 799       // chi2 distance between TPC refitted constrained and TPC-ITS tracks
 800       // chi2 distance between ITSout and InnerParams tracks
 801       // MC information
 802
 803       Double_t x[3]; track->GetXYZ(x);
 804       Double_t b[3]; AliTracker::GetBxByBz(x,b);
 805
 806       //
 807       // Transform TPC inner params to track reference frame
 808       //
 809       Bool_t isOKtpcInner = kFALSE;
 810       AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
 811       if (tpcInner) {
 812         // transform to the track reference frame
 813         isOKtpcInner = tpcInner->Rotate(track->GetAlpha());
 814         isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
 815       }
 816
 817       //
 818       // Constrain TPC inner to vertex
 819       // clone TPCinner has to be deleted
 820       //
 821       Bool_t isOKtpcInnerC = kFALSE;
 822       AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
 823       if (tpcInnerC) {
 824         isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);
 825         isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());
 826         isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
 827       }
 828
 829       //
 830       // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex
 831       // Clone track InnerParams has to be deleted
 832       //
 833       Bool_t isOKtrackInnerC = kFALSE;
 834       AliExternalTrackParam * trackInnerC =  new AliExternalTrackParam(*(track->GetInnerParam()));
 835       if (trackInnerC) {
 836         isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);
 837         isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());
 838         isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
 839       }
 840
 841       //
 842       // calculate chi2 between vi and vj vectors
 843       // with covi and covj covariance matrices
 844       // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
 845       //
 846       TMatrixD deltaT(5,1), deltaTtrackC(5,1);
 847       TMatrixD delta(1,5),  deltatrackC(1,5);
 848       TMatrixD covarM(5,5), covarMtrackC(5,5);
 849       TMatrixD chi2(1,1);
 850       TMatrixD chi2trackC(1,1);
 851
 852       if(isOKtpcInnerC && isOKtrackInnerC)
 853       {
 854         for (Int_t ipar=0; ipar<5; ipar++) {
 855           deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
 856           delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
 857
 858           deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
 859           deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
 860
 861           for (Int_t jpar=0; jpar<5; jpar++) {
 862             Int_t index=track->GetIndex(ipar,jpar);
 863             covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
 864             covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];
 865           }
 866         }
 867
 868         // chi2 distance TPC constrained and TPC+ITS
 869         TMatrixD covarMInv = covarM.Invert();
 870         TMatrixD mat2 = covarMInv*deltaT;
 871         chi2 = delta*mat2;
 872         //chi2.Print();
 873
 874         // chi2 distance TPC refitted constrained and TPC+ITS
 875         TMatrixD covarMInvtrackC = covarMtrackC.Invert();
 876         TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;
 877         chi2trackC = deltatrackC*mat2trackC;
 878         //chi2trackC.Print();
 879       }
 880
 881
 882       //
 883       // Propagate ITSout to TPC inner wall
 884       // and calculate chi2 distance to track (InnerParams)
 885       //
 886       const Double_t kTPCRadius=85;
 887       const Double_t kStep=3;
 888
 889       // clone track InnerParams has to be deleted
 890       Bool_t isOKtrackInnerC2 = kFALSE;
 891       AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));
 892       if (trackInnerC2) {
 893         isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);
 894       }
 895
 896       Bool_t isOKouterITSc = kFALSE;
 897       AliExternalTrackParam *outerITSc = NULL;
 898       TMatrixD chi2OuterITS(1,1);
 899
 900       if(esdFriend && esdFriend->TestSkipBit()==kFALSE)
 901       {
 902         // propagate ITSout to TPC inner wall
 903         AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);
 904
 905         if(friendTrack)
 906         {
 907           outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut());
 908           if(outerITSc)
 909           {
 910             isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);
 911             isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());
 912             isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());
 913
 914             //
 915             // calculate chi2 between outerITS and innerParams
 916             // cov without z-coordinate at the moment
 917             //
 918             TMatrixD deltaTouterITS(4,1);
 919             TMatrixD deltaouterITS(1,4);
 920             TMatrixD covarMouterITS(4,4);
 921
 922             if(isOKtrackInnerC2 && isOKouterITSc) {
 923               Int_t kipar = 0;
 924               Int_t kjpar = 0;
 925               for (Int_t ipar=0; ipar<5; ipar++) {
 926                 if(ipar!=1) {
 927                   deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
 928                   deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
 929                 }
 930
 931                 kjpar=0;
 932                 for (Int_t jpar=0; jpar<5; jpar++) {
 933                   Int_t index=outerITSc->GetIndex(ipar,jpar);
 934                   if(ipar !=1 || jpar!=1) {
 935                     covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];
 936                   }
 937                   if(jpar!=1)  kjpar++;
 938                 }
 939                 if(ipar!=1) kipar++;
 940               }
 941
 942               // chi2 distance ITSout and InnerParams
 943               TMatrixD covarMInvouterITS = covarMouterITS.Invert();
 944               TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;
 945               chi2OuterITS = deltaouterITS*mat2outerITS;
 946               //chi2OuterITS.Print();
 947             }
 948           }
 949         }
 950       }
 951
 952       //
 953       // MC info
 954       //
 955       TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;
 956       TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;
 957       Int_t mech=-1, mechTPC=-1, mechITS=-1;
 958       Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;
 959       Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;
 960       Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;
 961       Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;
 962
 963       AliTrackReference *refTPCIn = NULL;
 964       AliTrackReference *refITS = NULL;
 965
 966       Bool_t isOKtrackInnerC3 = kFALSE;
 967       AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
 968
 969       if(IsUseMCInfo())
 970       {
 971         if(!stack) return;
 972
 973         //
 974         // global track
 975         //
 976         Int_t label = TMath::Abs(track->GetLabel());
 977         particle = stack->Particle(label);
 978         if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)
 979         {
 980           particleMother = GetMother(particle,stack);
 981           mech = particle->GetUniqueID();
 982           isPrim = stack->IsPhysicalPrimary(label);
 983           isFromStrangess  = IsFromStrangeness(label,stack);
 984           isFromConversion = IsFromConversion(label,stack);
 985           isFromMaterial   = IsFromMaterial(label,stack);
 986         }
 987
 988         //
 989         // TPC track
 990         //
 991         Int_t labelTPC = TMath::Abs(track->GetTPCLabel());
 992         particleTPC = stack->Particle(labelTPC);
 993         if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)
 994         {
 995           particleMotherTPC = GetMother(particleTPC,stack);
 996           mechTPC = particleTPC->GetUniqueID();
 997           isPrimTPC = stack->IsPhysicalPrimary(labelTPC);
 998           isFromStrangessTPC  = IsFromStrangeness(labelTPC,stack);
 999           isFromConversionTPC = IsFromConversion(labelTPC,stack);
1000           isFromMaterialTPC   = IsFromMaterial(labelTPC,stack);
1001         }
1002
1003         //
1004         // store first track reference
1005         // for TPC track
1006         //
1007         TParticle *part=0;
1008         TClonesArray *trefs=0;
1009         Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);
1010
1011         if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.)
1012         {
1013           Int_t nTrackRef = trefs->GetEntries();
1014           //printf("nTrackRef %d \n",nTrackRef);
1015
1016           Int_t countITS = 0;
1017           for (Int_t iref = 0; iref < nTrackRef; iref++)
1018           {
1019             AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);
1020
1021              // Ref. in the middle ITS
1022             if(ref && ref->DetectorId()==AliTrackReference::kITS)
1023             {
1024               if(!refITS && countITS==2) {
1025                 refITS = ref;
1026                 //printf("refITS %p \n",refITS);
1027               }
1028               countITS++;
1029             }
1030
1031             // TPC
1032             if(ref && ref->DetectorId()==AliTrackReference::kTPC)
1033             {
1034               if(!refTPCIn) {
1035                 refTPCIn = ref;
1036                 //printf("refTPCIn %p \n",refTPCIn);
1037                 //break;
1038               }
1039             }
1040           }
1041
1042           // transform inner params to TrackRef
1043           // reference frame
1044           if(refTPCIn && trackInnerC3)
1045           {
1046             Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());
1047             isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);
1048             isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);
1049           }
1050         }
1051
1052         //
1053         // ITS track
1054         //
1055         Int_t labelITS = TMath::Abs(track->GetITSLabel());
1056         particleITS = stack->Particle(labelITS);
1057         if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)
1058         {
1059           particleMotherITS = GetMother(particleITS,stack);
1060           mechITS = particleITS->GetUniqueID();
1061           isPrimITS = stack->IsPhysicalPrimary(labelITS);
1062           isFromStrangessITS  = IsFromStrangeness(labelITS,stack);
1063           isFromConversionITS = IsFromConversion(labelITS,stack);
1064           isFromMaterialITS   = IsFromMaterial(labelITS,stack);
1065         }
1066       }
1067
1068       //
1069       Bool_t dumpToTree=kFALSE;
1070
1071       if(isOKtpcInnerC  && isOKtrackInnerC) dumpToTree = kTRUE;
1072       if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1073       if(fUseMCInfo     && isOKtrackInnerC3) dumpToTree = kTRUE;
1074
1075       //
1076       if(fTreeSRedirector && dumpToTree)
1077       {
1078         (*fTreeSRedirector)<<"highPt"<<
1079           "fileName.="<<&fileName<<
1080           "runNumber="<<runNumber<<
1081           "evtTimeStamp="<<evtTimeStamp<<
1082           "evtNumberInFile="<<evtNumberInFile<<
1083           "Bz="<<bz<<
1084           "vtxESD.="<<vtxESD<<
1085           "IRtot="<<ir1<<
1086           "IRint2="<<ir2<<
1087           "mult="<<mult<<
1088           "esdTrack.="<<track<<
1089           "extTPCInnerC.="<<tpcInnerC<<
1090           "extInnerParamC.="<<trackInnerC<<
1091           "extInnerParam.="<<trackInnerC2<<
1092           "extOuterITS.="<<outerITSc<<
1093           "extInnerParamRef.="<<trackInnerC3<<
1094           "refTPCIn.="<<refTPCIn<<
1095           "refITS.="<<refITS<<
1096           "chi2TPCInnerC="<<chi2(0,0)<<
1097           "chi2InnerC="<<chi2trackC(0,0)<<
1098           "chi2OuterITS="<<chi2OuterITS(0,0)<<
1099           "centralityF="<<centralityF<<
1100           "particle.="<<particle<<
1101           "particleMother.="<<particleMother<<
1102           "mech="<<mech<<
1103           "isPrim="<<isPrim<<
1104           "isFromStrangess="<<isFromStrangess<<
1105           "isFromConversion="<<isFromConversion<<
1106           "isFromMaterial="<<isFromMaterial<<
1107           "particleTPC.="<<particleTPC<<
1108           "particleMotherTPC.="<<particleMotherTPC<<
1109           "mechTPC="<<mechTPC<<
1110           "isPrimTPC="<<isPrimTPC<<
1111           "isFromStrangessTPC="<<isFromStrangessTPC<<
1112           "isFromConversionTPC="<<isFromConversionTPC<<
1113           "isFromMaterialTPC="<<isFromMaterialTPC<<
1114           "particleITS.="<<particleITS<<
1115           "particleMotherITS.="<<particleMotherITS<<
1116           "mechITS="<<mechITS<<
1117           "isPrimITS="<<isPrimITS<<
1118           "isFromStrangessITS="<<isFromStrangessITS<<
1119           "isFromConversionITS="<<isFromConversionITS<<
1120           "isFromMaterialITS="<<isFromMaterialITS<<
1121           "\n";
1122         }
1123
1124         if(tpcInnerC) delete tpcInnerC;
1125         if(trackInnerC) delete trackInnerC;
1126         if(trackInnerC2) delete trackInnerC2;
1127         if(outerITSc) delete outerITSc;
1128         if(trackInnerC3) delete trackInnerC3;
1129     }
1130   }
1131
1132 }
1133
1134
1135 //_____________________________________________________________________________
1136 void AlidNdPtTrackDumpTask::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1137 {
1138   //
1139   // Fill tree for efficiency studies MC only
1140
1141   if(!esdEvent) {
1142     AliDebug(AliLog::kError, "esdEvent not available");
1143     return;
1144   }
1145
1146    if(!mcEvent) {
1147     AliDebug(AliLog::kError, "mcEvent not available");
1148     return;
1149   }
1150
1151   // get selection cuts
1152   AlidNdPtEventCuts *evtCuts = GetEventCuts();
1153   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
1154   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1155
1156   if(!evtCuts || !accCuts  || !esdTrackCuts) {
1157     AliDebug(AliLog::kError, "cuts not available");
1158     return;
1159   }
1160
1161   // trigger selection
1162   Bool_t isEventTriggered = kTRUE;
1163   AliPhysicsSelection *physicsSelection = NULL;
1164   AliTriggerAnalysis* triggerAnalysis = NULL;
1165
1166   //
1167   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1168   if (!inputHandler)
1169   {
1170     Printf("ERROR: Could not receive input handler");
1171     return;
1172   }
1173
1174   // get file name
1175   TTree *chain = (TChain*)GetInputData(0);
1176   if(!chain) {
1177     Printf("ERROR: Could not receive input chain");
1178     return;
1179   }
1180   TObjString fileName(chain->GetCurrentFile()->GetName());
1181
1182   // trigger
1183   if(evtCuts->IsTriggerRequired())
1184   {
1185     // always MB
1186     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1187
1188     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1189     if(!physicsSelection) return;
1190
1191     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1192       // set trigger (V0AND)
1193       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1194       if(!triggerAnalysis) return;
1195       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1196     }
1197   }
1198
1199   // centrality determination
1200   Float_t centralityF = -1;
1201   AliCentrality *esdCentrality = esdEvent->GetCentrality();
1202   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1203
1204   // use MC information
1205   AliHeader* header = 0;
1206   AliGenEventHeader* genHeader = 0;
1207   AliStack* stack = 0;
1208   TArrayF vtxMC(3);
1209
1210   Int_t multMCTrueTracks = 0;
1211   if(IsUseMCInfo())
1212   {
1213     //
1214     if(!mcEvent) {
1215       AliDebug(AliLog::kError, "mcEvent not available");
1216       return;
1217     }
1218     // get MC event header
1219     header = mcEvent->Header();
1220     if (!header) {
1221       AliDebug(AliLog::kError, "Header not available");
1222       return;
1223     }
1224     // MC particle stack
1225     stack = mcEvent->Stack();
1226     if (!stack) {
1227       AliDebug(AliLog::kError, "Stack not available");
1228       return;
1229     }
1230
1231     // get MC vertex
1232     genHeader = header->GenEventHeader();
1233     if (!genHeader) {
1234       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1235       return;
1236     }
1237     genHeader->PrimaryVertex(vtxMC);
1238
1239     // multipliticy of all MC primary tracks
1240     // in Zv, pt and eta ranges)
1241     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1242
1243   } // end bUseMC
1244
1245   // get reconstructed vertex
1246   //const AliESDVertex* vtxESD = 0;
1247   AliESDVertex* vtxESD = 0;
1248   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
1249         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1250   }
1251   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
1252      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1253   }
1254   else {
1255         return;
1256   }
1257
1258   if(!vtxESD) return;
1259
1260   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1261   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1262   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1263
1264   // check event cuts
1265   if(isEventOK && isEventTriggered)
1266   {
1267     if(IsUseMCInfo())
1268     {
1269       if(!stack) return;
1270
1271       //
1272       // MC info
1273       //
1274       TParticle *particle=NULL;
1275       TParticle *particleMother=NULL;
1276       Int_t mech=-1;
1277
1278       // reco event info
1279       Double_t vert[3] = {0};
1280       vert[0] = vtxESD->GetXv();
1281       vert[1] = vtxESD->GetYv();
1282       vert[2] = vtxESD->GetZv();
1283       Int_t mult = vtxESD->GetNContributors();
1284       Double_t bz = esdEvent->GetMagneticField();
1285       Double_t runNumber = esdEvent->GetRunNumber();
1286       Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1287       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1288
1289       // loop over MC stack
1290       for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
1291       {
1292          particle = stack->Particle(iMc);
1293          if (!particle)
1294          continue;
1295
1296          // only charged particles
1297          if(!particle->GetPDG()) continue;
1298          Double_t charge = particle->GetPDG()->Charge()/3.;
1299          if (TMath::Abs(charge) < 0.001)
1300          continue;
1301
1302          // only primary particles
1303          Bool_t prim = stack->IsPhysicalPrimary(iMc);
1304          if(!prim) continue;
1305
1306          // downscale low-pT particles
1307          Double_t scalempt= TMath::Min(particle->Pt(),10.);
1308          Double_t downscaleF = gRandom->Rndm();
1309          downscaleF *= fLowPtTrackDownscaligF;
1310          if(TMath::Exp(2*scalempt)<downscaleF) continue;
1311
1312          // is particle in acceptance
1313          if(!accCuts->AcceptTrack(particle)) continue;
1314
1315          // check if particle reconstructed
1316          Bool_t isRec = kFALSE;
1317          Int_t  trackIndex = -1;
1318          for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1319          {
1320
1321            AliESDtrack *track = esdEvent->GetTrack(iTrack);
1322            if(!track) continue;
1323            if(track->Charge()==0) continue;
1324            if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track))
1325            {
1326              Int_t label =  TMath::Abs(track->GetLabel());
1327              if(label == iMc) {
1328                isRec = kTRUE;
1329                trackIndex = iTrack;
1330                break;
1331              }
1332            }
1333          }
1334
1335          // Store information in the output tree
1336          AliESDtrack *recTrack = NULL;
1337          if(trackIndex>-1)  {
1338            recTrack = esdEvent->GetTrack(trackIndex);
1339          } else {
1340            recTrack = new AliESDtrack();
1341          }
1342
1343          particleMother = GetMother(particle,stack);
1344          mech = particle->GetUniqueID();
1345
1346          //MC particle track length
1347          Double_t tpcTrackLength = 0.;
1348          AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
1349          if(mcParticle) {
1350            Int_t counter;
1351            tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
1352          }
1353
1354
1355          //
1356          if(fTreeSRedirector) {
1357            (*fTreeSRedirector)<<"MCEffTree"<<
1358            "fileName.="<<&fileName<<
1359            "runNumber="<<runNumber<<
1360            "evtTimeStamp="<<evtTimeStamp<<
1361            "evtNumberInFile="<<evtNumberInFile<<
1362            "Bz="<<bz<<
1363            "vtxESD.="<<vtxESD<<
1364            "mult="<<mult<<
1365            "esdTrack.="<<recTrack<<
1366            "isRec="<<isRec<<
1367            "tpcTrackLength="<<tpcTrackLength<<
1368            "particle.="<<particle<<
1369            "particleMother.="<<particleMother<<
1370            "mech="<<mech<<
1371            "\n";
1372          }
1373
1374          if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
1375       }
1376     }
1377   }
1378
1379 }
1380
1381 //_____________________________________________________________________________
1382 Bool_t AlidNdPtTrackDumpTask::IsHighDeDxParticle(AliESDtrack * track) {
1383   //
1384   // check if particle is Z > 1
1385   //
1386   if (track->GetTPCNcls() < 60) return kFALSE;
1387   Double_t mom = track->GetInnerParam()->GetP();
1388   if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
1389   Float_t dca[2], bCov[3];
1390   track->GetImpactParameters(dca,bCov);
1391   //
1392
1393   Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
1394
1395   if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
1396
1397   return kFALSE;
1398 }
1399
1400 //_____________________________________________________________________________
1401 void AlidNdPtTrackDumpTask::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1402 {
1403   //
1404   // Select real events with V0 (K0s and Lambda) high-pT candidates
1405   //
1406   if(!esdEvent) {
1407     AliDebug(AliLog::kError, "esdEvent not available");
1408     return;
1409   }
1410
1411   // get selection cuts
1412   AlidNdPtEventCuts *evtCuts = GetEventCuts();
1413   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
1414   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1415
1416   if(!evtCuts || !accCuts  || !esdTrackCuts) {
1417     AliDebug(AliLog::kError, "cuts not available");
1418     return;
1419   }
1420
1421   // trigger selection
1422   Bool_t isEventTriggered = kTRUE;
1423   AliPhysicsSelection *physicsSelection = NULL;
1424   AliTriggerAnalysis* triggerAnalysis = NULL;
1425
1426   //
1427   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1428   if (!inputHandler)
1429   {
1430     Printf("ERROR: Could not receive input handler");
1431     return;
1432   }
1433
1434   // get file name
1435   TTree *chain = (TChain*)GetInputData(0);
1436   if(!chain) {
1437     Printf("ERROR: Could not receive input chain");
1438     return;
1439   }
1440   TObjString fileName(chain->GetCurrentFile()->GetName());
1441
1442   // trigger
1443   if(evtCuts->IsTriggerRequired())
1444   {
1445     // always MB
1446     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1447
1448     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1449     if(!physicsSelection) return;
1450     //SetPhysicsTriggerSelection(physicsSelection);
1451
1452     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1453       // set trigger (V0AND)
1454       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1455       if(!triggerAnalysis) return;
1456       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1457     }
1458   }
1459
1460   // centrality determination
1461   Float_t centralityF = -1;
1462   AliCentrality *esdCentrality = esdEvent->GetCentrality();
1463   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1464
1465
1466   // get reconstructed vertex
1467   //const AliESDVertex* vtxESD = 0;
1468   AliESDVertex* vtxESD = 0;
1469   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
1470         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1471   }
1472   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
1473      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1474   }
1475   else {
1476         return;
1477   }
1478
1479   if(!vtxESD) return;
1480
1481   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1482   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1483   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1484
1485   // check event cuts
1486   if(isEventOK && isEventTriggered) {
1487   //
1488   // Dump the pt downscaled V0 into the tree
1489   //
1490   Int_t ntracks = esdEvent->GetNumberOfTracks();
1491   Int_t nV0s = esdEvent->GetNumberOfV0s();
1492   Int_t run = esdEvent->GetRunNumber();
1493   Int_t time= esdEvent->GetTimeStamp();
1494   Int_t evNr=esdEvent->GetEventNumberInFile();
1495   Double_t bz = esdEvent->GetMagneticField();
1496
1497
1498   for (Int_t iv0=0; iv0<nV0s; iv0++){
1499     AliESDv0 * v0 = esdEvent->GetV0(iv0);
1500     if (!v0) continue;
1501     AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
1502     AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
1503     if (!track0) continue;
1504     if (!track1) continue;
1505     if (track0->GetSign()<0) {
1506       track1 = esdEvent->GetTrack(v0->GetIndex(0));
1507       track0 = esdEvent->GetTrack(v0->GetIndex(1));
1508     }
1509     //
1510     Bool_t isDownscaled = IsV0Downscaled(v0);
1511     if (isDownscaled) continue;
1512     AliKFParticle kfparticle; //
1513     Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
1514     if (type==0) continue;
1515
1516     if(!fTreeSRedirector) return;
1517     (*fTreeSRedirector)<<"V0s"<<
1518       "isDownscaled="<<isDownscaled<<
1519       "Bz="<<bz<<
1520       "fileName.="<<&fileName<<
1521       "runNumber="<<run<<
1522       "evtTimeStamp="<<time<<
1523       "evtNumberInFile="<<evNr<<
1524       "type="<<type<<
1525       "ntracks="<<ntracks<<
1526       "v0.="<<v0<<
1527       "kf.="<<&kfparticle<<
1528       "track0.="<<track0<<
1529       "track1.="<<track1<<
1530       "centralityF="<<centralityF<<
1531       "\n";
1532   }
1533   }
1534 }
1535
1536 //_____________________________________________________________________________
1537 void AlidNdPtTrackDumpTask::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1538 {
1539   //
1540   // Select real events with large TPC dEdx signal
1541   //
1542   if(!esdEvent) {
1543     AliDebug(AliLog::kError, "esdEvent not available");
1544     return;
1545   }
1546
1547   // get selection cuts
1548   AlidNdPtEventCuts *evtCuts = GetEventCuts();
1549   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
1550   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1551
1552   if(!evtCuts || !accCuts  || !esdTrackCuts) {
1553     AliDebug(AliLog::kError, "cuts not available");
1554     return;
1555   }
1556
1557   // get file name
1558   TTree *chain = (TChain*)GetInputData(0);
1559   if(!chain) {
1560     Printf("ERROR: Could not receive input chain");
1561     return;
1562   }
1563   TObjString fileName(chain->GetCurrentFile()->GetName());
1564
1565   // trigger
1566   Bool_t isEventTriggered = kTRUE;
1567   AliPhysicsSelection *physicsSelection = NULL;
1568   AliTriggerAnalysis* triggerAnalysis = NULL;
1569
1570   //
1571   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1572   if (!inputHandler)
1573   {
1574     Printf("ERROR: Could not receive input handler");
1575     return;
1576   }
1577
1578   if(evtCuts->IsTriggerRequired())
1579   {
1580     // always MB
1581     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1582
1583     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1584     if(!physicsSelection) return;
1585
1586     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1587       // set trigger (V0AND)
1588       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1589       if(!triggerAnalysis) return;
1590       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1591     }
1592   }
1593
1594   // get reconstructed vertex
1595   AliESDVertex* vtxESD = 0;
1596   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
1597         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1598   }
1599   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
1600      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1601   }
1602   else {
1603         return;
1604   }
1605   if(!vtxESD) return;
1606
1607   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1608   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1609   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1610
1611
1612   // check event cuts
1613   if(isEventOK && isEventTriggered)
1614   {
1615     Double_t vert[3] = {0};
1616     vert[0] = vtxESD->GetXv();
1617     vert[1] = vtxESD->GetYv();
1618     vert[2] = vtxESD->GetZv();
1619     Int_t mult = vtxESD->GetNContributors();
1620     Double_t bz = esdEvent->GetMagneticField();
1621     Double_t runNumber = esdEvent->GetRunNumber();
1622     Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1623     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1624
1625     // large dEdx
1626     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1627     {
1628       AliESDtrack *track = esdEvent->GetTrack(iTrack);
1629       if(!track) continue;
1630       if(track->Charge()==0) continue;
1631       if(!esdTrackCuts->AcceptTrack(track)) continue;
1632       if(!accCuts->AcceptTrack(track)) continue;
1633
1634       if(!IsHighDeDxParticle(track)) continue;
1635
1636       if(!fTreeSRedirector) return;
1637       (*fTreeSRedirector)<<"dEdx"<<
1638       "fileName.="<<&fileName<<
1639       "runNumber="<<runNumber<<
1640       "evtTimeStamp="<<evtTimeStamp<<
1641       "evtNumberInFile="<<evtNumberInFile<<
1642       "Bz="<<bz<<
1643       "vtxESD.="<<vtxESD<<
1644       "mult="<<mult<<
1645       "esdTrack.="<<track<<
1646       "\n";
1647     }
1648   }
1649 }
1650
1651 //_____________________________________________________________________________
1652 Int_t   AlidNdPtTrackDumpTask::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
1653 {
1654   //
1655   // Create KF particle in case the V0 fullfill selection criteria
1656   //
1657   // Selection criteria
1658   //  0. algorithm cut
1659   //  1. track cut
1660   //  3. chi2 cut
1661   //  4. rough mass cut
1662   //  5. Normalized pointing angle cut
1663   //
1664   const Double_t cutMass=0.2;
1665   const Double_t kSigmaDCACut=3;
1666   //
1667   // 0.) algo cut - accept only on the fly
1668   //
1669   if (v0->GetOnFlyStatus() ==kFALSE) return 0;
1670   //
1671   // 1.) track cut
1672   //
1673   AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
1674   AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
1675   /*
1676     TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
1677     TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
1678     TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
1679   */
1680   if (TMath::Abs(track0->GetTgl())>1) return 0;
1681   if (TMath::Abs(track1->GetTgl())>1) return 0;
1682   if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
1683   if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
1684   //if ((track0->GetITSclusters(0))<2) return 0;
1685   //if ((track1->GetITSclusters(0))<2) return 0;
1686   Float_t pos0[2]={0}, cov0[3]={0};
1687   Float_t pos1[2]={0}, cov1[3]={0};
1688   track0->GetImpactParameters(pos0,cov0);
1689   track0->GetImpactParameters(pos1,cov1);
1690   //
1691   if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
1692   if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
1693   //
1694   //
1695   // 3.) Chi2 cut
1696   //
1697   Double_t chi2KF = v0->GetKFInfo(2,2,2);
1698   if (chi2KF>25) return 0;
1699   //
1700   // 4.) Rough mass cut - 0.200 GeV
1701   //
1702   static Double_t masses[2]={-1};
1703   if (masses[0]<0){
1704     masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
1705     masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
1706   }
1707   Double_t mass00=  v0->GetEffMass(0,0);
1708   Double_t mass22=  v0->GetEffMass(2,2);
1709   Double_t mass42=  v0->GetEffMass(4,2);
1710   Double_t mass24=  v0->GetEffMass(2,4);
1711   Bool_t massOK=kFALSE;
1712   Int_t type=0;
1713   Int_t ptype=0;
1714   Double_t dmass=1;
1715   Int_t p1=0, p2=0;
1716   if (TMath::Abs(mass00-0)<cutMass) {
1717     massOK=kTRUE; type+=1;
1718     if (TMath::Abs(mass00-0)<dmass) {
1719       ptype=1;
1720       dmass=TMath::Abs(mass00-0);
1721       p1=0; p2=0;
1722     }
1723   }
1724   if (TMath::Abs(mass24-masses[1])<cutMass) {
1725     massOK=kTRUE; type+=2;
1726     if (TMath::Abs(mass24-masses[1])<dmass){
1727       dmass = TMath::Abs(mass24-masses[1]);
1728       ptype=2;
1729       p1=2; p2=4;
1730     }
1731   }
1732   if (TMath::Abs(mass42-masses[1])<cutMass) {
1733     massOK=kTRUE; type+=4;
1734     if (TMath::Abs(mass42-masses[1])<dmass){
1735       dmass = TMath::Abs(mass42-masses[1]);
1736       ptype=4;
1737       p1=4; p2=2;
1738     }
1739   }
1740   if (TMath::Abs(mass22-masses[0])<cutMass) {
1741     massOK=kTRUE; type+=8;
1742     if (TMath::Abs(mass22-masses[0])<dmass){
1743       dmass = TMath::Abs(mass22-masses[0]);
1744       ptype=8;
1745       p1=2; p2=2;
1746     }
1747   }
1748   if (type==0) return 0;
1749   //
1750   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
1751   const AliExternalTrackParam *paramP = v0->GetParamP();
1752   const AliExternalTrackParam *paramN = v0->GetParamN();
1753   if (paramP->GetSign()<0){
1754     paramP=v0->GetParamP();
1755     paramN=v0->GetParamN();
1756   }
1757   //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
1758   //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
1759   //
1760   AliKFParticle kfp1( *paramP, spdg[p1]  );
1761   AliKFParticle kfp2( *paramN, -1 *spdg[p2]  );
1762   AliKFParticle V0KF;
1763   (V0KF)+=kfp1;
1764   (V0KF)+=kfp2;
1765   kfparticle=V0KF;
1766   //
1767   // Pointing angle
1768   //
1769   Double_t  errPhi    = V0KF.GetErrPhi();
1770   Double_t  pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
1771   if (pointAngle/errPhi>10) return 0;
1772   //
1773   return ptype;
1774 }
1775
1776 //_____________________________________________________________________________
1777 Bool_t AlidNdPtTrackDumpTask::IsV0Downscaled(AliESDv0 *const v0)
1778 {
1779   //
1780   // Downscale randomly low pt V0
1781   //
1782   //return kFALSE;
1783   Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
1784   Double_t scalempt= TMath::Min(maxPt,10.);
1785   Double_t downscaleF = gRandom->Rndm();
1786   downscaleF *= fLowPtV0DownscaligF;
1787
1788   //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
1789   if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;
1790   return kFALSE;
1791
1792   /*
1793     TH1F his1("his1","his1",100,0,10);
1794     TH1F his2("his2","his2",100,0,10);
1795     {for (Int_t i=0; i<10000; i++){
1796        Double_t rnd=gRandom->Exp(1);
1797        Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
1798        his1->Fill(rnd);
1799        if (!isDownscaled) his2->Fill(rnd);
1800     }}
1801
1802    */
1803
1804 }
1805
1806
1807
1808 //_____________________________________________________________________________
1809 Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
1810 {
1811  // Constrain TPC inner params constrained
1812  //
1813       if(!tpcInnerC) return kFALSE;
1814       if(!vtx) return kFALSE;
1815
1816       Double_t dz[2],cov[3];
1817       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
1818       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
1819       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
1820       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
1821
1822
1823       Double_t covar[6]; vtx->GetCovMatrix(covar);
1824       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
1825       Double_t c[3]={covar[2],0.,covar[5]};
1826       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
1827       if (chi2C>kVeryBig) return kFALSE;
1828
1829       if(!tpcInnerC->Update(p,c)) return kFALSE;
1830
1831   return kTRUE;
1832 }
1833
1834 //_____________________________________________________________________________
1835 Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
1836 {
1837  // Constrain TPC inner params constrained
1838  //
1839       if(!trackInnerC) return kFALSE;
1840       if(!vtx) return kFALSE;
1841
1842       const Double_t kRadius  = 2.8;
1843       const Double_t kMaxStep = 1.0;
1844
1845       Double_t dz[2],cov[3];
1846
1847       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
1848       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
1849       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
1850
1851       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
1852       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
1853
1854       //
1855       Double_t covar[6]; vtx->GetCovMatrix(covar);
1856       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
1857       Double_t c[3]={covar[2],0.,covar[5]};
1858       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
1859       if (chi2C>kVeryBig) return kFALSE;
1860
1861       if(!trackInnerC->Update(p,c)) return kFALSE;
1862
1863   return kTRUE;
1864 }
1865
1866
1867 //_____________________________________________________________________________
1868 TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack)
1869 {
1870   if(!particle) return NULL;
1871   if(!stack) return NULL;
1872
1873   Int_t motherLabel = TMath::Abs(particle->GetMother(0));
1874   TParticle* mother = NULL;
1875   mother = stack->Particle(motherLabel);
1876
1877 return mother;
1878 }
1879
1880 //_____________________________________________________________________________
1881 Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack)
1882 {
1883   Bool_t isFromConversion = kFALSE;
1884
1885   if(stack) {
1886     TParticle* particle = stack->Particle(label);
1887
1888     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
1889     {
1890        Int_t mech = particle->GetUniqueID(); // production mechanism
1891        Bool_t isPrim = stack->IsPhysicalPrimary(label);
1892
1893        Int_t motherLabel = TMath::Abs(particle->GetMother(0));
1894        TParticle* mother = stack->Particle(motherLabel);
1895        if(mother) {
1896           Int_t motherPdg = mother->GetPdgCode();
1897
1898           if(!isPrim && mech==5 && motherPdg==kGamma) {
1899             isFromConversion=kTRUE;
1900           }
1901        }
1902     }
1903   }
1904
1905   return isFromConversion;
1906 }
1907
1908 //_____________________________________________________________________________
1909 Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack)
1910 {
1911   Bool_t isFromMaterial = kFALSE;
1912
1913   if(stack) {
1914     TParticle* particle = stack->Particle(label);
1915
1916     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
1917     {
1918        Int_t mech = particle->GetUniqueID(); // production mechanism
1919        Bool_t isPrim = stack->IsPhysicalPrimary(label);
1920
1921        Int_t motherLabel = TMath::Abs(particle->GetMother(0));
1922        TParticle* mother = stack->Particle(motherLabel);
1923        if(mother) {
1924           if(!isPrim && mech==13) {
1925             isFromMaterial=kTRUE;
1926           }
1927        }
1928      }
1929   }
1930
1931   return isFromMaterial;
1932 }
1933
1934 //_____________________________________________________________________________
1935 Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack)
1936 {
1937   Bool_t isFromStrangeness = kFALSE;
1938
1939   if(stack) {
1940     TParticle* particle = stack->Particle(label);
1941
1942     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
1943     {
1944        Int_t mech = particle->GetUniqueID(); // production mechanism
1945        Bool_t isPrim = stack->IsPhysicalPrimary(label);
1946
1947        Int_t motherLabel = TMath::Abs(particle->GetMother(0));
1948        TParticle* mother = stack->Particle(motherLabel);
1949        if(mother) {
1950           Int_t motherPdg = mother->GetPdgCode();
1951
1952           // K+-, lambda, antilambda, K0s decays
1953           if(!isPrim && mech==4 &&
1954               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
1955           {
1956             isFromStrangeness = kTRUE;
1957           }
1958        }
1959     }
1960   }
1961
1962   return isFromStrangeness;
1963 }
1964
1965
1966 //_____________________________________________________________________________
1967 void AlidNdPtTrackDumpTask::FinishTaskOutput()
1968 {
1969   //
1970   // Called one at the end
1971   // locally on working node
1972   //
1973
1974   // must be deleted to store trees
1975   if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;
1976
1977   // open temporary file and copy trees to the ouptut container
1978
1979   TChain* chain = 0;
1980   /*
1981   TTree* tree1 = 0;
1982   TTree* tree2 = 0;
1983   TTree* tree3 = 0;
1984   TTree* tree4 = 0;
1985   TTree* tree5 = 0;
1986   TTree* tree6 = 0;
1987   */
1988   //
1989   chain = new TChain("highPt");
1990   if(chain) {
1991     chain->Add("jotwinow_Temp_Trees.root");
1992     fTree1 = chain->CopyTree("1");
1993     delete chain; chain=0;
1994   }
1995   if(fTree1) fTree1->Print();
1996
1997   //
1998   chain = new TChain("V0s");
1999   if(chain) {
2000     chain->Add("jotwinow_Temp_Trees.root");
2001     fTree2 = chain->CopyTree("1");
2002     delete chain; chain=0;
2003   }
2004   if(fTree2) fTree2->Print();
2005
2006   //
2007   chain = new TChain("dEdx");
2008   if(chain) {
2009     chain->Add("jotwinow_Temp_Trees.root");
2010     fTree3 = chain->CopyTree("1");
2011     delete chain; chain=0;
2012   }
2013   if(fTree3) fTree3->Print();
2014
2015   //
2016   chain = new TChain("Laser");
2017   if(chain) {
2018     chain->Add("jotwinow_Temp_Trees.root");
2019     fTree4 = chain->CopyTree("1");
2020     delete chain; chain=0;
2021   }
2022   if(fTree4) fTree4->Print();
2023
2024   //
2025   chain = new TChain("MCEffTree");
2026   if(chain) {
2027     chain->Add("jotwinow_Temp_Trees.root");
2028     fTree5 = chain->CopyTree("1");
2029     delete chain; chain=0;
2030   }
2031   if(fTree5) fTree5->Print();
2032
2033   //
2034   chain = new TChain("CosmicPairs");
2035   if(chain) {
2036     chain->Add("jotwinow_Temp_Trees.root");
2037     fTree6 = chain->CopyTree("1");
2038     delete chain; chain=0;
2039   }
2040   if(fTree6) fTree6->Print();
2041
2042
2043   OpenFile(1);
2044
2045   // Post output data.
2046   PostData(1, fTree1);
2047   PostData(2, fTree2);
2048   PostData(3, fTree3);
2049   PostData(4, fTree4);
2050   PostData(5, fTree5);
2051   PostData(6, fTree6);
2052 }
2053
2054 //_____________________________________________________________________________
2055 void AlidNdPtTrackDumpTask::Terminate(Option_t *)
2056 {
2057   // Called one at the end
2058   /*
2059   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));
2060   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;
2061   TChain* chain = new TChain("highPt");
2062   if(!chain) return;
2063   chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");
2064   TTree *tree = chain->CopyTree("1");
2065   if (chain) { delete chain; chain=0; }
2066   if(!tree) return;
2067   tree->Print();
2068   fOutputSummary = tree;
2069
2070   if (!fOutputSummary) {
2071     Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable %p \n", GetOutputData(1));
2072     return;
2073   }
2074   */
2075 }