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 "TFile.h"\r
  28 #include "TMatrixD.h"\r
  29 #include "TRandom3.h"\r
  30 \r
  31 #include "AliHeader.h"  \r
  32 #include "AliGenEventHeader.h"  \r
  33 #include "AliInputEventHandler.h"  \r
  34 #include "AliStack.h"  \r
  35 #include "AliTrackReference.h"  \r
  36 \r
  37 #include "AliPhysicsSelection.h"\r
  38 #include "AliAnalysisTask.h"\r
  39 #include "AliAnalysisManager.h"\r
  40 #include "AliESDEvent.h"\r
  41 #include "AliESDfriend.h"\r
  42 #include "AliMCEvent.h"\r
  43 #include "AliESDInputHandler.h"\r
  44 #include "AliESDVertex.h"\r
  45 #include "AliTracker.h"\r
  46 #include "AliGeomManager.h"\r
  47 \r
  48 #include "AliCentrality.h"\r
  49 #include "AliESDVZERO.h"\r
  50 #include "AliMultiplicity.h"\r
  51 \r
  52 #include "AliESDtrackCuts.h"\r
  53 #include "AliMCEventHandler.h"\r
  54 #include "AlidNdPt.h"\r
  55 #include "AlidNdPtEventCuts.h"\r
  56 #include "AlidNdPtAcceptanceCuts.h"\r
  57 \r
  58 #include "AlidNdPtTrackDumpTask.h"\r
  59 #include "AliKFParticle.h"\r
  60 #include "AliESDv0.h"\r
  61 \r
  62 using namespace std;\r
  63 \r
  64 ClassImp(AlidNdPtTrackDumpTask)\r
  65 \r
  66 //_____________________________________________________________________________\r
  67 AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name) \r
  68   : AliAnalysisTaskSE(name)\r
  69   , fESD(0)\r
  70   , fMC(0)\r
  71   , fESDfriend(0)\r
  72   , fOutput(0)\r
  73   , fPitList(0)\r
  74   , fUseMCInfo(kFALSE)\r
  75   , fdNdPtEventCuts(0)\r
  76   , fdNdPtAcceptanceCuts(0)\r
  77   , fdNdPtRecAcceptanceCuts(0)\r
  78   , fEsdTrackCuts(0)\r
  79   , fTrigger(AliTriggerAnalysis::kMB1) \r
  80   , fAnalysisMode(AlidNdPtHelper::kTPC) \r
  81   , fOutputSummary(0)\r
  82   , fTreeSRedirector(0)\r
  83   , fCentralityEstimator(0)\r
  84   , fLowPtTrackDownscaligF(0)\r
  85   , fLowPtV0DownscaligF(0)\r
  86 {\r
  87   // Constructor\r
  88 \r
  89   // Define input and output slots here\r
  90   DefineOutput(0, TTree::Class());\r
  91   //DefineOutput(1, TList::Class());\r
  92 \r
  93 }\r
  94 \r
  95 //_____________________________________________________________________________\r
  96 AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()\r
  97 {\r
  98   if(fOutput) delete fOutput;  fOutput =0; \r
  99   //if(fOutputSummary) delete fOutputSummary;  fOutputSummary =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 output tree\r
 132   //\r
 133   fTreeSRedirector = new TTreeSRedirector("jotwinow_HighPt_TrackAndV0_Trees.root");\r
 134 \r
 135   PostData(0, fOutputSummary);\r
 136   //PostData(1, fOutput);\r
 137 }\r
 138 \r
 139 //_____________________________________________________________________________\r
 140 void AlidNdPtTrackDumpTask::UserExec(Option_t *) \r
 141 {\r
 142   //\r
 143   // Called for each event\r
 144   //\r
 145 \r
 146   // ESD event\r
 147   fESD = (AliESDEvent*) (InputEvent());\r
 148   if (!fESD) {\r
 149     Printf("ERROR: ESD event not available");\r
 150     return;\r
 151   }\r
 152 \r
 153   // MC event\r
 154   if(fUseMCInfo) {\r
 155     fMC = MCEvent();\r
 156     if (!fMC) {\r
 157       Printf("ERROR: MC event not available");\r
 158       return;\r
 159     }\r
 160   }\r
 161 \r
 162   fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
 163   if(!fESDfriend) {\r
 164     Printf("ERROR: ESD friends not available");\r
 165   }\r
 166 \r
 167   //\r
 168   Process(fESD,fMC,fESDfriend);\r
 169   ProcessV0(fESD,fMC,fESDfriend);\r
 170 \r
 171   // Post output data.\r
 172   PostData(0, fOutputSummary);\r
 173   //PostData(1, fOutput);\r
 174 }\r
 175 \r
 176 //_____________________________________________________________________________\r
 177 void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
 178 {\r
 179   //\r
 180   // Process real and/or simulated events\r
 181   //\r
 182   if(!esdEvent) {\r
 183     AliDebug(AliLog::kError, "esdEvent not available");\r
 184     return;\r
 185   }\r
 186 \r
 187   // get selection cuts\r
 188   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
 189   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
 190   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
 191 \r
 192   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
 193     AliDebug(AliLog::kError, "cuts not available");\r
 194     return;\r
 195   }\r
 196 \r
 197   // trigger selection\r
 198   Bool_t isEventTriggered = kTRUE;\r
 199   AliPhysicsSelection *physicsSelection = NULL;\r
 200   AliTriggerAnalysis* triggerAnalysis = NULL;\r
 201 \r
 202   // \r
 203   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
 204   if (!inputHandler)\r
 205   {\r
 206     Printf("ERROR: Could not receive input handler");\r
 207     return;\r
 208   }\r
 209    \r
 210   // get file name\r
 211   TTree *chain = (TChain*)GetInputData(0);\r
 212   if(!chain) { \r
 213     Printf("ERROR: Could not receive input chain");\r
 214     return;\r
 215   }\r
 216   TObjString fileName(chain->GetCurrentFile()->GetName());\r
 217 \r
 218   // trigger\r
 219   if(evtCuts->IsTriggerRequired())  \r
 220   {\r
 221     // always MB\r
 222     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
 223 \r
 224     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
 225     if(!physicsSelection) return;\r
 226     //SetPhysicsTriggerSelection(physicsSelection);\r
 227 \r
 228     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
 229       // set trigger (V0AND)\r
 230       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
 231       if(!triggerAnalysis) return;\r
 232       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
 233     }\r
 234   }\r
 235 \r
 236   // centrality determination\r
 237   Float_t centralityF = -1;\r
 238   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
 239   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 240 \r
 241   // use MC information\r
 242   AliHeader* header = 0;\r
 243   AliGenEventHeader* genHeader = 0;\r
 244   AliStack* stack = 0;\r
 245   TArrayF vtxMC(3);\r
 246 \r
 247   Int_t multMCTrueTracks = 0;\r
 248   if(IsUseMCInfo())\r
 249   {\r
 250     //\r
 251     if(!mcEvent) {\r
 252       AliDebug(AliLog::kError, "mcEvent not available");\r
 253       return;\r
 254     }\r
 255     // get MC event header\r
 256     header = mcEvent->Header();\r
 257     if (!header) {\r
 258       AliDebug(AliLog::kError, "Header not available");\r
 259       return;\r
 260     }\r
 261     // MC particle stack\r
 262     stack = mcEvent->Stack();\r
 263     if (!stack) {\r
 264       AliDebug(AliLog::kError, "Stack not available");\r
 265       return;\r
 266     }\r
 267 \r
 268     // get MC vertex\r
 269     genHeader = header->GenEventHeader();\r
 270     if (!genHeader) {\r
 271       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
 272       return;\r
 273     }\r
 274     genHeader->PrimaryVertex(vtxMC);\r
 275 \r
 276     // multipliticy of all MC primary tracks\r
 277     // in Zv, pt and eta ranges)\r
 278     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
 279 \r
 280   } // end bUseMC\r
 281 \r
 282   // get reconstructed vertex  \r
 283   //const AliESDVertex* vtxESD = 0; \r
 284   const AliESDVertex* vtxESD = 0; \r
 285   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
 286         vtxESD = esdEvent->GetPrimaryVertexTPC();\r
 287   }\r
 288   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
 289      vtxESD = esdEvent->GetPrimaryVertexTracks();\r
 290   }\r
 291   else {\r
 292         return;\r
 293   }\r
 294 \r
 295   if(!vtxESD) return;\r
 296 \r
 297   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
 298   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
 299   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
 300 \r
 301   // check event cuts\r
 302   if(isEventOK && isEventTriggered)\r
 303   {\r
 304     TRandom3 random;\r
 305 \r
 306     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
 307     {\r
 308       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
 309       if(!track) continue;\r
 310       if(track->Charge()==0) continue;\r
 311       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
 312       if(!accCuts->AcceptTrack(track)) continue;\r
 313 \r
 314       // downscale low-pT tracks\r
 315       Double_t scalempt= TMath::Min(track->Pt(),10.);\r
 316       if(TMath::Exp(2*scalempt)<fLowPtTrackDownscaligF*random.Rndm()) continue;\r
 317 \r
 318       // Dump to the tree \r
 319       // vertex\r
 320       // TPC constrained tracks\r
 321       // InnerParams constrained tracks\r
 322       // TPC-ITS tracks\r
 323       // ITSout-InnerParams tracks\r
 324       // chi2 distance between TPC constrained and TPC-ITS tracks\r
 325       // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
 326       // chi2 distance between ITSout and InnerParams tracks\r
 327       // MC information\r
 328       \r
 329       Double_t x[3]; track->GetXYZ(x);\r
 330       Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
 331       Bool_t isOK = kFALSE;\r
 332 \r
 333       //\r
 334       // Constrain TPC-only tracks (TPCinner) to vertex\r
 335       //\r
 336       AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
 337       if (!tpcInner) continue;\r
 338       // transform to the track reference frame \r
 339       isOK = tpcInner->Rotate(track->GetAlpha());\r
 340       isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 341       if(!isOK) continue;\r
 342 \r
 343       // clone TPCinner has to be deleted\r
 344       AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
 345       if (!tpcInnerC) continue;\r
 346  \r
 347       // constrain TPCinner \r
 348       //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());\r
 349       isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
 350 \r
 351       // transform to the track reference frame \r
 352       isOK = tpcInnerC->Rotate(track->GetAlpha());\r
 353       isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 354 \r
 355       if(!isOK) {\r
 356         if(tpcInnerC) delete tpcInnerC;\r
 357         continue;\r
 358       }\r
 359 \r
 360 \r
 361       //\r
 362       // Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
 363       // to vertex\r
 364       //\r
 365       // clone track InnerParams has to be deleted\r
 366       AliExternalTrackParam * trackInnerC =  new AliExternalTrackParam(*(track->GetInnerParam()));\r
 367       if (!trackInnerC) continue;\r
 368  \r
 369       // constrain track InnerParams \r
 370       isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
 371 \r
 372       // transform to the track reference frame \r
 373       isOK = trackInnerC->Rotate(track->GetAlpha());\r
 374       isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 375 \r
 376       if(!isOK) {\r
 377         if(trackInnerC) delete trackInnerC;\r
 378         continue;\r
 379       }\r
 380       \r
 381       //\r
 382       // calculate chi2 between vi and vj vectors\r
 383       // with covi and covj covariance matrices\r
 384       // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
 385       //\r
 386       TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
 387       TMatrixD delta(1,5),  deltatrackC(1,5);\r
 388       TMatrixD covarM(5,5), covarMtrackC(5,5);\r
 389 \r
 390       for (Int_t ipar=0; ipar<5; ipar++) {\r
 391         deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 392         delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 393 \r
 394         deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 395         deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 396 \r
 397         for (Int_t jpar=0; jpar<5; jpar++) {\r
 398           Int_t index=track->GetIndex(ipar,jpar);\r
 399           covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
 400           covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
 401         }\r
 402       }\r
 403       // chi2 distance TPC constrained and TPC+ITS\r
 404       TMatrixD covarMInv = covarM.Invert();\r
 405       TMatrixD mat2 = covarMInv*deltaT;\r
 406       TMatrixD chi2 = delta*mat2; \r
 407       //chi2.Print();\r
 408 \r
 409       // chi2 distance TPC refitted constrained and TPC+ITS\r
 410       TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
 411       TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
 412       TMatrixD chi2trackC = deltatrackC*mat2trackC; \r
 413       //chi2trackC.Print();\r
 414 \r
 415 \r
 416       //\r
 417       // Propagate ITSout to TPC inner wall \r
 418       // and calculate chi2 distance to track (InnerParams)\r
 419       //\r
 420       const Double_t kTPCRadius=85; \r
 421       const Double_t kStep=3; \r
 422 \r
 423       // clone track InnerParams has to be deleted\r
 424       AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
 425       if (!trackInnerC2) continue;\r
 426       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
 427       {\r
 428           if(trackInnerC2) delete trackInnerC2;\r
 429           continue;\r
 430       }\r
 431 \r
 432       AliExternalTrackParam *outerITSc = new AliExternalTrackParam();\r
 433       if(!outerITSc) continue;\r
 434 \r
 435       TMatrixD chi2OuterITS(1,1);\r
 436       chi2OuterITS(0,0) = 0;\r
 437 \r
 438       if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
 439       {\r
 440         // propagate ITSout to TPC inner wall\r
 441         AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
 442 \r
 443         if(friendTrack) \r
 444         {\r
 445           if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) ) \r
 446           {\r
 447             if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
 448             {\r
 449               // transform ITSout to the track InnerParams reference frame \r
 450               isOK = kFALSE;\r
 451               isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
 452               isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
 453 \r
 454               if(!isOK) {\r
 455                 if(outerITSc) delete outerITSc;\r
 456                 if(trackInnerC2) delete trackInnerC2;\r
 457                 continue;\r
 458               }\r
 459               \r
 460               //\r
 461               // calculate chi2 between outerITS and innerParams\r
 462               // cov without z-coordinate at the moment\r
 463               //\r
 464               TMatrixD deltaTouterITS(4,1);\r
 465               TMatrixD deltaouterITS(1,4);\r
 466               TMatrixD covarMouterITS(4,4);\r
 467 \r
 468               Int_t kipar = 0;\r
 469               Int_t kjpar = 0;\r
 470               for (Int_t ipar=0; ipar<5; ipar++) {\r
 471                 if(ipar!=1) {\r
 472                   deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
 473                   deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
 474                 }\r
 475 \r
 476                 kjpar=0;\r
 477                 for (Int_t jpar=0; jpar<5; jpar++) {\r
 478                   Int_t index=outerITSc->GetIndex(ipar,jpar);\r
 479                   if(ipar !=1 || jpar!=1) {\r
 480                     covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
 481                   }\r
 482                   if(jpar!=1)  kjpar++;\r
 483                 }\r
 484                 if(ipar!=1) kipar++;\r
 485               }\r
 486 \r
 487               // chi2 distance ITSout and InnerParams\r
 488               TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
 489               TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
 490               chi2OuterITS = deltaouterITS*mat2outerITS; \r
 491               //chi2OuterITS.Print();\r
 492             } \r
 493           }\r
 494         }\r
 495       }\r
 496 \r
 497       //\r
 498       // MC info\r
 499       //\r
 500       TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
 501       TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
 502       Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
 503       Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
 504       Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
 505       Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
 506       Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
 507 \r
 508       AliTrackReference *refTPCIn = NULL;\r
 509       AliTrackReference *refITS = NULL;\r
 510       AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
 511       if(!trackInnerC3) continue;\r
 512 \r
 513       if(IsUseMCInfo()) \r
 514       {\r
 515         if(!stack) return;\r
 516 \r
 517         //\r
 518         // global track\r
 519         //\r
 520         Int_t label = TMath::Abs(track->GetLabel()); \r
 521         particle = stack->Particle(label);\r
 522         if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
 523         {\r
 524           particleMother = GetMother(particle,stack);\r
 525           mech = particle->GetUniqueID();\r
 526           isPrim = stack->IsPhysicalPrimary(label);\r
 527           isFromStrangess  = IsFromStrangeness(label,stack);\r
 528           isFromConversion = IsFromConversion(label,stack);\r
 529           isFromMaterial   = IsFromMaterial(label,stack);\r
 530         }\r
 531 \r
 532         //\r
 533         // TPC track\r
 534         //\r
 535         Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
 536         particleTPC = stack->Particle(labelTPC);\r
 537         if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
 538         {\r
 539           particleMotherTPC = GetMother(particleTPC,stack);\r
 540           mechTPC = particleTPC->GetUniqueID();\r
 541           isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
 542           isFromStrangessTPC  = IsFromStrangeness(labelTPC,stack);\r
 543           isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
 544           isFromMaterialTPC   = IsFromMaterial(labelTPC,stack);\r
 545         }\r
 546 \r
 547         //\r
 548         // store first track reference\r
 549         // for TPC track\r
 550         //\r
 551         TParticle *part=0;\r
 552         TClonesArray *trefs=0;\r
 553         Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
 554 \r
 555         if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
 556         {\r
 557           Int_t nTrackRef = trefs->GetEntries();\r
 558           //printf("nTrackRef %d \n",nTrackRef);\r
 559 \r
 560           Int_t countITS = 0;\r
 561           for (Int_t iref = 0; iref < nTrackRef; iref++) \r
 562           {\r
 563             AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
 564             //printf("ref %p \n",ref);\r
 565             //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
 566             //printf("AliTrackReference::kTPC  %d \n",AliTrackReference::kTPC);\r
 567 \r
 568 \r
 569              // Ref. in the middle ITS \r
 570             if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
 571             {\r
 572               if(!refITS && countITS==2) {\r
 573                 refITS = ref;\r
 574                 //printf("refITS %p \n",refITS);\r
 575               }\r
 576               countITS++;\r
 577             }\r
 578 \r
 579             // TPC\r
 580             if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
 581             {\r
 582               if(!refTPCIn) {\r
 583                 refTPCIn = ref;\r
 584                 //printf("refTPCIn %p \n",refTPCIn);\r
 585                 //break;\r
 586               }\r
 587             }\r
 588           }\r
 589 \r
 590           // transform inner params to TrackRef\r
 591           // reference frame\r
 592           if(refTPCIn && trackInnerC3) \r
 593           {\r
 594             Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
 595             isOK = kFALSE;\r
 596             isOK = trackInnerC3->Rotate(kRefPhi);\r
 597             isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
 598 \r
 599             if(!isOK){\r
 600               if(trackInnerC3) delete trackInnerC3;\r
 601               if(refTPCIn) delete refTPCIn;\r
 602             }\r
 603           }\r
 604         }\r
 605 \r
 606         //\r
 607         // ITS track\r
 608         //\r
 609         Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
 610         particleITS = stack->Particle(labelITS);\r
 611         if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
 612         {\r
 613           particleMotherITS = GetMother(particleITS,stack);\r
 614           mechITS = particleITS->GetUniqueID();\r
 615           isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
 616           isFromStrangessITS  = IsFromStrangeness(labelITS,stack);\r
 617           isFromConversionITS = IsFromConversion(labelITS,stack);\r
 618           isFromMaterialITS   = IsFromMaterial(labelITS,stack);\r
 619         }\r
 620       }\r
 621 \r
 622       //\r
 623       Double_t vert[3] = {0}; \r
 624       vert[0] = vtxESD->GetXv();\r
 625       vert[1] = vtxESD->GetYv();\r
 626       vert[2] = vtxESD->GetZv();\r
 627       Int_t mult = vtxESD->GetNContributors();\r
 628       Double_t bz = esdEvent->GetMagneticField();\r
 629       Double_t runNumber = esdEvent->GetRunNumber();\r
 630       Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
 631       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
 632 \r
 633 \r
 634       //\r
 635       if(!fTreeSRedirector) return;\r
 636       (*fTreeSRedirector)<<"dNdPtTree"<<\r
 637         "fileName.="<<&fileName<<\r
 638         "runNumber="<<runNumber<<\r
 639         "evtTimeStamp="<<evtTimeStamp<<\r
 640         "evtNumberInFile="<<evtNumberInFile<<\r
 641         "Bz="<<bz<<\r
 642         "vertX="<<vert[0]<<\r
 643         "vertY="<<vert[1]<<\r
 644         "vertZ="<<vert[2]<<\r
 645         "mult="<<mult<<\r
 646         "esdTrack.="<<track<<\r
 647         "extTPCInnerC.="<<tpcInnerC<<\r
 648         "extInnerParamC.="<<trackInnerC<<\r
 649         "extInnerParam.="<<trackInnerC2<<\r
 650         "extOuterITS.="<<outerITSc<<\r
 651         "extInnerParamRef.="<<trackInnerC3<<\r
 652         "refTPCIn.="<<refTPCIn<<\r
 653         "refITS.="<<refITS<<\r
 654         "chi2TPCInnerC="<<chi2(0,0)<<\r
 655         "chi2InnerC="<<chi2trackC(0,0)<<\r
 656         "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
 657         "centralityF="<<centralityF<<\r
 658         "particle.="<<particle<<\r
 659         "particleMother.="<<particleMother<<\r
 660         "mech="<<mech<<\r
 661         "isPrim="<<isPrim<<\r
 662         "isFromStrangess="<<isFromStrangess<<\r
 663         "isFromConversion="<<isFromConversion<<\r
 664         "isFromMaterial="<<isFromMaterial<<\r
 665         "particleTPC.="<<particleTPC<<\r
 666         "particleMotherTPC.="<<particleMotherTPC<<\r
 667         "mechTPC="<<mechTPC<<\r
 668         "isPrimTPC="<<isPrimTPC<<\r
 669         "isFromStrangessTPC="<<isFromStrangessTPC<<\r
 670         "isFromConversionTPC="<<isFromConversionTPC<<\r
 671         "isFromMaterialTPC="<<isFromMaterialTPC<<\r
 672         "particleITS.="<<particleITS<<\r
 673         "particleMotherITS.="<<particleMotherITS<<\r
 674         "mechITS="<<mechITS<<\r
 675         "isPrimITS="<<isPrimITS<<\r
 676         "isFromStrangessITS="<<isFromStrangessITS<<\r
 677         "isFromConversionITS="<<isFromConversionITS<<\r
 678         "isFromMaterialITS="<<isFromMaterialITS<<\r
 679         "\n";\r
 680 \r
 681         if(tpcInnerC) delete tpcInnerC;\r
 682         if(trackInnerC) delete trackInnerC;\r
 683         if(trackInnerC2) delete trackInnerC2;\r
 684         if(outerITSc) delete outerITSc;\r
 685 \r
 686         if(trackInnerC3) delete trackInnerC3;\r
 687     }\r
 688   }\r
 689 \r
 690   PostData(0, fOutputSummary);\r
 691   //PostData(1, fOutput);\r
 692 }\r
 693 \r
 694 //_____________________________________________________________________________\r
 695 void AlidNdPtTrackDumpTask::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
 696 {\r
 697   //\r
 698   // Process real and/or simulated events\r
 699   //\r
 700   if(!esdEvent) {\r
 701     AliDebug(AliLog::kError, "esdEvent not available");\r
 702     return;\r
 703   }\r
 704 \r
 705   // get selection cuts\r
 706   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
 707   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
 708   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
 709 \r
 710   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
 711     AliDebug(AliLog::kError, "cuts not available");\r
 712     return;\r
 713   }\r
 714 \r
 715 \r
 716 \r
 717 \r
 718   // trigger selection\r
 719   Bool_t isEventTriggered = kTRUE;\r
 720   AliPhysicsSelection *physicsSelection = NULL;\r
 721   AliTriggerAnalysis* triggerAnalysis = NULL;\r
 722 \r
 723   // \r
 724   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
 725   if (!inputHandler)\r
 726   {\r
 727     Printf("ERROR: Could not receive input handler");\r
 728     return;\r
 729   }\r
 730    \r
 731   // get file name\r
 732   TTree *chain = (TChain*)GetInputData(0);\r
 733   if(!chain) { \r
 734     Printf("ERROR: Could not receive input chain");\r
 735     return;\r
 736   }\r
 737   TObjString fileName(chain->GetCurrentFile()->GetName());\r
 738 \r
 739   // trigger\r
 740   if(evtCuts->IsTriggerRequired())  \r
 741   {\r
 742     // always MB\r
 743     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
 744 \r
 745     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
 746     if(!physicsSelection) return;\r
 747     //SetPhysicsTriggerSelection(physicsSelection);\r
 748 \r
 749     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
 750       // set trigger (V0AND)\r
 751       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
 752       if(!triggerAnalysis) return;\r
 753       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
 754     }\r
 755   }\r
 756 \r
 757   // centrality determination\r
 758   Float_t centralityF = -1;\r
 759   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
 760   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 761 \r
 762 \r
 763   // get reconstructed vertex  \r
 764   //const AliESDVertex* vtxESD = 0; \r
 765   const AliESDVertex* vtxESD = 0; \r
 766   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
 767         vtxESD = esdEvent->GetPrimaryVertexTPC();\r
 768   }\r
 769   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
 770      vtxESD = esdEvent->GetPrimaryVertexTracks();\r
 771   }\r
 772   else {\r
 773         return;\r
 774   }\r
 775 \r
 776   if(!vtxESD) return;\r
 777 \r
 778   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
 779   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
 780   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
 781 \r
 782   // check event cuts\r
 783   if(isEventOK && isEventTriggered) {\r
 784   //\r
 785   // Dump the pt downscaled V0 into the tree\r
 786   // \r
 787   //\r
 788   Int_t ntracks = esdEvent->GetNumberOfTracks();\r
 789   Int_t nV0s = esdEvent->GetNumberOfV0s();\r
 790   Int_t run = esdEvent->GetRunNumber();\r
 791   Int_t time= esdEvent->GetTimeStamp();\r
 792   Int_t evNr=esdEvent->GetEventNumberInFile();\r
 793   \r
 794   for (Int_t iv0=0; iv0<nV0s; iv0++){\r
 795     AliESDv0 * v0 = esdEvent->GetV0(iv0);\r
 796     if (!v0) continue;\r
 797     AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));\r
 798     AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));\r
 799     if (!track0) continue;\r
 800     if (!track1) continue;\r
 801     if (track0->GetSign()<0) {\r
 802       track1 = esdEvent->GetTrack(v0->GetIndex(0));\r
 803       track0 = esdEvent->GetTrack(v0->GetIndex(1));\r
 804     }\r
 805     //\r
 806     Bool_t isDownscaled = IsV0Downscaled(v0);\r
 807     if (isDownscaled) continue;\r
 808     AliKFParticle kfparticle; //\r
 809     Int_t type=GetKFParticle(v0,esdEvent,kfparticle);\r
 810     if (type==0) continue;   \r
 811 \r
 812     if(!fTreeSRedirector) return;\r
 813     (*fTreeSRedirector)<<"V0s"<<\r
 814       "isDownscaled="<<isDownscaled<<\r
 815       "run="<<run<<\r
 816       "fname="<<&fileName<<\r
 817       "time="<<time<<\r
 818       "evNr="<<evNr<<\r
 819       "type="<<type<<\r
 820       "ntracks="<<ntracks<<\r
 821       "v0.="<<v0<<\r
 822       "kf.="<<&kfparticle<<\r
 823       "track0.="<<track0<<\r
 824       "track1.="<<track1<<\r
 825       "centralityF="<<centralityF<<\r
 826       "\n";\r
 827   }\r
 828   }\r
 829   PostData(0, fOutputSummary);\r
 830 }\r
 831 \r
 832 //_____________________________________________________________________________\r
 833 Int_t   AlidNdPtTrackDumpTask::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)\r
 834 {\r
 835   //\r
 836   // Create KF particle in case the V0 fullfill selection criteria\r
 837   //\r
 838   // Selection criteria\r
 839   //  0. algorithm cut\r
 840   //  1. track cut\r
 841   //  3. chi2 cut\r
 842   //  4. rough mass cut\r
 843   //  5. Normalized pointing angle cut\r
 844   //\r
 845   const Double_t cutMass=0.2;\r
 846   const Double_t kSigmaDCACut=3;\r
 847   //\r
 848   // 0.) algo cut - accept only on the fly\r
 849   //\r
 850   if (v0->GetOnFlyStatus() ==kFALSE) return 0;     \r
 851   //\r
 852   // 1.) track cut\r
 853   // \r
 854   AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));\r
 855   AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));\r
 856   /*\r
 857     TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";\r
 858     TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";\r
 859     TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";\r
 860   */  \r
 861   if (TMath::Abs(track0->GetTgl())>1) return 0;\r
 862   if (TMath::Abs(track1->GetTgl())>1) return 0;\r
 863   if ((track0->GetTPCClusterInfo(2,1))<100) return 0;\r
 864   if ((track1->GetTPCClusterInfo(2,1))<100) return 0;\r
 865   //if ((track0->GetITSclusters(0))<2) return 0;\r
 866   //if ((track1->GetITSclusters(0))<2) return 0; \r
 867   Float_t pos0[2]={0}, cov0[3]={0};\r
 868   Float_t pos1[2]={0}, cov1[3]={0};\r
 869   track0->GetImpactParameters(pos0,cov0);\r
 870   track0->GetImpactParameters(pos1,cov1);\r
 871   //\r
 872   if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;\r
 873   if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;\r
 874   // \r
 875   //\r
 876   // 3.) Chi2 cut\r
 877   //\r
 878   Double_t chi2KF = v0->GetKFInfo(2,2,2);\r
 879   if (chi2KF>25) return 0;\r
 880   //\r
 881   // 4.) Rough mass cut - 0.200 GeV\r
 882   //\r
 883   static Double_t masses[2]={-1};\r
 884   if (masses[0]<0){\r
 885     masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();\r
 886     masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();\r
 887   }\r
 888   Double_t mass00=  v0->GetEffMass(0,0);\r
 889   Double_t mass22=  v0->GetEffMass(2,2);\r
 890   Double_t mass42=  v0->GetEffMass(4,2);\r
 891   Double_t mass24=  v0->GetEffMass(2,4);\r
 892   Bool_t massOK=kFALSE;\r
 893   Int_t type=0;\r
 894   Int_t ptype=0;\r
 895   Double_t dmass=1;\r
 896   Int_t p1=0, p2=0;\r
 897   if (TMath::Abs(mass00-0)<cutMass) {\r
 898     massOK=kTRUE; type+=1; \r
 899     if (TMath::Abs(mass00-0)<dmass) {\r
 900       ptype=1;\r
 901       dmass=TMath::Abs(mass00-0);      \r
 902       p1=0; p2=0;\r
 903     } \r
 904   }\r
 905   if (TMath::Abs(mass24-masses[1])<cutMass) {\r
 906     massOK=kTRUE; type+=2; \r
 907     if (TMath::Abs(mass24-masses[1])<dmass){\r
 908       dmass = TMath::Abs(mass24-masses[1]);\r
 909       ptype=2;\r
 910       p1=2; p2=4;\r
 911     }\r
 912   }\r
 913   if (TMath::Abs(mass42-masses[1])<cutMass) {\r
 914     massOK=kTRUE; type+=4;\r
 915     if (TMath::Abs(mass42-masses[1])<dmass){\r
 916       dmass = TMath::Abs(mass42-masses[1]);\r
 917       ptype=4;\r
 918       p1=4; p2=2;\r
 919     }\r
 920   }\r
 921   if (TMath::Abs(mass22-masses[0])<cutMass) {\r
 922     massOK=kTRUE; type+=8;\r
 923     if (TMath::Abs(mass22-masses[0])<dmass){\r
 924       dmass = TMath::Abs(mass22-masses[0]);\r
 925       ptype=8;\r
 926       p1=2; p2=2;\r
 927     }\r
 928   }\r
 929   if (type==0) return 0;\r
 930   //\r
 931   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};\r
 932   const AliExternalTrackParam *paramP = v0->GetParamP();\r
 933   const AliExternalTrackParam *paramN = v0->GetParamN();\r
 934   if (paramP->GetSign()<0){\r
 935     paramP=v0->GetParamP();\r
 936     paramN=v0->GetParamN();\r
 937   }\r
 938   //Double_t *pparam1 = (Double_t*)paramP->GetParameter();\r
 939   //Double_t *pparam2 = (Double_t*)paramN->GetParameter();\r
 940   //\r
 941   AliKFParticle kfp1( *paramP, spdg[p1]  );\r
 942   AliKFParticle kfp2( *paramN, -1 *spdg[p2]  );\r
 943   AliKFParticle V0KF;\r
 944   (V0KF)+=kfp1;\r
 945   (V0KF)+=kfp2;\r
 946   kfparticle=V0KF;\r
 947   //\r
 948   // Pointing angle\r
 949   //\r
 950   Double_t  errPhi    = V0KF.GetErrPhi();\r
 951   Double_t  pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());\r
 952   if (pointAngle/errPhi>10) return 0;  \r
 953   //\r
 954   return ptype;  \r
 955 }\r
 956 \r
 957 //_____________________________________________________________________________\r
 958 Bool_t AlidNdPtTrackDumpTask::IsV0Downscaled(AliESDv0 *const v0)\r
 959 {\r
 960   //\r
 961   // Downscale randomly low pt V0\r
 962   //\r
 963   //return kFALSE;\r
 964   Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());\r
 965   Double_t scalempt= TMath::Min(maxPt,10.);\r
 966   if (TMath::Exp(2*scalempt)<fLowPtV0DownscaligF*gRandom->Rndm()) return kTRUE;\r
 967   return kFALSE;\r
 968   /*\r
 969   \r
 970     TH1F his1("his1","his1",100,0,10);\r
 971     TH1F his2("his2","his2",100,0,10);\r
 972     {for (Int_t i=0; i<10000; i++){\r
 973        Double_t rnd=gRandom->Exp(1);\r
 974        Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();\r
 975        his1->Fill(rnd); \r
 976        if (!isDownscaled) his2->Fill(rnd); \r
 977     }}\r
 978 \r
 979    */\r
 980 \r
 981 }\r
 982 \r
 983 \r
 984 \r
 985 \r
 986 \r
 987 \r
 988 \r
 989 \r
 990 \r
 991 //_____________________________________________________________________________\r
 992 Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
 993 {\r
 994  // Constrain TPC inner params constrained\r
 995  //\r
 996       if(!tpcInnerC) return kFALSE; \r
 997       if(!vtx) return kFALSE;\r
 998 \r
 999       Double_t dz[2],cov[3];\r
1000       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1001       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1002       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1003       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1004 \r
1005 \r
1006       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1007       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
1008       Double_t c[3]={covar[2],0.,covar[5]};\r
1009       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
1010       if (chi2C>kVeryBig) return kFALSE; \r
1011 \r
1012       if(!tpcInnerC->Update(p,c)) return kFALSE;\r
1013 \r
1014   return kTRUE;\r
1015 }\r
1016 \r
1017 //_____________________________________________________________________________\r
1018 Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
1019 {\r
1020  // Constrain TPC inner params constrained\r
1021  //\r
1022       if(!trackInnerC) return kFALSE; \r
1023       if(!vtx) return kFALSE;\r
1024 \r
1025       const Double_t kRadius  = 2.8; \r
1026       const Double_t kMaxStep = 1.0; \r
1027 \r
1028       Double_t dz[2],cov[3];\r
1029 \r
1030       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1031       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1032       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1033 \r
1034       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
1035       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1036 \r
1037       //\r
1038       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1039       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
1040       Double_t c[3]={covar[2],0.,covar[5]};\r
1041       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
1042       if (chi2C>kVeryBig) return kFALSE; \r
1043 \r
1044       if(!trackInnerC->Update(p,c)) return kFALSE;\r
1045 \r
1046   return kTRUE;\r
1047 }\r
1048 \r
1049 \r
1050 //_____________________________________________________________________________\r
1051 TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
1052 {\r
1053   if(!particle) return NULL;\r
1054   if(!stack) return NULL;\r
1055 \r
1056   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1057   TParticle* mother = NULL; \r
1058   mother = stack->Particle(motherLabel); \r
1059 \r
1060 return mother;\r
1061 }\r
1062 \r
1063 //_____________________________________________________________________________\r
1064 Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
1065 {\r
1066   Bool_t isFromConversion = kFALSE;\r
1067 \r
1068   if(stack) {\r
1069     TParticle* particle = stack->Particle(label);\r
1070 \r
1071     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1072     {\r
1073        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1074        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1075 \r
1076        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1077        TParticle* mother = stack->Particle(motherLabel); \r
1078        if(mother) {\r
1079           Int_t motherPdg = mother->GetPdgCode();\r
1080 \r
1081           if(!isPrim && mech==5 && motherPdg==kGamma) { \r
1082             isFromConversion=kTRUE; \r
1083           }\r
1084        }\r
1085     } \r
1086   } \r
1087 \r
1088   return isFromConversion;\r
1089 }\r
1090 \r
1091 //_____________________________________________________________________________\r
1092 Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
1093 {\r
1094   Bool_t isFromMaterial = kFALSE;\r
1095 \r
1096   if(stack) {\r
1097     TParticle* particle = stack->Particle(label);\r
1098 \r
1099     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1100     {\r
1101        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1102        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1103 \r
1104        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1105        TParticle* mother = stack->Particle(motherLabel); \r
1106        if(mother) {\r
1107           if(!isPrim && mech==13) { \r
1108             isFromMaterial=kTRUE; \r
1109           }\r
1110        }\r
1111      } \r
1112   } \r
1113 \r
1114   return isFromMaterial;\r
1115 }\r
1116 \r
1117 //_____________________________________________________________________________\r
1118 Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
1119 {\r
1120   Bool_t isFromStrangeness = kFALSE;\r
1121 \r
1122   if(stack) {\r
1123     TParticle* particle = stack->Particle(label);\r
1124 \r
1125     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1126     {\r
1127        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1128        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1129 \r
1130        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1131        TParticle* mother = stack->Particle(motherLabel); \r
1132        if(mother) {\r
1133           Int_t motherPdg = mother->GetPdgCode();\r
1134 \r
1135           // K+-, lambda, antilambda, K0s decays\r
1136           if(!isPrim && mech==4 && \r
1137               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
1138           {\r
1139             isFromStrangeness = kTRUE;\r
1140           } \r
1141        }\r
1142     } \r
1143   } \r
1144 \r
1145   return isFromStrangeness;\r
1146 }\r
1147 \r
1148 \r
1149 //_____________________________________________________________________________\r
1150 void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
1151 {\r
1152   //\r
1153   // Called one at the end \r
1154   // locally on working node\r
1155   //\r
1156 \r
1157   // Post output data.\r
1158   PostData(1, fOutput);\r
1159   //PostData(0, fOutputSummary);\r
1160 }\r
1161 \r
1162 //_____________________________________________________________________________\r
1163 void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
1164 {\r
1165   // Called one at the end \r
1166   if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;\r
1167   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(0));\r
1168   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
1169 \r
1170   TChain* chain = new TChain("dNdPtTree");\r
1171   if(!chain) return;\r
1172   chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");\r
1173   TTree *tree = chain->CopyTree("1");\r
1174   if (chain) { delete chain; chain=0; }\r
1175   if(!tree) return;\r
1176   tree->Print();\r
1177   fOutputSummary = tree;\r
1178 \r
1179   if (!fOutputSummary) {\r
1180     Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable GetOutputData(0)==0x0 ..." );\r
1181     return;\r
1182   }\r
1183 \r
1184   PostData(0, fOutputSummary);\r
1185   //PostData(1, fOutput);\r
1186 }\r