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