PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtTrackDumpTask.cxx

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