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       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
 542       if(!fTreeSRedirector) return;\r
 543       (*fTreeSRedirector)<<"highPt"<<\r
 544         "fileName.="<<&fileName<<\r
 545         "runNumber="<<runNumber<<\r
 546         "evtTimeStamp="<<evtTimeStamp<<\r
 547         "evtNumberInFile="<<evtNumberInFile<<\r
 548         "triggerClass="<<&triggerClass<<      //  trigger\r
 549         "Bz="<<bz<<\r
 550         "vtxESD.="<<vtxESD<<\r
 551         "IRtot="<<ir1<<\r
 552         "IRint2="<<ir2<<\r
 553         "mult="<<mult<<\r
 554         "esdTrack.="<<track<<\r
 555         "centralityF="<<centralityF<<\r
 556         "\n";\r
 557     }\r
 558   }\r
 559   \r
 560 }\r
 561 \r
 562 \r
 563 //_____________________________________________________________________________\r
 564 void AliAnalysisTaskFilteredTree::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const /*esdFriend*/)\r
 565 {\r
 566   //\r
 567   // Process laser events\r
 568   //\r
 569   if(!esdEvent) {\r
 570     AliDebug(AliLog::kError, "esdEvent not available");\r
 571     return;\r
 572   }\r
 573 \r
 574   // get file name\r
 575   TTree *chain = (TChain*)GetInputData(0);\r
 576   if(!chain) { \r
 577     Printf("ERROR: Could not receive input chain");\r
 578     return;\r
 579   }\r
 580   TObjString fileName(chain->GetCurrentFile()->GetName());\r
 581 \r
 582   // laser events \r
 583   const AliESDHeader* esdHeader = esdEvent->GetHeader();\r
 584   if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) \r
 585   {\r
 586     Int_t countLaserTracks = 0;\r
 587     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
 588     {\r
 589       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
 590       if(!track) continue;\r
 591 \r
 592       if(track->GetTPCInnerParam()) countLaserTracks++;\r
 593     }\r
 594        \r
 595     if(countLaserTracks > 100) \r
 596     {      \r
 597       Double_t runNumber = esdEvent->GetRunNumber();\r
 598       Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
 599       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
 600       Double_t bz = esdEvent->GetMagneticField();\r
 601       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
 602 \r
 603       if(!fTreeSRedirector) return;\r
 604       (*fTreeSRedirector)<<"Laser"<<\r
 605         "fileName.="<<&fileName<<\r
 606         "runNumber="<<runNumber<<\r
 607         "evtTimeStamp="<<evtTimeStamp<<\r
 608         "evtNumberInFile="<<evtNumberInFile<<\r
 609         "triggerClass="<<&triggerClass<<      //  trigger\r
 610         "Bz="<<bz<<\r
 611         "multTPCtracks="<<countLaserTracks<<\r
 612         "\n";\r
 613     }\r
 614   }\r
 615 }\r
 616 \r
 617 //_____________________________________________________________________________\r
 618 void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
 619 {\r
 620   //\r
 621   // Select real events with high-pT tracks\r
 622   // Calculate and stor in the output tree:\r
 623   //  TPC constrained tracks\r
 624   //  InnerParams constrained tracks\r
 625   //  TPC-ITS tracks\r
 626   //  ITSout-InnerParams tracks\r
 627   //  chi2 distance between TPC constrained and TPC-ITS tracks\r
 628   //  chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
 629   //  chi2 distance between ITSout and InnerParams tracks\r
 630   //  MC information: \r
 631   //   track references at ITSin, TPCin; InnerParam at first TPC track reference, \r
 632   //   particle ID, mother ID, production mechanism ...\r
 633   // \r
 634 \r
 635   if(!esdEvent) {\r
 636     AliDebug(AliLog::kError, "esdEvent not available");\r
 637     return;\r
 638   }\r
 639 \r
 640   // get selection cuts\r
 641   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
 642   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
 643   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
 644 \r
 645   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
 646     AliDebug(AliLog::kError, "cuts not available");\r
 647     return;\r
 648   }\r
 649 \r
 650   // trigger selection\r
 651   Bool_t isEventTriggered = kTRUE;\r
 652   AliPhysicsSelection *physicsSelection = NULL;\r
 653   AliTriggerAnalysis* triggerAnalysis = NULL;\r
 654 \r
 655   // \r
 656   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
 657   if (!inputHandler)\r
 658   {\r
 659     Printf("ERROR: Could not receive input handler");\r
 660     return;\r
 661   }\r
 662    \r
 663   // get file name\r
 664   TTree *chain = (TChain*)GetInputData(0);\r
 665   if(!chain) { \r
 666     Printf("ERROR: Could not receive input chain");\r
 667     return;\r
 668   }\r
 669   TObjString fileName(chain->GetCurrentFile()->GetName());\r
 670 \r
 671   // trigger\r
 672   if(evtCuts->IsTriggerRequired())  \r
 673   {\r
 674     // always MB\r
 675     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
 676 \r
 677     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
 678     if(!physicsSelection) return;\r
 679     //SetPhysicsTriggerSelection(physicsSelection);\r
 680 \r
 681     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
 682       // set trigger (V0AND)\r
 683       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
 684       if(!triggerAnalysis) return;\r
 685       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
 686     }\r
 687   }\r
 688 \r
 689   // centrality determination\r
 690   Float_t centralityF = -1;\r
 691   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
 692   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 693 \r
 694   // use MC information\r
 695   AliHeader* header = 0;\r
 696   AliGenEventHeader* genHeader = 0;\r
 697   AliStack* stack = 0;\r
 698   TArrayF vtxMC(3);\r
 699 \r
 700   Int_t multMCTrueTracks = 0;\r
 701   if(IsUseMCInfo())\r
 702   {\r
 703     //\r
 704     if(!mcEvent) {\r
 705       AliDebug(AliLog::kError, "mcEvent not available");\r
 706       return;\r
 707     }\r
 708     // get MC event header\r
 709     header = mcEvent->Header();\r
 710     if (!header) {\r
 711       AliDebug(AliLog::kError, "Header not available");\r
 712       return;\r
 713     }\r
 714     // MC particle stack\r
 715     stack = mcEvent->Stack();\r
 716     if (!stack) {\r
 717       AliDebug(AliLog::kError, "Stack not available");\r
 718       return;\r
 719     }\r
 720 \r
 721     // get MC vertex\r
 722     genHeader = header->GenEventHeader();\r
 723     if (!genHeader) {\r
 724       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
 725       return;\r
 726     }\r
 727     genHeader->PrimaryVertex(vtxMC);\r
 728 \r
 729     // multipliticy of all MC primary tracks\r
 730     // in Zv, pt and eta ranges)\r
 731     multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
 732 \r
 733   } // end bUseMC\r
 734 \r
 735   // get reconstructed vertex  \r
 736   //const AliESDVertex* vtxESD = 0; \r
 737   AliESDVertex* vtxESD = 0; \r
 738   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
 739         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
 740   }\r
 741   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
 742      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
 743   }\r
 744   else {\r
 745         return;\r
 746   }\r
 747 \r
 748   if(!vtxESD) return;\r
 749 \r
 750   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
 751   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
 752   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
 753 \r
 754 \r
 755   // check event cuts\r
 756   if(isEventOK && isEventTriggered)\r
 757   {\r
 758     //\r
 759     // get IR information\r
 760     //\r
 761     AliESDHeader *esdHeader = 0;\r
 762     esdHeader = esdEvent->GetHeader();\r
 763     if(!esdHeader) return;\r
 764     //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s\r
 765     //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval\r
 766     //\r
 767     Int_t ir1 = 0;\r
 768     Int_t ir2 = 0;\r
 769 \r
 770     //\r
 771     Double_t vert[3] = {0}; \r
 772     vert[0] = vtxESD->GetXv();\r
 773     vert[1] = vtxESD->GetYv();\r
 774     vert[2] = vtxESD->GetZv();\r
 775     Int_t mult = vtxESD->GetNContributors();\r
 776     Double_t bz = esdEvent->GetMagneticField();\r
 777     Double_t runNumber = esdEvent->GetRunNumber();\r
 778     Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
 779     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
 780 \r
 781     // high pT tracks\r
 782     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
 783     {\r
 784       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
 785       if(!track) continue;\r
 786       if(track->Charge()==0) continue;\r
 787       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
 788       if(!accCuts->AcceptTrack(track)) continue;\r
 789       \r
 790       // downscale low-pT tracks\r
 791       Double_t scalempt= TMath::Min(track->Pt(),10.);\r
 792       Double_t downscaleF = gRandom->Rndm();\r
 793       downscaleF *= fLowPtTrackDownscaligF;\r
 794       if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
 795       //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
 796 \r
 797       // Dump to the tree \r
 798       // vertex\r
 799       // TPC constrained tracks\r
 800       // InnerParams constrained tracks\r
 801       // TPC-ITS tracks\r
 802       // ITSout-InnerParams tracks\r
 803       // chi2 distance between TPC constrained and TPC-ITS tracks\r
 804       // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
 805       // chi2 distance between ITSout and InnerParams tracks\r
 806       // MC information\r
 807       \r
 808       Double_t x[3]; track->GetXYZ(x);\r
 809       Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
 810 \r
 811       //\r
 812       // Transform TPC inner params to track reference frame\r
 813       //\r
 814       Bool_t isOKtpcInner = kFALSE;\r
 815       AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
 816       if (tpcInner) {\r
 817         // transform to the track reference frame \r
 818         isOKtpcInner = tpcInner->Rotate(track->GetAlpha());\r
 819         isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 820       }\r
 821 \r
 822       //\r
 823       // Constrain TPC inner to vertex\r
 824       // clone TPCinner has to be deleted\r
 825       //\r
 826       Bool_t isOKtpcInnerC = kFALSE;\r
 827       AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
 828       if (tpcInnerC) {\r
 829         isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
 830         isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());\r
 831         isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 832       }\r
 833 \r
 834       //\r
 835       // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex  \r
 836       // Clone track InnerParams has to be deleted\r
 837       //\r
 838       Bool_t isOKtrackInnerC = kFALSE;\r
 839       AliExternalTrackParam * trackInnerC =  new AliExternalTrackParam(*(track->GetInnerParam()));\r
 840       if (trackInnerC) {\r
 841         isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
 842         isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());\r
 843         isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
 844       } \r
 845       \r
 846       //\r
 847       // calculate chi2 between vi and vj vectors\r
 848       // with covi and covj covariance matrices\r
 849       // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
 850       //\r
 851       TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
 852       TMatrixD delta(1,5),  deltatrackC(1,5);\r
 853       TMatrixD covarM(5,5), covarMtrackC(5,5);\r
 854       TMatrixD chi2(1,1);\r
 855       TMatrixD chi2trackC(1,1);\r
 856 \r
 857       if(isOKtpcInnerC && isOKtrackInnerC) \r
 858       {\r
 859         for (Int_t ipar=0; ipar<5; ipar++) {\r
 860           deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 861           delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 862 \r
 863           deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 864           deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
 865 \r
 866           for (Int_t jpar=0; jpar<5; jpar++) {\r
 867             Int_t index=track->GetIndex(ipar,jpar);\r
 868             covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
 869             covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
 870           }\r
 871         }\r
 872 \r
 873         // chi2 distance TPC constrained and TPC+ITS\r
 874         TMatrixD covarMInv = covarM.Invert();\r
 875         TMatrixD mat2 = covarMInv*deltaT;\r
 876         chi2 = delta*mat2; \r
 877         //chi2.Print();\r
 878 \r
 879         // chi2 distance TPC refitted constrained and TPC+ITS\r
 880         TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
 881         TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
 882         chi2trackC = deltatrackC*mat2trackC; \r
 883         //chi2trackC.Print();\r
 884       }\r
 885 \r
 886 \r
 887       //\r
 888       // Propagate ITSout to TPC inner wall \r
 889       // and calculate chi2 distance to track (InnerParams)\r
 890       //\r
 891       const Double_t kTPCRadius=85; \r
 892       const Double_t kStep=3; \r
 893 \r
 894       // clone track InnerParams has to be deleted\r
 895       Bool_t isOKtrackInnerC2 = kFALSE;\r
 896       AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
 897       if (trackInnerC2) {\r
 898         isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);\r
 899       }\r
 900 \r
 901       Bool_t isOKouterITSc = kFALSE;\r
 902       AliExternalTrackParam *outerITSc = NULL;\r
 903       TMatrixD chi2OuterITS(1,1);\r
 904 \r
 905       if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
 906       {\r
 907         // propagate ITSout to TPC inner wall\r
 908         AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
 909 \r
 910         if(friendTrack) \r
 911         {\r
 912           outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut());\r
 913           if(outerITSc) \r
 914           {\r
 915             isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);\r
 916             isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
 917             isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
 918 \r
 919             //\r
 920             // calculate chi2 between outerITS and innerParams\r
 921             // cov without z-coordinate at the moment\r
 922             //\r
 923             TMatrixD deltaTouterITS(4,1);\r
 924             TMatrixD deltaouterITS(1,4);\r
 925             TMatrixD covarMouterITS(4,4);\r
 926 \r
 927             if(isOKtrackInnerC2 && isOKouterITSc) {\r
 928               Int_t kipar = 0;\r
 929               Int_t kjpar = 0;\r
 930               for (Int_t ipar=0; ipar<5; ipar++) {\r
 931                 if(ipar!=1) {\r
 932                   deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
 933                   deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
 934                 }\r
 935 \r
 936                 kjpar=0;\r
 937                 for (Int_t jpar=0; jpar<5; jpar++) {\r
 938                   Int_t index=outerITSc->GetIndex(ipar,jpar);\r
 939                   if(ipar !=1 || jpar!=1) {\r
 940                     covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
 941                   }\r
 942                   if(jpar!=1)  kjpar++;\r
 943                 }\r
 944                 if(ipar!=1) kipar++;\r
 945               }\r
 946 \r
 947               // chi2 distance ITSout and InnerParams\r
 948               TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
 949               TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
 950               chi2OuterITS = deltaouterITS*mat2outerITS; \r
 951               //chi2OuterITS.Print();\r
 952             } \r
 953           }\r
 954         }\r
 955       }\r
 956 \r
 957       //\r
 958       // MC info\r
 959       //\r
 960       TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
 961       TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
 962       Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
 963       Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
 964       Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
 965       Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
 966       Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
 967 \r
 968       AliTrackReference *refTPCIn = NULL;\r
 969       AliTrackReference *refITS = NULL;\r
 970 \r
 971       Bool_t isOKtrackInnerC3 = kFALSE;\r
 972       AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
 973 \r
 974       if(IsUseMCInfo()) \r
 975       {\r
 976         if(!stack) return;\r
 977 \r
 978         //\r
 979         // global track\r
 980         //\r
 981         Int_t label = TMath::Abs(track->GetLabel()); \r
 982         particle = stack->Particle(label);\r
 983         if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
 984         {\r
 985           particleMother = GetMother(particle,stack);\r
 986           mech = particle->GetUniqueID();\r
 987           isPrim = stack->IsPhysicalPrimary(label);\r
 988           isFromStrangess  = IsFromStrangeness(label,stack);\r
 989           isFromConversion = IsFromConversion(label,stack);\r
 990           isFromMaterial   = IsFromMaterial(label,stack);\r
 991         }\r
 992 \r
 993         //\r
 994         // TPC track\r
 995         //\r
 996         Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
 997         particleTPC = stack->Particle(labelTPC);\r
 998         if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
 999         {\r
1000           particleMotherTPC = GetMother(particleTPC,stack);\r
1001           mechTPC = particleTPC->GetUniqueID();\r
1002           isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
1003           isFromStrangessTPC  = IsFromStrangeness(labelTPC,stack);\r
1004           isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
1005           isFromMaterialTPC   = IsFromMaterial(labelTPC,stack);\r
1006         }\r
1007 \r
1008         //\r
1009         // store first track reference\r
1010         // for TPC track\r
1011         //\r
1012         TParticle *part=0;\r
1013         TClonesArray *trefs=0;\r
1014         Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
1015 \r
1016         if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
1017         {\r
1018           Int_t nTrackRef = trefs->GetEntries();\r
1019           //printf("nTrackRef %d \n",nTrackRef);\r
1020 \r
1021           Int_t countITS = 0;\r
1022           for (Int_t iref = 0; iref < nTrackRef; iref++) \r
1023           {\r
1024             AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
1025 \r
1026              // Ref. in the middle ITS \r
1027             if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
1028             {\r
1029               if(!refITS && countITS==2) {\r
1030                 refITS = ref;\r
1031                 //printf("refITS %p \n",refITS);\r
1032               }\r
1033               countITS++;\r
1034             }\r
1035 \r
1036             // TPC\r
1037             if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
1038             {\r
1039               if(!refTPCIn) {\r
1040                 refTPCIn = ref;\r
1041                 //printf("refTPCIn %p \n",refTPCIn);\r
1042                 //break;\r
1043               }\r
1044             }\r
1045           }\r
1046 \r
1047           // transform inner params to TrackRef\r
1048           // reference frame\r
1049           if(refTPCIn && trackInnerC3) \r
1050           {\r
1051             Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
1052             isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);\r
1053             isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
1054           }\r
1055         }\r
1056 \r
1057         //\r
1058         // ITS track\r
1059         //\r
1060         Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
1061         particleITS = stack->Particle(labelITS);\r
1062         if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
1063         {\r
1064           particleMotherITS = GetMother(particleITS,stack);\r
1065           mechITS = particleITS->GetUniqueID();\r
1066           isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
1067           isFromStrangessITS  = IsFromStrangeness(labelITS,stack);\r
1068           isFromConversionITS = IsFromConversion(labelITS,stack);\r
1069           isFromMaterialITS   = IsFromMaterial(labelITS,stack);\r
1070         }\r
1071       }\r
1072 \r
1073       //\r
1074       Bool_t dumpToTree=kFALSE;\r
1075       \r
1076       if(isOKtpcInnerC  && isOKtrackInnerC) dumpToTree = kTRUE;\r
1077       if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;\r
1078       if(fUseMCInfo     && isOKtrackInnerC3) dumpToTree = kTRUE;\r
1079       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1080 \r
1081       //\r
1082       if(fTreeSRedirector && dumpToTree) \r
1083       {\r
1084         (*fTreeSRedirector)<<"highPt"<<\r
1085           "fileName.="<<&fileName<<\r
1086           "runNumber="<<runNumber<<\r
1087           "evtTimeStamp="<<evtTimeStamp<<\r
1088           "evtNumberInFile="<<evtNumberInFile<<\r
1089           "triggerClass="<<&triggerClass<<      //  trigger\r
1090           "Bz="<<bz<<\r
1091           "vtxESD.="<<vtxESD<<\r
1092           "IRtot="<<ir1<<\r
1093           "IRint2="<<ir2<<\r
1094           "mult="<<mult<<\r
1095           "esdTrack.="<<track<<\r
1096           "extTPCInnerC.="<<tpcInnerC<<\r
1097           "extInnerParamC.="<<trackInnerC<<\r
1098           "extInnerParam.="<<trackInnerC2<<\r
1099           "extOuterITS.="<<outerITSc<<\r
1100           "extInnerParamRef.="<<trackInnerC3<<\r
1101           "refTPCIn.="<<refTPCIn<<\r
1102           "refITS.="<<refITS<<\r
1103           "chi2TPCInnerC="<<chi2(0,0)<<\r
1104           "chi2InnerC="<<chi2trackC(0,0)<<\r
1105           "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
1106           "centralityF="<<centralityF<<\r
1107           "particle.="<<particle<<\r
1108           "particleMother.="<<particleMother<<\r
1109           "mech="<<mech<<\r
1110           "isPrim="<<isPrim<<\r
1111           "isFromStrangess="<<isFromStrangess<<\r
1112           "isFromConversion="<<isFromConversion<<\r
1113           "isFromMaterial="<<isFromMaterial<<\r
1114           "particleTPC.="<<particleTPC<<\r
1115           "particleMotherTPC.="<<particleMotherTPC<<\r
1116           "mechTPC="<<mechTPC<<\r
1117           "isPrimTPC="<<isPrimTPC<<\r
1118           "isFromStrangessTPC="<<isFromStrangessTPC<<\r
1119           "isFromConversionTPC="<<isFromConversionTPC<<\r
1120           "isFromMaterialTPC="<<isFromMaterialTPC<<\r
1121           "particleITS.="<<particleITS<<\r
1122           "particleMotherITS.="<<particleMotherITS<<\r
1123           "mechITS="<<mechITS<<\r
1124           "isPrimITS="<<isPrimITS<<\r
1125           "isFromStrangessITS="<<isFromStrangessITS<<\r
1126           "isFromConversionITS="<<isFromConversionITS<<\r
1127           "isFromMaterialITS="<<isFromMaterialITS<<\r
1128           "\n";\r
1129         }\r
1130       \r
1131         if(tpcInnerC) delete tpcInnerC;\r
1132         if(trackInnerC) delete trackInnerC;\r
1133         if(trackInnerC2) delete trackInnerC2;\r
1134         if(outerITSc) delete outerITSc;\r
1135         if(trackInnerC3) delete trackInnerC3;\r
1136     }\r
1137   }\r
1138   \r
1139 }\r
1140 \r
1141 \r
1142 //_____________________________________________________________________________\r
1143 void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1144 {\r
1145   //\r
1146   // Fill tree for efficiency studies MC only\r
1147 \r
1148   if(!esdEvent) {\r
1149     AliDebug(AliLog::kError, "esdEvent not available");\r
1150     return;\r
1151   }\r
1152 \r
1153    if(!mcEvent) {\r
1154     AliDebug(AliLog::kError, "mcEvent not available");\r
1155     return;\r
1156   }\r
1157 \r
1158   // get selection cuts\r
1159   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1160   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1161   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1162 \r
1163   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1164     AliDebug(AliLog::kError, "cuts not available");\r
1165     return;\r
1166   }\r
1167 \r
1168   // trigger selection\r
1169   Bool_t isEventTriggered = kTRUE;\r
1170   AliPhysicsSelection *physicsSelection = NULL;\r
1171   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1172 \r
1173   // \r
1174   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1175   if (!inputHandler)\r
1176   {\r
1177     Printf("ERROR: Could not receive input handler");\r
1178     return;\r
1179   }\r
1180    \r
1181   // get file name\r
1182   TTree *chain = (TChain*)GetInputData(0);\r
1183   if(!chain) { \r
1184     Printf("ERROR: Could not receive input chain");\r
1185     return;\r
1186   }\r
1187   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1188 \r
1189   // trigger\r
1190   if(evtCuts->IsTriggerRequired())  \r
1191   {\r
1192     // always MB\r
1193     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1194 \r
1195     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1196     if(!physicsSelection) return;\r
1197 \r
1198     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1199       // set trigger (V0AND)\r
1200       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1201       if(!triggerAnalysis) return;\r
1202       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1203     }\r
1204   }\r
1205 \r
1206   // centrality determination\r
1207   Float_t centralityF = -1;\r
1208   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1209   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1210 \r
1211   // use MC information\r
1212   AliHeader* header = 0;\r
1213   AliGenEventHeader* genHeader = 0;\r
1214   AliStack* stack = 0;\r
1215   TArrayF vtxMC(3);\r
1216 \r
1217   Int_t multMCTrueTracks = 0;\r
1218   if(IsUseMCInfo())\r
1219   {\r
1220     //\r
1221     if(!mcEvent) {\r
1222       AliDebug(AliLog::kError, "mcEvent not available");\r
1223       return;\r
1224     }\r
1225     // get MC event header\r
1226     header = mcEvent->Header();\r
1227     if (!header) {\r
1228       AliDebug(AliLog::kError, "Header not available");\r
1229       return;\r
1230     }\r
1231     // MC particle stack\r
1232     stack = mcEvent->Stack();\r
1233     if (!stack) {\r
1234       AliDebug(AliLog::kError, "Stack not available");\r
1235       return;\r
1236     }\r
1237 \r
1238     // get MC vertex\r
1239     genHeader = header->GenEventHeader();\r
1240     if (!genHeader) {\r
1241       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
1242       return;\r
1243     }\r
1244     genHeader->PrimaryVertex(vtxMC);\r
1245 \r
1246     // multipliticy of all MC primary tracks\r
1247     // in Zv, pt and eta ranges)\r
1248     multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
1249 \r
1250   } // end bUseMC\r
1251 \r
1252   // get reconstructed vertex  \r
1253   //const AliESDVertex* vtxESD = 0; \r
1254   AliESDVertex* vtxESD = 0; \r
1255   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1256         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1257   }\r
1258   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1259      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1260   }\r
1261   else {\r
1262         return;\r
1263   }\r
1264 \r
1265   if(!vtxESD) return;\r
1266 \r
1267   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1268   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1269   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1270 \r
1271   TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1272 \r
1273   // check event cuts\r
1274   if(isEventOK && isEventTriggered)\r
1275   {\r
1276     if(IsUseMCInfo()) \r
1277     {\r
1278       if(!stack) return;\r
1279 \r
1280       //\r
1281       // MC info\r
1282       //\r
1283       TParticle *particle=NULL;\r
1284       TParticle *particleMother=NULL;\r
1285       Int_t mech=-1;\r
1286 \r
1287       // reco event info\r
1288       Double_t vert[3] = {0}; \r
1289       vert[0] = vtxESD->GetXv();\r
1290       vert[1] = vtxESD->GetYv();\r
1291       vert[2] = vtxESD->GetZv();\r
1292       Int_t mult = vtxESD->GetNContributors();\r
1293       Double_t bz = esdEvent->GetMagneticField();\r
1294       Double_t runNumber = esdEvent->GetRunNumber();\r
1295       Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1296       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
1297 \r
1298       // loop over MC stack\r
1299       for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc) \r
1300       {\r
1301          particle = stack->Particle(iMc);\r
1302          if (!particle)\r
1303          continue;\r
1304 \r
1305          // only charged particles\r
1306          if(!particle->GetPDG()) continue;\r
1307          Double_t charge = particle->GetPDG()->Charge()/3.;\r
1308          if (TMath::Abs(charge) < 0.001)\r
1309          continue;\r
1310 \r
1311          // only primary particles\r
1312          Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
1313          if(!prim) continue;\r
1314 \r
1315          // downscale low-pT particles\r
1316          Double_t scalempt= TMath::Min(particle->Pt(),10.);\r
1317          Double_t downscaleF = gRandom->Rndm();\r
1318          downscaleF *= fLowPtTrackDownscaligF;\r
1319          if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
1320 \r
1321          // is particle in acceptance\r
1322          if(!accCuts->AcceptTrack(particle)) continue;\r
1323        \r
1324          // check if particle reconstructed\r
1325          Bool_t isRec = kFALSE;\r
1326          Int_t  trackIndex = -1;\r
1327          for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
1328          {\r
1329            \r
1330            AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
1331            if(!track) continue;\r
1332            if(track->Charge()==0) continue;\r
1333            if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track)) \r
1334            {\r
1335              Int_t label =  TMath::Abs(track->GetLabel());\r
1336              if(label == iMc) {\r
1337                isRec = kTRUE;\r
1338                trackIndex = iTrack;\r
1339                break;\r
1340              }\r
1341            } \r
1342          }\r
1343 \r
1344          // Store information in the output tree\r
1345          AliESDtrack *recTrack = NULL; \r
1346          if(trackIndex>-1)  { \r
1347            recTrack = esdEvent->GetTrack(trackIndex); \r
1348          } else {\r
1349            recTrack = new AliESDtrack(); \r
1350          } \r
1351 \r
1352          particleMother = GetMother(particle,stack);\r
1353          mech = particle->GetUniqueID();\r
1354 \r
1355          //MC particle track length\r
1356          Double_t tpcTrackLength = 0.;\r
1357          AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);\r
1358          if(mcParticle) {\r
1359            Int_t counter;\r
1360            tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);\r
1361          } \r
1362 \r
1363 \r
1364          //\r
1365          if(fTreeSRedirector) {\r
1366            (*fTreeSRedirector)<<"MCEffTree"<<\r
1367            "fileName.="<<&fileName<<\r
1368            "triggerClass.="<<&triggerClass<<\r
1369            "runNumber="<<runNumber<<\r
1370            "evtTimeStamp="<<evtTimeStamp<<\r
1371            "evtNumberInFile="<<evtNumberInFile<<\r
1372            "Bz="<<bz<<\r
1373            "vtxESD.="<<vtxESD<<\r
1374            "mult="<<mult<<\r
1375            "esdTrack.="<<recTrack<<\r
1376            "isRec="<<isRec<<\r
1377            "tpcTrackLength="<<tpcTrackLength<<\r
1378            "particle.="<<particle<<\r
1379            "particleMother.="<<particleMother<<\r
1380            "mech="<<mech<<\r
1381            "\n";\r
1382          }\r
1383 \r
1384          if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;\r
1385       }\r
1386     }\r
1387   }\r
1388   \r
1389 }\r
1390 \r
1391 //_____________________________________________________________________________\r
1392 Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {\r
1393   //\r
1394   // check if particle is Z > 1 \r
1395   //\r
1396   if (track->GetTPCNcls() < 60) return kFALSE;\r
1397   Double_t mom = track->GetInnerParam()->GetP();\r
1398   if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization\r
1399   Float_t dca[2], bCov[3];\r
1400   track->GetImpactParameters(dca,bCov);\r
1401   //\r
1402 \r
1403   Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);\r
1404 \r
1405   if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;\r
1406 \r
1407   return kFALSE;\r
1408 }\r
1409 \r
1410 //_____________________________________________________________________________\r
1411 void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1412 {\r
1413   //\r
1414   // Select real events with V0 (K0s and Lambda) high-pT candidates\r
1415   //\r
1416   if(!esdEvent) {\r
1417     AliDebug(AliLog::kError, "esdEvent not available");\r
1418     return;\r
1419   }\r
1420 \r
1421   // get selection cuts\r
1422   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1423   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1424   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1425 \r
1426   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1427     AliDebug(AliLog::kError, "cuts not available");\r
1428     return;\r
1429   }\r
1430 \r
1431   // trigger selection\r
1432   Bool_t isEventTriggered = kTRUE;\r
1433   AliPhysicsSelection *physicsSelection = NULL;\r
1434   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1435 \r
1436   // \r
1437   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1438   if (!inputHandler)\r
1439   {\r
1440     Printf("ERROR: Could not receive input handler");\r
1441     return;\r
1442   }\r
1443    \r
1444   // get file name\r
1445   TTree *chain = (TChain*)GetInputData(0);\r
1446   if(!chain) { \r
1447     Printf("ERROR: Could not receive input chain");\r
1448     return;\r
1449   }\r
1450   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1451 \r
1452   // trigger\r
1453   if(evtCuts->IsTriggerRequired())  \r
1454   {\r
1455     // always MB\r
1456     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1457 \r
1458     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1459     if(!physicsSelection) return;\r
1460     //SetPhysicsTriggerSelection(physicsSelection);\r
1461 \r
1462     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1463       // set trigger (V0AND)\r
1464       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1465       if(!triggerAnalysis) return;\r
1466       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1467     }\r
1468   }\r
1469 \r
1470   // centrality determination\r
1471   Float_t centralityF = -1;\r
1472   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1473   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1474 \r
1475 \r
1476   // get reconstructed vertex  \r
1477   //const AliESDVertex* vtxESD = 0; \r
1478   AliESDVertex* vtxESD = 0; \r
1479   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1480         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1481   }\r
1482   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1483      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1484   }\r
1485   else {\r
1486         return;\r
1487   }\r
1488 \r
1489   if(!vtxESD) return;\r
1490 \r
1491   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1492   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1493   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1494 \r
1495   // check event cuts\r
1496   if(isEventOK && isEventTriggered) {\r
1497   //\r
1498   // Dump the pt downscaled V0 into the tree\r
1499   // \r
1500   Int_t ntracks = esdEvent->GetNumberOfTracks();\r
1501   Int_t nV0s = esdEvent->GetNumberOfV0s();\r
1502   Int_t run = esdEvent->GetRunNumber();\r
1503   Int_t time= esdEvent->GetTimeStamp();\r
1504   Int_t evNr=esdEvent->GetEventNumberInFile();\r
1505   Double_t bz = esdEvent->GetMagneticField();\r
1506 \r
1507 \r
1508   for (Int_t iv0=0; iv0<nV0s; iv0++){\r
1509     AliESDv0 * v0 = esdEvent->GetV0(iv0);\r
1510     if (!v0) continue;\r
1511     AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));\r
1512     AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));\r
1513     if (!track0) continue;\r
1514     if (!track1) continue;\r
1515     if (track0->GetSign()<0) {\r
1516       track1 = esdEvent->GetTrack(v0->GetIndex(0));\r
1517       track0 = esdEvent->GetTrack(v0->GetIndex(1));\r
1518     }\r
1519     //\r
1520     Bool_t isDownscaled = IsV0Downscaled(v0);\r
1521     if (isDownscaled) continue;\r
1522     AliKFParticle kfparticle; //\r
1523     Int_t type=GetKFParticle(v0,esdEvent,kfparticle);\r
1524     if (type==0) continue;   \r
1525     TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1526 \r
1527     if(!fTreeSRedirector) return;\r
1528     (*fTreeSRedirector)<<"V0s"<<\r
1529       "isDownscaled="<<isDownscaled<<\r
1530       "triggerClass="<<&triggerClass<<      //  trigger\r
1531       "Bz="<<bz<<\r
1532       "fileName.="<<&fileName<<\r
1533       "runNumber="<<run<<\r
1534       "evtTimeStamp="<<time<<\r
1535       "evtNumberInFile="<<evNr<<\r
1536       "type="<<type<<\r
1537       "ntracks="<<ntracks<<\r
1538       "v0.="<<v0<<\r
1539       "kf.="<<&kfparticle<<\r
1540       "track0.="<<track0<<\r
1541       "track1.="<<track1<<\r
1542       "centralityF="<<centralityF<<\r
1543       "\n";\r
1544   }\r
1545   }\r
1546 }\r
1547 \r
1548 //_____________________________________________________________________________\r
1549 void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1550 {\r
1551   //\r
1552   // Select real events with large TPC dEdx signal\r
1553   //\r
1554   if(!esdEvent) {\r
1555     AliDebug(AliLog::kError, "esdEvent not available");\r
1556     return;\r
1557   }\r
1558 \r
1559   // get selection cuts\r
1560   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1561   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1562   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1563 \r
1564   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1565     AliDebug(AliLog::kError, "cuts not available");\r
1566     return;\r
1567   }\r
1568 \r
1569   // get file name\r
1570   TTree *chain = (TChain*)GetInputData(0);\r
1571   if(!chain) { \r
1572     Printf("ERROR: Could not receive input chain");\r
1573     return;\r
1574   }\r
1575   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1576 \r
1577   // trigger\r
1578   Bool_t isEventTriggered = kTRUE;\r
1579   AliPhysicsSelection *physicsSelection = NULL;\r
1580   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1581 \r
1582   // \r
1583   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1584   if (!inputHandler)\r
1585   {\r
1586     Printf("ERROR: Could not receive input handler");\r
1587     return;\r
1588   }\r
1589 \r
1590   if(evtCuts->IsTriggerRequired())  \r
1591   {\r
1592     // always MB\r
1593     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1594 \r
1595     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1596     if(!physicsSelection) return;\r
1597 \r
1598     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1599       // set trigger (V0AND)\r
1600       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1601       if(!triggerAnalysis) return;\r
1602       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1603     }\r
1604   }\r
1605 \r
1606   // get reconstructed vertex  \r
1607   AliESDVertex* vtxESD = 0; \r
1608   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1609         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1610   }\r
1611   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1612      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1613   }\r
1614   else {\r
1615         return;\r
1616   }\r
1617   if(!vtxESD) return;\r
1618 \r
1619   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1620   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1621   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1622 \r
1623 \r
1624   // check event cuts\r
1625   if(isEventOK && isEventTriggered)\r
1626   {\r
1627     Double_t vert[3] = {0}; \r
1628     vert[0] = vtxESD->GetXv();\r
1629     vert[1] = vtxESD->GetYv();\r
1630     vert[2] = vtxESD->GetZv();\r
1631     Int_t mult = vtxESD->GetNContributors();\r
1632     Double_t bz = esdEvent->GetMagneticField();\r
1633     Double_t runNumber = esdEvent->GetRunNumber();\r
1634     Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1635     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
1636 \r
1637     // large dEdx\r
1638     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
1639     {\r
1640       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
1641       if(!track) continue;\r
1642       if(track->Charge()==0) continue;\r
1643       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
1644       if(!accCuts->AcceptTrack(track)) continue;\r
1645 \r
1646       if(!IsHighDeDxParticle(track)) continue;\r
1647       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1648 \r
1649       if(!fTreeSRedirector) return;\r
1650       (*fTreeSRedirector)<<"dEdx"<<\r
1651       "fileName.="<<&fileName<<\r
1652       "runNumber="<<runNumber<<\r
1653       "evtTimeStamp="<<evtTimeStamp<<\r
1654       "evtNumberInFile="<<evtNumberInFile<<\r
1655         "triggerClass="<<&triggerClass<<      //  trigger\r
1656       "Bz="<<bz<<\r
1657       "vtxESD.="<<vtxESD<<\r
1658       "mult="<<mult<<\r
1659       "esdTrack.="<<track<<\r
1660       "\n";\r
1661     }\r
1662   }\r
1663 }\r
1664 \r
1665 //_____________________________________________________________________________\r
1666 Int_t   AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)\r
1667 {\r
1668   //\r
1669   // Create KF particle in case the V0 fullfill selection criteria\r
1670   //\r
1671   // Selection criteria\r
1672   //  0. algorithm cut\r
1673   //  1. track cut\r
1674   //  3. chi2 cut\r
1675   //  4. rough mass cut\r
1676   //  5. Normalized pointing angle cut\r
1677   //\r
1678   const Double_t cutMass=0.2;\r
1679   const Double_t kSigmaDCACut=3;\r
1680   //\r
1681   // 0.) algo cut - accept only on the fly\r
1682   //\r
1683   if (v0->GetOnFlyStatus() ==kFALSE) return 0;     \r
1684   //\r
1685   // 1.) track cut\r
1686   // \r
1687   AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));\r
1688   AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));\r
1689   /*\r
1690     TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";\r
1691     TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";\r
1692     TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";\r
1693   */  \r
1694   if (TMath::Abs(track0->GetTgl())>1) return 0;\r
1695   if (TMath::Abs(track1->GetTgl())>1) return 0;\r
1696   if ((track0->GetTPCClusterInfo(2,1))<100) return 0;\r
1697   if ((track1->GetTPCClusterInfo(2,1))<100) return 0;\r
1698   //if ((track0->GetITSclusters(0))<2) return 0;\r
1699   //if ((track1->GetITSclusters(0))<2) return 0; \r
1700   Float_t pos0[2]={0}, cov0[3]={0};\r
1701   Float_t pos1[2]={0}, cov1[3]={0};\r
1702   track0->GetImpactParameters(pos0,cov0);\r
1703   track0->GetImpactParameters(pos1,cov1);\r
1704   //\r
1705   if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;\r
1706   if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;\r
1707   // \r
1708   //\r
1709   // 3.) Chi2 cut\r
1710   //\r
1711   Double_t chi2KF = v0->GetKFInfo(2,2,2);\r
1712   if (chi2KF>25) return 0;\r
1713   //\r
1714   // 4.) Rough mass cut - 0.200 GeV\r
1715   //\r
1716   static Double_t masses[2]={-1};\r
1717   if (masses[0]<0){\r
1718     masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();\r
1719     masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();\r
1720   }\r
1721   Double_t mass00=  v0->GetEffMass(0,0);\r
1722   Double_t mass22=  v0->GetEffMass(2,2);\r
1723   Double_t mass42=  v0->GetEffMass(4,2);\r
1724   Double_t mass24=  v0->GetEffMass(2,4);\r
1725   Bool_t massOK=kFALSE;\r
1726   Int_t type=0;\r
1727   Int_t ptype=0;\r
1728   Double_t dmass=1;\r
1729   Int_t p1=0, p2=0;\r
1730   if (TMath::Abs(mass00-0)<cutMass) {\r
1731     massOK=kTRUE; type+=1; \r
1732     if (TMath::Abs(mass00-0)<dmass) {\r
1733       ptype=1;\r
1734       dmass=TMath::Abs(mass00-0);      \r
1735       p1=0; p2=0;\r
1736     } \r
1737   }\r
1738   if (TMath::Abs(mass24-masses[1])<cutMass) {\r
1739     massOK=kTRUE; type+=2; \r
1740     if (TMath::Abs(mass24-masses[1])<dmass){\r
1741       dmass = TMath::Abs(mass24-masses[1]);\r
1742       ptype=2;\r
1743       p1=2; p2=4;\r
1744     }\r
1745   }\r
1746   if (TMath::Abs(mass42-masses[1])<cutMass) {\r
1747     massOK=kTRUE; type+=4;\r
1748     if (TMath::Abs(mass42-masses[1])<dmass){\r
1749       dmass = TMath::Abs(mass42-masses[1]);\r
1750       ptype=4;\r
1751       p1=4; p2=2;\r
1752     }\r
1753   }\r
1754   if (TMath::Abs(mass22-masses[0])<cutMass) {\r
1755     massOK=kTRUE; type+=8;\r
1756     if (TMath::Abs(mass22-masses[0])<dmass){\r
1757       dmass = TMath::Abs(mass22-masses[0]);\r
1758       ptype=8;\r
1759       p1=2; p2=2;\r
1760     }\r
1761   }\r
1762   if (type==0) return 0;\r
1763   //\r
1764   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};\r
1765   const AliExternalTrackParam *paramP = v0->GetParamP();\r
1766   const AliExternalTrackParam *paramN = v0->GetParamN();\r
1767   if (paramP->GetSign()<0){\r
1768     paramP=v0->GetParamP();\r
1769     paramN=v0->GetParamN();\r
1770   }\r
1771   //Double_t *pparam1 = (Double_t*)paramP->GetParameter();\r
1772   //Double_t *pparam2 = (Double_t*)paramN->GetParameter();\r
1773   //\r
1774   AliKFParticle kfp1( *paramP, spdg[p1]  );\r
1775   AliKFParticle kfp2( *paramN, -1 *spdg[p2]  );\r
1776   AliKFParticle V0KF;\r
1777   (V0KF)+=kfp1;\r
1778   (V0KF)+=kfp2;\r
1779   kfparticle=V0KF;\r
1780   //\r
1781   // Pointing angle\r
1782   //\r
1783   Double_t  errPhi    = V0KF.GetErrPhi();\r
1784   Double_t  pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());\r
1785   if (pointAngle/errPhi>10) return 0;  \r
1786   //\r
1787   return ptype;  \r
1788 }\r
1789 \r
1790 //_____________________________________________________________________________\r
1791 Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)\r
1792 {\r
1793   //\r
1794   // Downscale randomly low pt V0\r
1795   //\r
1796   //return kFALSE;\r
1797   Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());\r
1798   Double_t scalempt= TMath::Min(maxPt,10.);\r
1799   Double_t downscaleF = gRandom->Rndm();\r
1800   downscaleF *= fLowPtV0DownscaligF;\r
1801   \r
1802   //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
1803   if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;\r
1804   return kFALSE;\r
1805 \r
1806   /*\r
1807     TH1F his1("his1","his1",100,0,10);\r
1808     TH1F his2("his2","his2",100,0,10);\r
1809     {for (Int_t i=0; i<10000; i++){\r
1810        Double_t rnd=gRandom->Exp(1);\r
1811        Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();\r
1812        his1->Fill(rnd); \r
1813        if (!isDownscaled) his2->Fill(rnd); \r
1814     }}\r
1815 \r
1816    */\r
1817 \r
1818 }\r
1819 \r
1820 \r
1821 \r
1822 //_____________________________________________________________________________\r
1823 Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
1824 {\r
1825  // Constrain TPC inner params constrained\r
1826  //\r
1827       if(!tpcInnerC) return kFALSE; \r
1828       if(!vtx) return kFALSE;\r
1829 \r
1830       Double_t dz[2],cov[3];\r
1831       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1832       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1833       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1834       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1835 \r
1836 \r
1837       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1838       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
1839       Double_t c[3]={covar[2],0.,covar[5]};\r
1840       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
1841       if (chi2C>kVeryBig) return kFALSE; \r
1842 \r
1843       if(!tpcInnerC->Update(p,c)) return kFALSE;\r
1844 \r
1845   return kTRUE;\r
1846 }\r
1847 \r
1848 //_____________________________________________________________________________\r
1849 Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
1850 {\r
1851  // Constrain TPC inner params constrained\r
1852  //\r
1853       if(!trackInnerC) return kFALSE; \r
1854       if(!vtx) return kFALSE;\r
1855 \r
1856       const Double_t kRadius  = 2.8; \r
1857       const Double_t kMaxStep = 1.0; \r
1858 \r
1859       Double_t dz[2],cov[3];\r
1860 \r
1861       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1862       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1863       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1864 \r
1865       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
1866       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1867 \r
1868       //\r
1869       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1870       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
1871       Double_t c[3]={covar[2],0.,covar[5]};\r
1872       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
1873       if (chi2C>kVeryBig) return kFALSE; \r
1874 \r
1875       if(!trackInnerC->Update(p,c)) return kFALSE;\r
1876 \r
1877   return kTRUE;\r
1878 }\r
1879 \r
1880 \r
1881 //_____________________________________________________________________________\r
1882 TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack) \r
1883 {\r
1884   if(!particle) return NULL;\r
1885   if(!stack) return NULL;\r
1886 \r
1887   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1888   TParticle* mother = NULL; \r
1889   mother = stack->Particle(motherLabel); \r
1890 \r
1891 return mother;\r
1892 }\r
1893 \r
1894 //_____________________________________________________________________________\r
1895 Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack) \r
1896 {\r
1897   Bool_t isFromConversion = kFALSE;\r
1898 \r
1899   if(stack) {\r
1900     TParticle* particle = stack->Particle(label);\r
1901 \r
1902     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1903     {\r
1904        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1905        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1906 \r
1907        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1908        TParticle* mother = stack->Particle(motherLabel); \r
1909        if(mother) {\r
1910           Int_t motherPdg = mother->GetPdgCode();\r
1911 \r
1912           if(!isPrim && mech==5 && motherPdg==kGamma) { \r
1913             isFromConversion=kTRUE; \r
1914           }\r
1915        }\r
1916     } \r
1917   } \r
1918 \r
1919   return isFromConversion;\r
1920 }\r
1921 \r
1922 //_____________________________________________________________________________\r
1923 Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack) \r
1924 {\r
1925   Bool_t isFromMaterial = kFALSE;\r
1926 \r
1927   if(stack) {\r
1928     TParticle* particle = stack->Particle(label);\r
1929 \r
1930     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1931     {\r
1932        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1933        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1934 \r
1935        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1936        TParticle* mother = stack->Particle(motherLabel); \r
1937        if(mother) {\r
1938           if(!isPrim && mech==13) { \r
1939             isFromMaterial=kTRUE; \r
1940           }\r
1941        }\r
1942      } \r
1943   } \r
1944 \r
1945   return isFromMaterial;\r
1946 }\r
1947 \r
1948 //_____________________________________________________________________________\r
1949 Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
1950 {\r
1951   Bool_t isFromStrangeness = kFALSE;\r
1952 \r
1953   if(stack) {\r
1954     TParticle* particle = stack->Particle(label);\r
1955 \r
1956     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1957     {\r
1958        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1959        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1960 \r
1961        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1962        TParticle* mother = stack->Particle(motherLabel); \r
1963        if(mother) {\r
1964           Int_t motherPdg = mother->GetPdgCode();\r
1965 \r
1966           // K+-, lambda, antilambda, K0s decays\r
1967           if(!isPrim && mech==4 && \r
1968               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
1969           {\r
1970             isFromStrangeness = kTRUE;\r
1971           } \r
1972        }\r
1973     } \r
1974   } \r
1975 \r
1976   return isFromStrangeness;\r
1977 }\r
1978 \r
1979 \r
1980 //_____________________________________________________________________________\r
1981 void AliAnalysisTaskFilteredTree::FinishTaskOutput() \r
1982 {\r
1983   //\r
1984   // Called one at the end \r
1985   // locally on working node\r
1986   //\r
1987 \r
1988   // must be deleted to store trees\r
1989   if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;\r
1990 \r
1991   // open temporary file and copy trees to the ouptut container\r
1992 \r
1993   TChain* chain = 0;\r
1994   //\r
1995   chain = new TChain("highPt");\r
1996   if(chain) { \r
1997     chain->Add("jotwinow_Temp_Trees.root");\r
1998     fHighPtTree = chain->CopyTree("1");\r
1999     delete chain; chain=0; \r
2000   }\r
2001   if(fHighPtTree) fHighPtTree->Print();\r
2002 \r
2003   //\r
2004   chain = new TChain("V0s");\r
2005   if(chain) { \r
2006     chain->Add("jotwinow_Temp_Trees.root");\r
2007     fV0Tree = chain->CopyTree("1");\r
2008     delete chain; chain=0; \r
2009   }\r
2010   if(fV0Tree) fV0Tree->Print();\r
2011 \r
2012   //\r
2013   chain = new TChain("dEdx");\r
2014   if(chain) { \r
2015     chain->Add("jotwinow_Temp_Trees.root");\r
2016     fdEdxTree = chain->CopyTree("1");\r
2017     delete chain; chain=0; \r
2018   }\r
2019   if(fdEdxTree) fdEdxTree->Print();\r
2020 \r
2021   //\r
2022   chain = new TChain("Laser");\r
2023   if(chain) { \r
2024     chain->Add("jotwinow_Temp_Trees.root");\r
2025     fLaserTree = chain->CopyTree("1");\r
2026     delete chain; chain=0; \r
2027   }\r
2028   if(fLaserTree) fLaserTree->Print();\r
2029 \r
2030   //\r
2031   chain = new TChain("MCEffTree");\r
2032   if(chain) { \r
2033     chain->Add("jotwinow_Temp_Trees.root");\r
2034     fMCEffTree = chain->CopyTree("1");\r
2035     delete chain; chain=0; \r
2036   }\r
2037   if(fMCEffTree) fMCEffTree->Print();\r
2038 \r
2039   //\r
2040   chain = new TChain("CosmicPairs");\r
2041   if(chain) { \r
2042     chain->Add("jotwinow_Temp_Trees.root");\r
2043     fCosmicPairsTree = chain->CopyTree("1");\r
2044     delete chain; chain=0; \r
2045   }\r
2046   if(fCosmicPairsTree) fCosmicPairsTree->Print();  \r
2047 \r
2048 \r
2049   OpenFile(1);\r
2050 \r
2051   // Post output data.\r
2052   PostData(1, fHighPtTree);\r
2053   PostData(2, fV0Tree);\r
2054   PostData(3, fdEdxTree);\r
2055   PostData(4, fLaserTree);\r
2056   PostData(5, fMCEffTree);\r
2057   PostData(6, fCosmicPairsTree);\r
2058 }\r
2059 \r
2060 //_____________________________________________________________________________\r
2061 void AliAnalysisTaskFilteredTree::Terminate(Option_t *) \r
2062 {\r
2063   // Called one at the end \r
2064   /*\r
2065   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));\r
2066   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
2067   TChain* chain = new TChain("highPt");\r
2068   if(!chain) return;\r
2069   chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");\r
2070   TTree *tree = chain->CopyTree("1");\r
2071   if (chain) { delete chain; chain=0; }\r
2072   if(!tree) return;\r
2073   tree->Print();\r
2074   fOutputSummary = tree;\r
2075 \r
2076   if (!fOutputSummary) {\r
2077     Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1));\r
2078     return;\r
2079   }\r
2080   */\r
2081 }\r
2082 \r
2083 //_____________________________________________________________________________\r
2084 Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)\r
2085 {\r
2086   //\r
2087   // calculate mc event true track multiplicity\r
2088   //\r
2089   if(!mcEvent) return 0;\r
2090 \r
2091   AliStack* stack = 0;\r
2092   Int_t mult = 0;\r
2093 \r
2094   // MC particle stack\r
2095   stack = mcEvent->Stack();\r
2096   if (!stack) return 0;\r
2097 \r
2098   //\r
2099   //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());\r
2100   //\r
2101 \r
2102   Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);\r
2103   if(!isEventOK) return 0; \r
2104 \r
2105   Int_t nPart  = stack->GetNtrack();\r
2106   for (Int_t iMc = 0; iMc < nPart; ++iMc) \r
2107   {\r
2108      TParticle* particle = stack->Particle(iMc);\r
2109      if (!particle)\r
2110      continue;\r
2111 \r
2112      // only charged particles\r
2113      if(!particle->GetPDG()) continue;\r
2114      Double_t charge = particle->GetPDG()->Charge()/3.;\r
2115      if (TMath::Abs(charge) < 0.001)\r
2116      continue;\r
2117       \r
2118      // physical primary\r
2119      Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
2120      if(!prim) continue;\r
2121 \r
2122      // checked accepted without pt cut\r
2123      //if(accCuts->AcceptTrack(particle)) \r
2124      if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) \r
2125      {\r
2126        mult++;\r
2127      }\r
2128   }\r
2129 \r
2130 return mult;  \r
2131 }\r
2132 \r
2133 \r