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