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