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 \r
  20 #include "TChain.h"\r
  21 #include "TTreeStream.h"\r
  22 #include "TTree.h"\r
  23 #include "TH1F.h"\r
  24 #include "TCanvas.h"\r
  25 #include "TList.h"\r
  26 #include "TFile.h"\r
  27 #include "TMatrixD.h"\r
  28 #include "TRandom.h"\r
  29 \r
  30 #include "AliHeader.h"  \r
  31 #include "AliGenEventHeader.h"  \r
  32 #include "AliInputEventHandler.h"  \r
  33 #include "AliStack.h"  \r
  34 #include "AliTrackReference.h"  \r
  35 \r
  36 #include "AliPhysicsSelection.h"\r
  37 #include "AliAnalysisTask.h"\r
  38 #include "AliAnalysisManager.h"\r
  39 #include "AliESDEvent.h"\r
  40 #include "AliESDfriend.h"\r
  41 #include "AliMCEvent.h"\r
  42 #include "AliESDInputHandler.h"\r
  43 #include "AliESDVertex.h"\r
  44 #include "AliTracker.h"\r
  45 #include "AliGeomManager.h"\r
  46 \r
  47 #include "AliCentrality.h"\r
  48 #include "AliESDVZERO.h"\r
  49 #include "AliMultiplicity.h"\r
  50 \r
  51 #include "AliESDtrackCuts.h"\r
  52 #include "AliMCEventHandler.h"\r
  53 #include "dNdPt/AlidNdPt.h"\r
  54 #include "dNdPt/AlidNdPtEventCuts.h"\r
  55 #include "dNdPt/AlidNdPtAcceptanceCuts.h"\r
  56 \r
  57 #include "dNdPt/AlidNdPtTrackDumpTask.h"\r
  58 \r
  59 using namespace std;\r
  60 \r
  61 ClassImp(AlidNdPtTrackDumpTask)\r
  62 \r
  63 //_____________________________________________________________________________\r
  64 AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name) \r
  65   : AliAnalysisTaskSE(name)\r
  66   , fESD(0)\r
  67   , fMC(0)\r
  68   , fESDfriend(0)\r
  69   , fOutput(0)\r
  70   , fPitList(0)\r
  71   , fUseMCInfo(kFALSE)\r
  72   , fdNdPtEventCuts(0)\r
  73   , fdNdPtAcceptanceCuts(0)\r
  74   , fdNdPtRecAcceptanceCuts(0)\r
  75   , fEsdTrackCuts(0)\r
  76   , fTrigger(AliTriggerAnalysis::kMB1) \r
  77   , fAnalysisMode(AlidNdPtHelper::kTPC) \r
  78   , fOutputSummary(0)\r
  79   , fTreeSRedirector(0)\r
  80   , fCentralityEstimator(0)\r
  81 {\r
  82   // Constructor\r
  83 \r
  84   // Define input and output slots here\r
  85   DefineOutput(0, TTree::Class());\r
  86   //DefineOutput(1, TList::Class());\r
  87 \r
  88 }\r
  89 \r
  90 //_____________________________________________________________________________\r
  91 AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()\r
  92 {\r
  93   if(fOutput) delete fOutput;  fOutput =0; \r
  94   //if(fOutputSummary) delete fOutputSummary;  fOutputSummary =0; \r
  95   if(fTreeSRedirector) delete fTreeSRedirector;  fTreeSRedirector =0; \r
  96 \r
  97   if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL; \r
  98   if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;\r
  99   if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL;  \r
 100   if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;\r
 101 }\r
 102 \r
 103 //____________________________________________________________________________\r
 104 Bool_t AlidNdPtTrackDumpTask::Notify()\r
 105 {\r
 106   static Int_t count = 0;\r
 107   count++;\r
 108   TTree *chain = (TChain*)GetInputData(0);\r
 109   if(chain)\r
 110   {\r
 111     Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());\r
 112   }\r
 113 \r
 114 return kTRUE;\r
 115 }\r
 116 \r
 117 //_____________________________________________________________________________\r
 118 void AlidNdPtTrackDumpTask::UserCreateOutputObjects()\r
 119 {\r
 120   // Create histograms\r
 121   // Called once\r
 122   fOutput = new TList;\r
 123   fOutput->SetOwner();\r
 124 \r
 125   //\r
 126   // create output tree\r
 127   //\r
 128   fTreeSRedirector = new TTreeSRedirector("dNdPtOutliersAnalysisPbPb.root");\r
 129 \r
 130   PostData(0, fOutputSummary);\r
 131   //PostData(1, fOutput);\r
 132 }\r
 133 \r
 134 //_____________________________________________________________________________\r
 135 void AlidNdPtTrackDumpTask::UserExec(Option_t *) \r
 136 {\r
 137   //\r
 138   // Called for each event\r
 139   //\r
 140 \r
 141   // ESD event\r
 142   fESD = (AliESDEvent*) (InputEvent());\r
 143   if (!fESD) {\r
 144     Printf("ERROR: ESD event not available");\r
 145     return;\r
 146   }\r
 147 \r
 148   // MC event\r
 149   if(fUseMCInfo) {\r
 150     fMC = MCEvent();\r
 151     if (!fMC) {\r
 152       Printf("ERROR: MC event not available");\r
 153       return;\r
 154     }\r
 155   }\r
 156 \r
 157   fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
 158   if(!fESDfriend) {\r
 159     Printf("ERROR: ESD friends not available");\r
 160   }\r
 161 \r
 162   //\r
 163   Process(fESD,fMC,fESDfriend);\r
 164 \r
 165   // Post output data.\r
 166   PostData(0, fOutputSummary);\r
 167   //PostData(1, fOutput);\r
 168 }\r
 169 \r
 170 //_____________________________________________________________________________\r
 171 void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
 172 {\r
 173   //\r
 174   // Process real and/or simulated events\r
 175   //\r
 176   if(!esdEvent) {\r
 177     AliDebug(AliLog::kError, "esdEvent not available");\r
 178     return;\r
 179   }\r
 180 \r
 181   // get selection cuts\r
 182   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
 183   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
 184   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
 185 \r
 186   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
 187     AliDebug(AliLog::kError, "cuts not available");\r
 188     return;\r
 189   }\r
 190 \r
 191   // trigger selection\r
 192   Bool_t isEventTriggered = kTRUE;\r
 193   AliPhysicsSelection *physicsSelection = NULL;\r
 194   AliTriggerAnalysis* triggerAnalysis = NULL;\r
 195 \r
 196   // \r
 197   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
 198   if (!inputHandler)\r
 199   {\r
 200     Printf("ERROR: Could not receive input handler");\r
 201     return;\r
 202   }\r
 203 \r
 204   if(evtCuts->IsTriggerRequired())  \r
 205   {\r
 206     // always MB\r
 207     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
 208 \r
 209     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
 210     if(!physicsSelection) return;\r
 211     //SetPhysicsTriggerSelection(physicsSelection);\r
 212 \r
 213     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
 214       // set trigger (V0AND)\r
 215       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
 216       if(!triggerAnalysis) return;\r
 217       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
 218     }\r
 219   }\r
 220 \r
 221   // centrality determination\r
 222   Float_t centralityF = -1;\r
 223   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
 224   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 225 \r
 226   // use MC information\r
 227   AliHeader* header = 0;\r
 228   AliGenEventHeader* genHeader = 0;\r
 229   AliStack* stack = 0;\r
 230   TArrayF vtxMC(3);\r
 231 \r
 232   Int_t multMCTrueTracks = 0;\r
 233   if(IsUseMCInfo())\r
 234   {\r
 235     //\r
 236     if(!mcEvent) {\r
 237       AliDebug(AliLog::kError, "mcEvent not available");\r
 238       return;\r
 239     }\r
 240     // get MC event header\r
 241     header = mcEvent->Header();\r
 242     if (!header) {\r
 243       AliDebug(AliLog::kError, "Header not available");\r
 244       return;\r
 245     }\r
 246     // MC particle stack\r
 247     stack = mcEvent->Stack();\r
 248     if (!stack) {\r
 249       AliDebug(AliLog::kError, "Stack not available");\r
 250       return;\r
 251     }\r
 252 \r
 253     // get MC vertex\r
 254     genHeader = header->GenEventHeader();\r
 255     if (!genHeader) {\r
 256       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
 257       return;\r
 258     }\r
 259     genHeader->PrimaryVertex(vtxMC);\r
 260 \r
 261     // multipliticy of all MC primary tracks\r
 262     // in Zv, pt and eta ranges)\r
 263     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
 264 \r
 265   } // end bUseMC\r
 266 \r
 267   // get reconstructed vertex  \r
 268   //const AliESDVertex* vtxESD = 0; \r
 269   const AliESDVertex* vtxESD = 0; \r
 270   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
 271         vtxESD = esdEvent->GetPrimaryVertexTPC();\r
 272   }\r
 273   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
 274      vtxESD = esdEvent->GetPrimaryVertexTracks();\r
 275   }\r
 276   else {\r
 277         return;\r
 278   }\r
 279 \r
 280   if(!vtxESD) return;\r
 281 \r
 282   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
 283   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
 284   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
 285 \r
 286   // check event cuts\r
 287   if(isEventOK && isEventTriggered)\r
 288   {\r
 289     TRandom random;\r
 290 \r
 291     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
 292     {\r
 293       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
 294       if(!track) continue;\r
 295       if(track->Charge()==0) continue;\r
 296       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
 297       if(!accCuts->AcceptTrack(track)) continue;\r
 298 \r
 299       // downscale low-pT tracks\r
 300       if(TMath::Exp(2*track->Pt())<1000*random.Rndm()) continue;\r
 301 \r
 302       // Dump to the tree \r
 303       // vertex\r
 304       // TPC constrained tracks\r
 305       // InnerParams constrained tracks\r
 306       // TPC-ITS tracks\r
 307       // ITSout-InnerParams tracks\r
 308       // chi2 distance between TPC constrained and TPC-ITS tracks\r
 309       // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
 310       // chi2 distance between ITSout and InnerParams tracks\r
 311       // MC information\r
 312       \r
 313       Double_t x[3]; track->GetXYZ(x);\r
 314       Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
 315       Bool_t isOK = kFALSE;\r
 316 \r
 317       //\r
 318       // Constrain TPC-only tracks (TPCinner) to vertex\r
 319       //\r
 320       AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
 321       if (!tpcInner) continue;\r
 322       // transform to the track reference frame \r
 323       isOK = tpcInner->Rotate(track->GetAlpha());\r
 324       isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 325       if(!isOK) continue;\r
 326 \r
 327       // clone TPCinner has to be deleted\r
 328       AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
 329       if (!tpcInnerC) continue;\r
 330  \r
 331       // constrain TPCinner \r
 332       //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());\r
 333       isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
 334 \r
 335       // transform to the track reference frame \r
 336       isOK = tpcInnerC->Rotate(track->GetAlpha());\r
 337       isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 338 \r
 339       if(!isOK) {\r
 340         if(tpcInnerC) delete tpcInnerC;\r
 341         continue;\r
 342       }\r
 343 \r
 344 \r
 345       //\r
 346       // Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
 347       // to vertex\r
 348       //\r
 349       // clone track InnerParams has to be deleted\r
 350       AliExternalTrackParam * trackInnerC =  new AliExternalTrackParam(*(track->GetInnerParam()));\r
 351       if (!trackInnerC) continue;\r
 352  \r
 353       // constrain track InnerParams \r
 354       isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
 355 \r
 356       // transform to the track reference frame \r
 357       isOK = trackInnerC->Rotate(track->GetAlpha());\r
 358       isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 359 \r
 360       if(!isOK) {\r
 361         if(trackInnerC) delete trackInnerC;\r
 362         continue;\r
 363       }\r
 364       \r
 365       //\r
 366       // calculate chi2 between vi and vj vectors\r
 367       // with covi and covj covariance matrices\r
 368       // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
 369       //\r
 370       TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
 371       TMatrixD delta(1,5),  deltatrackC(1,5);\r
 372       TMatrixD covarM(5,5), covarMtrackC(5,5);\r
 373 \r
 374       for (Int_t ipar=0; ipar<5; ipar++) {\r
 375         deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 376         delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 377 \r
 378         deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 379         deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 380 \r
 381         for (Int_t jpar=0; jpar<5; jpar++) {\r
 382           Int_t index=track->GetIndex(ipar,jpar);\r
 383           covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
 384           covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
 385         }\r
 386       }\r
 387       // chi2 distance TPC constrained and TPC+ITS\r
 388       TMatrixD covarMInv = covarM.Invert();\r
 389       TMatrixD mat2 = covarMInv*deltaT;\r
 390       TMatrixD chi2 = delta*mat2; \r
 391       //chi2.Print();\r
 392 \r
 393       // chi2 distance TPC refitted constrained and TPC+ITS\r
 394       TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
 395       TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
 396       TMatrixD chi2trackC = deltatrackC*mat2trackC; \r
 397       //chi2trackC.Print();\r
 398 \r
 399 \r
 400       //\r
 401       // Propagate ITSout to TPC inner wall \r
 402       // and calculate chi2 distance to track (InnerParams)\r
 403       //\r
 404       const Double_t kTPCRadius=85; \r
 405       const Double_t kStep=3; \r
 406 \r
 407       // clone track InnerParams has to be deleted\r
 408       AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
 409       if (!trackInnerC2) continue;\r
 410       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
 411       {\r
 412           if(trackInnerC2) delete trackInnerC2;\r
 413           continue;\r
 414       }\r
 415 \r
 416       AliExternalTrackParam *outerITSc = new AliExternalTrackParam();\r
 417       if(!outerITSc) continue;\r
 418 \r
 419       TMatrixD chi2OuterITS(1,1);\r
 420       chi2OuterITS(0,0) = 0;\r
 421 \r
 422       if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
 423       {\r
 424         // propagate ITSout to TPC inner wall\r
 425         AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
 426 \r
 427         if(friendTrack) \r
 428         {\r
 429           if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) ) \r
 430           {\r
 431             if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
 432             {\r
 433               // transform ITSout to the track InnerParams reference frame \r
 434               isOK = kFALSE;\r
 435               isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
 436               isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
 437 \r
 438               if(!isOK) {\r
 439                 if(outerITSc) delete outerITSc;\r
 440                 if(trackInnerC2) delete trackInnerC2;\r
 441                 continue;\r
 442               }\r
 443               \r
 444               //\r
 445               // calculate chi2 between outerITS and innerParams\r
 446               // cov without z-coordinate at the moment\r
 447               //\r
 448               TMatrixD deltaTouterITS(4,1);\r
 449               TMatrixD deltaouterITS(1,4);\r
 450               TMatrixD covarMouterITS(4,4);\r
 451 \r
 452               Int_t kipar = 0;\r
 453               Int_t kjpar = 0;\r
 454               for (Int_t ipar=0; ipar<5; ipar++) {\r
 455                 if(ipar!=1) {\r
 456                   deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
 457                   deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
 458                 }\r
 459 \r
 460                 kjpar=0;\r
 461                 for (Int_t jpar=0; jpar<5; jpar++) {\r
 462                   Int_t index=outerITSc->GetIndex(ipar,jpar);\r
 463                   if(ipar !=1 || jpar!=1) {\r
 464                     covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
 465                   }\r
 466                   if(jpar!=1)  kjpar++;\r
 467                 }\r
 468                 if(ipar!=1) kipar++;\r
 469               }\r
 470 \r
 471               // chi2 distance ITSout and InnerParams\r
 472               TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
 473               TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
 474               chi2OuterITS = deltaouterITS*mat2outerITS; \r
 475               //chi2OuterITS.Print();\r
 476             } \r
 477           }\r
 478         }\r
 479       }\r
 480 \r
 481       //\r
 482       // MC info\r
 483       //\r
 484       TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
 485       TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
 486       Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
 487       Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
 488       Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
 489       Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
 490       Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
 491 \r
 492       AliTrackReference *refTPCIn = NULL;\r
 493       AliTrackReference *refITS = NULL;\r
 494       AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
 495       if(!trackInnerC3) continue;\r
 496 \r
 497       if(IsUseMCInfo()) \r
 498       {\r
 499         if(!stack) return;\r
 500 \r
 501         //\r
 502         // global track\r
 503         //\r
 504         Int_t label = TMath::Abs(track->GetLabel()); \r
 505         particle = stack->Particle(label);\r
 506         if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
 507         {\r
 508           particleMother = GetMother(particle,stack);\r
 509           mech = particle->GetUniqueID();\r
 510           isPrim = stack->IsPhysicalPrimary(label);\r
 511           isFromStrangess  = IsFromStrangeness(label,stack);\r
 512           isFromConversion = IsFromConversion(label,stack);\r
 513           isFromMaterial   = IsFromMaterial(label,stack);\r
 514         }\r
 515 \r
 516         //\r
 517         // TPC track\r
 518         //\r
 519         Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
 520         particleTPC = stack->Particle(labelTPC);\r
 521         if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
 522         {\r
 523           particleMotherTPC = GetMother(particleTPC,stack);\r
 524           mechTPC = particleTPC->GetUniqueID();\r
 525           isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
 526           isFromStrangessTPC  = IsFromStrangeness(labelTPC,stack);\r
 527           isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
 528           isFromMaterialTPC   = IsFromMaterial(labelTPC,stack);\r
 529         }\r
 530 \r
 531         //\r
 532         // store first track reference\r
 533         // for TPC track\r
 534         //\r
 535         TParticle *part=0;\r
 536         TClonesArray *trefs=0;\r
 537         Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
 538 \r
 539         if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
 540         {\r
 541           Int_t nTrackRef = trefs->GetEntries();\r
 542           //printf("nTrackRef %d \n",nTrackRef);\r
 543 \r
 544           Int_t countITS = 0;\r
 545           for (Int_t iref = 0; iref < nTrackRef; iref++) \r
 546           {\r
 547             AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
 548             //printf("ref %p \n",ref);\r
 549             //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
 550             //printf("AliTrackReference::kTPC  %d \n",AliTrackReference::kTPC);\r
 551 \r
 552 \r
 553              // Ref. in the middle ITS \r
 554             if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
 555             {\r
 556               if(!refITS && countITS==2) {\r
 557                 refITS = ref;\r
 558                 //printf("refITS %p \n",refITS);\r
 559               }\r
 560               countITS++;\r
 561             }\r
 562 \r
 563             // TPC\r
 564             if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
 565             {\r
 566               if(!refTPCIn) {\r
 567                 refTPCIn = ref;\r
 568                 //printf("refTPCIn %p \n",refTPCIn);\r
 569                 //break;\r
 570               }\r
 571             }\r
 572           }\r
 573 \r
 574           // transform inner params to TrackRef\r
 575           // reference frame\r
 576           if(refTPCIn && trackInnerC3) \r
 577           {\r
 578             Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
 579             isOK = kFALSE;\r
 580             isOK = trackInnerC3->Rotate(kRefPhi);\r
 581             isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
 582 \r
 583             if(!isOK){\r
 584               if(trackInnerC3) delete trackInnerC3;\r
 585               if(refTPCIn) delete refTPCIn;\r
 586             }\r
 587           }\r
 588         }\r
 589 \r
 590         //\r
 591         // ITS track\r
 592         //\r
 593         Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
 594         particleITS = stack->Particle(labelITS);\r
 595         if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
 596         {\r
 597           particleMotherITS = GetMother(particleITS,stack);\r
 598           mechITS = particleITS->GetUniqueID();\r
 599           isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
 600           isFromStrangessITS  = IsFromStrangeness(labelITS,stack);\r
 601           isFromConversionITS = IsFromConversion(labelITS,stack);\r
 602           isFromMaterialITS   = IsFromMaterial(labelITS,stack);\r
 603         }\r
 604       }\r
 605 \r
 606       //\r
 607       Double_t vert[3] = {0}; \r
 608       vert[0] = vtxESD->GetXv();\r
 609       vert[1] = vtxESD->GetYv();\r
 610       vert[2] = vtxESD->GetZv();\r
 611       Int_t mult = vtxESD->GetNContributors();\r
 612       Double_t bz = esdEvent->GetMagneticField();\r
 613       Double_t runNumber = esdEvent->GetRunNumber();\r
 614       Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
 615 \r
 616       //\r
 617       if(!fTreeSRedirector) return;\r
 618       (*fTreeSRedirector)<<"dNdPtTree"<<\r
 619         "runNumber="<<runNumber<<\r
 620         "evtTimeStamp="<<evtTimeStamp<<\r
 621         "Bz="<<bz<<\r
 622         "vertX="<<vert[0]<<\r
 623         "vertY="<<vert[1]<<\r
 624         "vertZ="<<vert[2]<<\r
 625         "mult="<<mult<<\r
 626         "esdTrack.="<<track<<\r
 627         "extTPCInnerC.="<<tpcInnerC<<\r
 628         "extInnerParamC.="<<trackInnerC<<\r
 629         "extInnerParam.="<<trackInnerC2<<\r
 630         "extOuterITS.="<<outerITSc<<\r
 631         "extInnerParamRef.="<<trackInnerC3<<\r
 632         "refTPCIn.="<<refTPCIn<<\r
 633         "refITS.="<<refITS<<\r
 634         "chi2TPCInnerC="<<chi2(0,0)<<\r
 635         "chi2InnerC="<<chi2trackC(0,0)<<\r
 636         "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
 637         "centralityF="<<centralityF<<\r
 638         "particle.="<<particle<<\r
 639         "particleMother.="<<particleMother<<\r
 640         "mech="<<mech<<\r
 641         "isPrim="<<isPrim<<\r
 642         "isFromStrangess="<<isFromStrangess<<\r
 643         "isFromConversion="<<isFromConversion<<\r
 644         "isFromMaterial="<<isFromMaterial<<\r
 645         "particleTPC.="<<particleTPC<<\r
 646         "particleMotherTPC.="<<particleMotherTPC<<\r
 647         "mechTPC="<<mechTPC<<\r
 648         "isPrimTPC="<<isPrimTPC<<\r
 649         "isFromStrangessTPC="<<isFromStrangessTPC<<\r
 650         "isFromConversionTPC="<<isFromConversionTPC<<\r
 651         "isFromMaterialTPC="<<isFromMaterialTPC<<\r
 652         "particleITS.="<<particleITS<<\r
 653         "particleMotherITS.="<<particleMotherITS<<\r
 654         "mechITS="<<mechITS<<\r
 655         "isPrimITS="<<isPrimITS<<\r
 656         "isFromStrangessITS="<<isFromStrangessITS<<\r
 657         "isFromConversionITS="<<isFromConversionITS<<\r
 658         "isFromMaterialITS="<<isFromMaterialITS<<\r
 659         "\n";\r
 660 \r
 661         if(tpcInnerC) delete tpcInnerC;\r
 662         if(trackInnerC) delete trackInnerC;\r
 663         if(trackInnerC2) delete trackInnerC2;\r
 664         if(outerITSc) delete outerITSc;\r
 665 \r
 666         if(trackInnerC3) delete trackInnerC3;\r
 667     }\r
 668   }\r
 669 \r
 670   PostData(0, fOutputSummary);\r
 671   //PostData(1, fOutput);\r
 672 }\r
 673 \r
 674 \r
 675 //_____________________________________________________________________________\r
 676 Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
 677 {\r
 678  // Constrain TPC inner params constrained\r
 679  //\r
 680       if(!tpcInnerC) return kFALSE; \r
 681       if(!vtx) return kFALSE;\r
 682 \r
 683       Double_t dz[2],cov[3];\r
 684       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
 685       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
 686       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
 687       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
 688 \r
 689 \r
 690       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
 691       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
 692       Double_t c[3]={covar[2],0.,covar[5]};\r
 693       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
 694       if (chi2C>kVeryBig) return kFALSE; \r
 695 \r
 696       if(!tpcInnerC->Update(p,c)) return kFALSE;\r
 697 \r
 698   return kTRUE;\r
 699 }\r
 700 \r
 701 //_____________________________________________________________________________\r
 702 Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
 703 {\r
 704  // Constrain TPC inner params constrained\r
 705  //\r
 706       if(!trackInnerC) return kFALSE; \r
 707       if(!vtx) return kFALSE;\r
 708 \r
 709       const Double_t kRadius  = 2.8; \r
 710       const Double_t kMaxStep = 1.0; \r
 711 \r
 712       Double_t dz[2],cov[3];\r
 713 \r
 714       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
 715       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
 716       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
 717 \r
 718       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
 719       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
 720 \r
 721       //\r
 722       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
 723       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
 724       Double_t c[3]={covar[2],0.,covar[5]};\r
 725       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
 726       if (chi2C>kVeryBig) return kFALSE; \r
 727 \r
 728       if(!trackInnerC->Update(p,c)) return kFALSE;\r
 729 \r
 730   return kTRUE;\r
 731 }\r
 732 \r
 733 \r
 734 //_____________________________________________________________________________\r
 735 TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
 736 {\r
 737   if(!particle) return NULL;\r
 738   if(!stack) return NULL;\r
 739 \r
 740   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
 741   TParticle* mother = NULL; \r
 742   mother = stack->Particle(motherLabel); \r
 743 \r
 744 return mother;\r
 745 }\r
 746 \r
 747 //_____________________________________________________________________________\r
 748 Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
 749 {\r
 750   Bool_t isFromConversion = kFALSE;\r
 751 \r
 752   if(stack) {\r
 753     TParticle* particle = stack->Particle(label);\r
 754 \r
 755     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
 756     {\r
 757        Int_t mech = particle->GetUniqueID(); // production mechanism \r
 758        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
 759 \r
 760        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
 761        TParticle* mother = stack->Particle(motherLabel); \r
 762        if(mother) {\r
 763           Int_t motherPdg = mother->GetPdgCode();\r
 764 \r
 765           if(!isPrim && mech==5 && motherPdg==kGamma) { \r
 766             isFromConversion=kTRUE; \r
 767           }\r
 768        }\r
 769     } \r
 770   } \r
 771 \r
 772   return isFromConversion;\r
 773 }\r
 774 \r
 775 //_____________________________________________________________________________\r
 776 Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
 777 {\r
 778   Bool_t isFromMaterial = kFALSE;\r
 779 \r
 780   if(stack) {\r
 781     TParticle* particle = stack->Particle(label);\r
 782 \r
 783     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
 784     {\r
 785        Int_t mech = particle->GetUniqueID(); // production mechanism \r
 786        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
 787 \r
 788        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
 789        TParticle* mother = stack->Particle(motherLabel); \r
 790        if(mother) {\r
 791           if(!isPrim && mech==13) { \r
 792             isFromMaterial=kTRUE; \r
 793           }\r
 794        }\r
 795      } \r
 796   } \r
 797 \r
 798   return isFromMaterial;\r
 799 }\r
 800 \r
 801 //_____________________________________________________________________________\r
 802 Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
 803 {\r
 804   Bool_t isFromStrangeness = kFALSE;\r
 805 \r
 806   if(stack) {\r
 807     TParticle* particle = stack->Particle(label);\r
 808 \r
 809     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
 810     {\r
 811        Int_t mech = particle->GetUniqueID(); // production mechanism \r
 812        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
 813 \r
 814        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
 815        TParticle* mother = stack->Particle(motherLabel); \r
 816        if(mother) {\r
 817           Int_t motherPdg = mother->GetPdgCode();\r
 818 \r
 819           // K+-, lambda, antilambda, K0s decays\r
 820           if(!isPrim && mech==4 && \r
 821               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
 822           {\r
 823             isFromStrangeness = kTRUE;\r
 824           } \r
 825        }\r
 826     } \r
 827   } \r
 828 \r
 829   return isFromStrangeness;\r
 830 }\r
 831 \r
 832 \r
 833 //_____________________________________________________________________________\r
 834 void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
 835 {\r
 836   //\r
 837   // Called one at the end \r
 838   // locally on working node\r
 839   //\r
 840 \r
 841   // Post output data.\r
 842   PostData(1, fOutput);\r
 843   //PostData(0, fOutputSummary);\r
 844 }\r
 845 \r
 846 //_____________________________________________________________________________\r
 847 void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
 848 {\r
 849   // Called one at the end \r
 850   \r
 851   // check output data\r
 852   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(0));\r
 853   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
 854   if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;\r
 855 \r
 856   TChain* chain = new TChain("dNdPtTree");\r
 857   if(!chain) return;\r
 858   chain->Add("dNdPtOutliersAnalysisPbPb.root");\r
 859   TTree *tree = chain->CopyTree("1");\r
 860   if (chain) { delete chain; chain=0; }\r
 861   if(!tree) return;\r
 862   tree->Print();\r
 863   fOutputSummary = tree;\r
 864 \r
 865   if (!fOutputSummary) {\r
 866     Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable GetOutputData(0)==0x0 ..." );\r
 867     return;\r
 868   }\r
 869   \r
 870 \r
 871 \r
 872   PostData(0, fOutputSummary);\r
 873   //PostData(1, fOutput);\r
 874 }\r