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