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         "fileName.="<<&fileName<<\r
1195         "runNumber="<<runNumber<<\r
1196         "evtTimeStamp="<<evtTimeStamp<<\r
1197         "evtNumberInFile="<<evtNumberInFile<<\r
1198         "triggerClass="<<&triggerClass<<      //  trigger\r
1199         "Bz="<<bz<<\r
1200         "vtxESD.="<<pvtxESD<<\r
1201         //"vtxESDx="<<vtxX<<\r
1202         //"vtxESDy="<<vtxY<<\r
1203         //"vtxESDz="<<vtxZ<<\r
1204         "IRtot="<<ir1<<\r
1205         "IRint2="<<ir2<<\r
1206         "mult="<<mult<<\r
1207         "esdTrack.="<<ptrack<<\r
1208         "extTPCInnerC.="<<ptpcInnerC<<\r
1209         "extInnerParamC.="<<ptrackInnerC<<\r
1210         "extInnerParam.="<<ptrackInnerC2<<\r
1211         "extOuterITS.="<<pouterITSc<<\r
1212         "extInnerParamRef.="<<ptrackInnerC3<<\r
1213         "chi2TPCInnerC="<<chi2(0,0)<<\r
1214         "chi2InnerC="<<chi2trackC(0,0)<<\r
1215         "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
1216         "centralityF="<<centralityF;\r
1217         if (fUseMCInfo)\r
1218         {\r
1219           (*fTreeSRedirector)<<"highPt"<<\r
1220           "refTPCIn.="<<refTPCIn<<\r
1221           "refITS.="<<refITS<<\r
1222           "particle.="<<particle<<\r
1223           "particleMother.="<<particleMother<<\r
1224           "mech="<<mech<<\r
1225           "isPrim="<<isPrim<<\r
1226           "isFromStrangess="<<isFromStrangess<<\r
1227           "isFromConversion="<<isFromConversion<<\r
1228           "isFromMaterial="<<isFromMaterial<<\r
1229           "particleTPC.="<<particleTPC<<\r
1230           "particleMotherTPC.="<<particleMotherTPC<<\r
1231           "mechTPC="<<mechTPC<<\r
1232           "isPrimTPC="<<isPrimTPC<<\r
1233           "isFromStrangessTPC="<<isFromStrangessTPC<<\r
1234           "isFromConversionTPC="<<isFromConversionTPC<<\r
1235           "isFromMaterialTPC="<<isFromMaterialTPC<<\r
1236           "particleITS.="<<particleITS<<\r
1237           "particleMotherITS.="<<particleMotherITS<<\r
1238           "mechITS="<<mechITS<<\r
1239           "isPrimITS="<<isPrimITS<<\r
1240           "isFromStrangessITS="<<isFromStrangessITS<<\r
1241           "isFromConversionITS="<<isFromConversionITS<<\r
1242           "isFromMaterialITS="<<isFromMaterialITS;\r
1243         }\r
1244         //finish writing the entry\r
1245         (*fTreeSRedirector)<<"highPt"<<"\n";\r
1246       }\r
1247 \r
1248       delete pvtxESD;\r
1249       delete ptrack;\r
1250       delete ptpcInnerC;\r
1251       delete ptrackInnerC;\r
1252       delete ptrackInnerC2;\r
1253       delete pouterITSc;\r
1254       delete ptrackInnerC3;\r
1255 \r
1256       ////////////////////\r
1257       //delete the dummy objects we might have created.\r
1258       if (newvtxESD) {delete vtxESD; vtxESD=NULL;}\r
1259       if (newtrack) {delete track; track=NULL;}\r
1260       if (newtpcInnerC) {delete tpcInnerC; tpcInnerC=NULL;}\r
1261       if (newtrackInnerC) {delete trackInnerC; trackInnerC=NULL;}\r
1262       if (newtrackInnerC2) {delete trackInnerC2; trackInnerC2=NULL;}\r
1263       if (newouterITSc) {delete outerITSc; outerITSc=NULL;}\r
1264       if (newtrackInnerC3) {delete trackInnerC3; trackInnerC3=NULL;}\r
1265       if (newrefTPCIn) {delete refTPCIn; refTPCIn=NULL;}\r
1266       if (newrefITS) {delete refITS; refITS=NULL;}\r
1267       if (newparticle) {delete particle; particle=NULL;}\r
1268       if (newparticleMother) {delete particleMother; particleMother=NULL;}\r
1269       if (newparticleTPC) {delete particleTPC; particleTPC=NULL;}\r
1270       if (newparticleMotherTPC) {delete particleMotherTPC; particleMotherTPC=NULL;}\r
1271       if (newparticleITS) {delete particleITS; particleITS=NULL;}\r
1272       if (newparticleMotherITS) {delete particleMotherITS; particleMotherITS=NULL;}\r
1273       ///////////////\r
1274 \r
1275       delete tpcInnerC;\r
1276       delete trackInnerC;\r
1277       delete trackInnerC2;\r
1278       delete outerITSc;\r
1279       delete trackInnerC3;\r
1280     }\r
1281   }\r
1282 \r
1283 }\r
1284 \r
1285 \r
1286 //_____________________________________________________________________________\r
1287 void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1288 {\r
1289   //\r
1290   // Fill tree for efficiency studies MC only\r
1291 \r
1292   if(!esdEvent) {\r
1293     AliDebug(AliLog::kError, "esdEvent not available");\r
1294     return;\r
1295   }\r
1296 \r
1297    if(!mcEvent) {\r
1298     AliDebug(AliLog::kError, "mcEvent not available");\r
1299     return;\r
1300   }\r
1301 \r
1302   // get selection cuts\r
1303   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1304   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1305   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1306 \r
1307   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1308     AliDebug(AliLog::kError, "cuts not available");\r
1309     return;\r
1310   }\r
1311 \r
1312   // trigger selection\r
1313   Bool_t isEventTriggered = kTRUE;\r
1314   AliPhysicsSelection *physicsSelection = NULL;\r
1315   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1316 \r
1317   // \r
1318   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1319   if (!inputHandler)\r
1320   {\r
1321     Printf("ERROR: Could not receive input handler");\r
1322     return;\r
1323   }\r
1324    \r
1325   // get file name\r
1326   TTree *chain = (TChain*)GetInputData(0);\r
1327   if(!chain) { \r
1328     Printf("ERROR: Could not receive input chain");\r
1329     return;\r
1330   }\r
1331   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1332 \r
1333   // trigger\r
1334   if(evtCuts->IsTriggerRequired())  \r
1335   {\r
1336     // always MB\r
1337     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1338 \r
1339     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1340     if(!physicsSelection) return;\r
1341 \r
1342     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1343       // set trigger (V0AND)\r
1344       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1345       if(!triggerAnalysis) return;\r
1346       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1347     }\r
1348   }\r
1349 \r
1350   // centrality determination\r
1351   Float_t centralityF = -1;\r
1352   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1353   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1354 \r
1355   // use MC information\r
1356   AliHeader* header = 0;\r
1357   AliGenEventHeader* genHeader = 0;\r
1358   AliStack* stack = 0;\r
1359   TArrayF vtxMC(3);\r
1360 \r
1361   Int_t multMCTrueTracks = 0;\r
1362   if(IsUseMCInfo())\r
1363   {\r
1364     //\r
1365     if(!mcEvent) {\r
1366       AliDebug(AliLog::kError, "mcEvent not available");\r
1367       return;\r
1368     }\r
1369     // get MC event header\r
1370     header = mcEvent->Header();\r
1371     if (!header) {\r
1372       AliDebug(AliLog::kError, "Header not available");\r
1373       return;\r
1374     }\r
1375     // MC particle stack\r
1376     stack = mcEvent->Stack();\r
1377     if (!stack) {\r
1378       AliDebug(AliLog::kError, "Stack not available");\r
1379       return;\r
1380     }\r
1381 \r
1382     // get MC vertex\r
1383     genHeader = header->GenEventHeader();\r
1384     if (!genHeader) {\r
1385       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
1386       return;\r
1387     }\r
1388     genHeader->PrimaryVertex(vtxMC);\r
1389 \r
1390     // multipliticy of all MC primary tracks\r
1391     // in Zv, pt and eta ranges)\r
1392     multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
1393 \r
1394   } // end bUseMC\r
1395 \r
1396   // get reconstructed vertex  \r
1397   //const AliESDVertex* vtxESD = 0; \r
1398   AliESDVertex* vtxESD = 0; \r
1399   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1400         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1401   }\r
1402   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1403      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1404   }\r
1405   else {\r
1406         return;\r
1407   }\r
1408 \r
1409   if(!vtxESD) return;\r
1410 \r
1411   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1412   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1413   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1414 \r
1415   TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1416 \r
1417   // check event cuts\r
1418   if(isEventOK && isEventTriggered)\r
1419   {\r
1420     if(IsUseMCInfo()) \r
1421     {\r
1422       if(!stack) return;\r
1423 \r
1424       //\r
1425       // MC info\r
1426       //\r
1427       TParticle *particle=NULL;\r
1428       TParticle *particleMother=NULL;\r
1429       Int_t mech=-1;\r
1430 \r
1431       // reco event info\r
1432       Double_t vert[3] = {0}; \r
1433       vert[0] = vtxESD->GetXv();\r
1434       vert[1] = vtxESD->GetYv();\r
1435       vert[2] = vtxESD->GetZv();\r
1436       Int_t mult = vtxESD->GetNContributors();\r
1437       Double_t bz = esdEvent->GetMagneticField();\r
1438       Double_t runNumber = esdEvent->GetRunNumber();\r
1439       Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1440       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
1441 \r
1442       // loop over MC stack\r
1443       for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc) \r
1444       {\r
1445          particle = stack->Particle(iMc);\r
1446          if (!particle)\r
1447          continue;\r
1448 \r
1449          // only charged particles\r
1450          if(!particle->GetPDG()) continue;\r
1451          Double_t charge = particle->GetPDG()->Charge()/3.;\r
1452          if (TMath::Abs(charge) < 0.001)\r
1453          continue;\r
1454 \r
1455          // only primary particles\r
1456          Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
1457          if(!prim) continue;\r
1458 \r
1459          // downscale low-pT particles\r
1460          Double_t scalempt= TMath::Min(particle->Pt(),10.);\r
1461          Double_t downscaleF = gRandom->Rndm();\r
1462          downscaleF *= fLowPtTrackDownscaligF;\r
1463          if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
1464 \r
1465          // is particle in acceptance\r
1466          if(!accCuts->AcceptTrack(particle)) continue;\r
1467        \r
1468          // check if particle reconstructed\r
1469          Bool_t isRec = kFALSE;\r
1470          Int_t  trackIndex = -1;\r
1471          for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
1472          {\r
1473            \r
1474            AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
1475            if(!track) continue;\r
1476            if(track->Charge()==0) continue;\r
1477            if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track)) \r
1478            {\r
1479              Int_t label =  TMath::Abs(track->GetLabel());\r
1480              if(label == iMc) {\r
1481                isRec = kTRUE;\r
1482                trackIndex = iTrack;\r
1483                break;\r
1484              }\r
1485            } \r
1486          }\r
1487 \r
1488          // Store information in the output tree\r
1489          AliESDtrack *recTrack = NULL; \r
1490          if(trackIndex>-1)  { \r
1491            recTrack = esdEvent->GetTrack(trackIndex); \r
1492          } else {\r
1493            recTrack = new AliESDtrack(); \r
1494          } \r
1495 \r
1496          particleMother = GetMother(particle,stack);\r
1497          mech = particle->GetUniqueID();\r
1498 \r
1499          //MC particle track length\r
1500          Double_t tpcTrackLength = 0.;\r
1501          AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);\r
1502          if(mcParticle) {\r
1503            Int_t counter;\r
1504            tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);\r
1505          } \r
1506 \r
1507 \r
1508          //\r
1509          if(fTreeSRedirector) {\r
1510            (*fTreeSRedirector)<<"MCEffTree"<<\r
1511            "fileName.="<<&fileName<<\r
1512            "triggerClass.="<<&triggerClass<<\r
1513            "runNumber="<<runNumber<<\r
1514            "evtTimeStamp="<<evtTimeStamp<<\r
1515            "evtNumberInFile="<<evtNumberInFile<<\r
1516            "Bz="<<bz<<\r
1517            "vtxESD.="<<vtxESD<<\r
1518            "mult="<<mult<<\r
1519            "esdTrack.="<<recTrack<<\r
1520            "isRec="<<isRec<<\r
1521            "tpcTrackLength="<<tpcTrackLength<<\r
1522            "particle.="<<particle<<\r
1523            "particleMother.="<<particleMother<<\r
1524            "mech="<<mech<<\r
1525            "\n";\r
1526          }\r
1527 \r
1528          if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;\r
1529       }\r
1530     }\r
1531   }\r
1532   \r
1533 }\r
1534 \r
1535 //_____________________________________________________________________________\r
1536 Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {\r
1537   //\r
1538   // check if particle is Z > 1 \r
1539   //\r
1540   if (track->GetTPCNcls() < 60) return kFALSE;\r
1541   Double_t mom = track->GetInnerParam()->GetP();\r
1542   if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization\r
1543   Float_t dca[2], bCov[3];\r
1544   track->GetImpactParameters(dca,bCov);\r
1545   //\r
1546 \r
1547   Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);\r
1548 \r
1549   if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;\r
1550 \r
1551   return kFALSE;\r
1552 }\r
1553 \r
1554 //_____________________________________________________________________________\r
1555 void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1556 {\r
1557   //\r
1558   // Select real events with V0 (K0s and Lambda) high-pT candidates\r
1559   //\r
1560   if(!esdEvent) {\r
1561     AliDebug(AliLog::kError, "esdEvent not available");\r
1562     return;\r
1563   }\r
1564 \r
1565   // get selection cuts\r
1566   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1567   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1568   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1569 \r
1570   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1571     AliDebug(AliLog::kError, "cuts not available");\r
1572     return;\r
1573   }\r
1574 \r
1575   // trigger selection\r
1576   Bool_t isEventTriggered = kTRUE;\r
1577   AliPhysicsSelection *physicsSelection = NULL;\r
1578   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1579 \r
1580   // \r
1581   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1582   if (!inputHandler)\r
1583   {\r
1584     Printf("ERROR: Could not receive input handler");\r
1585     return;\r
1586   }\r
1587    \r
1588   // get file name\r
1589   TTree *chain = (TChain*)GetInputData(0);\r
1590   if(!chain) { \r
1591     Printf("ERROR: Could not receive input chain");\r
1592     return;\r
1593   }\r
1594   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1595 \r
1596   // trigger\r
1597   if(evtCuts->IsTriggerRequired())  \r
1598   {\r
1599     // always MB\r
1600     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1601 \r
1602     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1603     if(!physicsSelection) return;\r
1604     //SetPhysicsTriggerSelection(physicsSelection);\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   // centrality determination\r
1615   Float_t centralityF = -1;\r
1616   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1617   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1618 \r
1619 \r
1620   // get reconstructed vertex  \r
1621   //const AliESDVertex* vtxESD = 0; \r
1622   AliESDVertex* vtxESD = 0; \r
1623   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1624         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1625   }\r
1626   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1627      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1628   }\r
1629   else {\r
1630         return;\r
1631   }\r
1632 \r
1633   if(!vtxESD) return;\r
1634 \r
1635   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1636   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1637   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1638 \r
1639   // check event cuts\r
1640   if(isEventOK && isEventTriggered) {\r
1641   //\r
1642   // Dump the pt downscaled V0 into the tree\r
1643   // \r
1644   Int_t ntracks = esdEvent->GetNumberOfTracks();\r
1645   Int_t nV0s = esdEvent->GetNumberOfV0s();\r
1646   Int_t run = esdEvent->GetRunNumber();\r
1647   Int_t time= esdEvent->GetTimeStamp();\r
1648   Int_t evNr=esdEvent->GetEventNumberInFile();\r
1649   Double_t bz = esdEvent->GetMagneticField();\r
1650 \r
1651 \r
1652   for (Int_t iv0=0; iv0<nV0s; iv0++){\r
1653     AliESDv0 * v0 = esdEvent->GetV0(iv0);\r
1654     if (!v0) continue;\r
1655     AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));\r
1656     AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));\r
1657     if (!track0) continue;\r
1658     if (!track1) continue;\r
1659     if (track0->GetSign()<0) {\r
1660       track1 = esdEvent->GetTrack(v0->GetIndex(0));\r
1661       track0 = esdEvent->GetTrack(v0->GetIndex(1));\r
1662     }\r
1663     //\r
1664     Bool_t isDownscaled = IsV0Downscaled(v0);\r
1665     if (isDownscaled) continue;\r
1666     AliKFParticle kfparticle; //\r
1667     Int_t type=GetKFParticle(v0,esdEvent,kfparticle);\r
1668     if (type==0) continue;   \r
1669     TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1670 \r
1671     if(!fTreeSRedirector) return;\r
1672     (*fTreeSRedirector)<<"V0s"<<\r
1673       "isDownscaled="<<isDownscaled<<\r
1674       "triggerClass="<<&triggerClass<<      //  trigger\r
1675       "Bz="<<bz<<\r
1676       "fileName.="<<&fileName<<\r
1677       "runNumber="<<run<<\r
1678       "evtTimeStamp="<<time<<\r
1679       "evtNumberInFile="<<evNr<<\r
1680       "type="<<type<<\r
1681       "ntracks="<<ntracks<<\r
1682       "v0.="<<v0<<\r
1683       "kf.="<<&kfparticle<<\r
1684       "track0.="<<track0<<\r
1685       "track1.="<<track1<<\r
1686       "centralityF="<<centralityF<<\r
1687       "\n";\r
1688   }\r
1689   }\r
1690 }\r
1691 \r
1692 //_____________________________________________________________________________\r
1693 void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1694 {\r
1695   //\r
1696   // Select real events with large TPC dEdx signal\r
1697   //\r
1698   if(!esdEvent) {\r
1699     AliDebug(AliLog::kError, "esdEvent not available");\r
1700     return;\r
1701   }\r
1702 \r
1703   // get selection cuts\r
1704   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
1705   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1706   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1707 \r
1708   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
1709     AliDebug(AliLog::kError, "cuts not available");\r
1710     return;\r
1711   }\r
1712 \r
1713   // get file name\r
1714   TTree *chain = (TChain*)GetInputData(0);\r
1715   if(!chain) { \r
1716     Printf("ERROR: Could not receive input chain");\r
1717     return;\r
1718   }\r
1719   TObjString fileName(chain->GetCurrentFile()->GetName());\r
1720 \r
1721   // trigger\r
1722   Bool_t isEventTriggered = kTRUE;\r
1723   AliPhysicsSelection *physicsSelection = NULL;\r
1724   AliTriggerAnalysis* triggerAnalysis = NULL;\r
1725 \r
1726   // \r
1727   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1728   if (!inputHandler)\r
1729   {\r
1730     Printf("ERROR: Could not receive input handler");\r
1731     return;\r
1732   }\r
1733 \r
1734   if(evtCuts->IsTriggerRequired())  \r
1735   {\r
1736     // always MB\r
1737     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1738 \r
1739     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1740     if(!physicsSelection) return;\r
1741 \r
1742     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1743       // set trigger (V0AND)\r
1744       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1745       if(!triggerAnalysis) return;\r
1746       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1747     }\r
1748   }\r
1749 \r
1750   // get reconstructed vertex  \r
1751   AliESDVertex* vtxESD = 0; \r
1752   if(GetAnalysisMode() == kTPCAnalysisMode) {\r
1753         vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
1754   }\r
1755   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
1756      vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
1757   }\r
1758   else {\r
1759         return;\r
1760   }\r
1761   if(!vtxESD) return;\r
1762 \r
1763   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1764   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1765   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1766 \r
1767 \r
1768   // check event cuts\r
1769   if(isEventOK && isEventTriggered)\r
1770   {\r
1771     Double_t vert[3] = {0}; \r
1772     vert[0] = vtxESD->GetXv();\r
1773     vert[1] = vtxESD->GetYv();\r
1774     vert[2] = vtxESD->GetZv();\r
1775     Int_t mult = vtxESD->GetNContributors();\r
1776     Double_t bz = esdEvent->GetMagneticField();\r
1777     Double_t runNumber = esdEvent->GetRunNumber();\r
1778     Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
1779     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
1780 \r
1781     // large dEdx\r
1782     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
1783     {\r
1784       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
1785       if(!track) continue;\r
1786       if(track->Charge()==0) continue;\r
1787       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
1788       if(!accCuts->AcceptTrack(track)) continue;\r
1789 \r
1790       if(!IsHighDeDxParticle(track)) continue;\r
1791       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
1792 \r
1793       if(!fTreeSRedirector) return;\r
1794       (*fTreeSRedirector)<<"dEdx"<<\r
1795       "fileName.="<<&fileName<<\r
1796       "runNumber="<<runNumber<<\r
1797       "evtTimeStamp="<<evtTimeStamp<<\r
1798       "evtNumberInFile="<<evtNumberInFile<<\r
1799         "triggerClass="<<&triggerClass<<      //  trigger\r
1800       "Bz="<<bz<<\r
1801       "vtxESD.="<<vtxESD<<\r
1802       "mult="<<mult<<\r
1803       "esdTrack.="<<track<<\r
1804       "\n";\r
1805     }\r
1806   }\r
1807 }\r
1808 \r
1809 //_____________________________________________________________________________\r
1810 Int_t   AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)\r
1811 {\r
1812   //\r
1813   // Create KF particle in case the V0 fullfill selection criteria\r
1814   //\r
1815   // Selection criteria\r
1816   //  0. algorithm cut\r
1817   //  1. track cut\r
1818   //  3. chi2 cut\r
1819   //  4. rough mass cut\r
1820   //  5. Normalized pointing angle cut\r
1821   //\r
1822   const Double_t cutMass=0.2;\r
1823   const Double_t kSigmaDCACut=3;\r
1824   //\r
1825   // 0.) algo cut - accept only on the fly\r
1826   //\r
1827   if (v0->GetOnFlyStatus() ==kFALSE) return 0;     \r
1828   //\r
1829   // 1.) track cut\r
1830   // \r
1831   AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));\r
1832   AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));\r
1833   /*\r
1834     TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";\r
1835     TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";\r
1836     TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";\r
1837   */  \r
1838   if (TMath::Abs(track0->GetTgl())>1) return 0;\r
1839   if (TMath::Abs(track1->GetTgl())>1) return 0;\r
1840   if ((track0->GetTPCClusterInfo(2,1))<100) return 0;\r
1841   if ((track1->GetTPCClusterInfo(2,1))<100) return 0;\r
1842   //if ((track0->GetITSclusters(0))<2) return 0;\r
1843   //if ((track1->GetITSclusters(0))<2) return 0; \r
1844   Float_t pos0[2]={0}, cov0[3]={0};\r
1845   Float_t pos1[2]={0}, cov1[3]={0};\r
1846   track0->GetImpactParameters(pos0,cov0);\r
1847   track0->GetImpactParameters(pos1,cov1);\r
1848   //\r
1849   if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;\r
1850   if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;\r
1851   // \r
1852   //\r
1853   // 3.) Chi2 cut\r
1854   //\r
1855   Double_t chi2KF = v0->GetKFInfo(2,2,2);\r
1856   if (chi2KF>25) return 0;\r
1857   //\r
1858   // 4.) Rough mass cut - 0.200 GeV\r
1859   //\r
1860   static Double_t masses[2]={-1};\r
1861   if (masses[0]<0){\r
1862     masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();\r
1863     masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();\r
1864   }\r
1865   Double_t mass00=  v0->GetEffMass(0,0);\r
1866   Double_t mass22=  v0->GetEffMass(2,2);\r
1867   Double_t mass42=  v0->GetEffMass(4,2);\r
1868   Double_t mass24=  v0->GetEffMass(2,4);\r
1869   Bool_t massOK=kFALSE;\r
1870   Int_t type=0;\r
1871   Int_t ptype=0;\r
1872   Double_t dmass=1;\r
1873   Int_t p1=0, p2=0;\r
1874   if (TMath::Abs(mass00-0)<cutMass) {\r
1875     massOK=kTRUE; type+=1; \r
1876     if (TMath::Abs(mass00-0)<dmass) {\r
1877       ptype=1;\r
1878       dmass=TMath::Abs(mass00-0);      \r
1879       p1=0; p2=0;\r
1880     } \r
1881   }\r
1882   if (TMath::Abs(mass24-masses[1])<cutMass) {\r
1883     massOK=kTRUE; type+=2; \r
1884     if (TMath::Abs(mass24-masses[1])<dmass){\r
1885       dmass = TMath::Abs(mass24-masses[1]);\r
1886       ptype=2;\r
1887       p1=2; p2=4;\r
1888     }\r
1889   }\r
1890   if (TMath::Abs(mass42-masses[1])<cutMass) {\r
1891     massOK=kTRUE; type+=4;\r
1892     if (TMath::Abs(mass42-masses[1])<dmass){\r
1893       dmass = TMath::Abs(mass42-masses[1]);\r
1894       ptype=4;\r
1895       p1=4; p2=2;\r
1896     }\r
1897   }\r
1898   if (TMath::Abs(mass22-masses[0])<cutMass) {\r
1899     massOK=kTRUE; type+=8;\r
1900     if (TMath::Abs(mass22-masses[0])<dmass){\r
1901       dmass = TMath::Abs(mass22-masses[0]);\r
1902       ptype=8;\r
1903       p1=2; p2=2;\r
1904     }\r
1905   }\r
1906   if (type==0) return 0;\r
1907   //\r
1908   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};\r
1909   const AliExternalTrackParam *paramP = v0->GetParamP();\r
1910   const AliExternalTrackParam *paramN = v0->GetParamN();\r
1911   if (paramP->GetSign()<0){\r
1912     paramP=v0->GetParamP();\r
1913     paramN=v0->GetParamN();\r
1914   }\r
1915   //Double_t *pparam1 = (Double_t*)paramP->GetParameter();\r
1916   //Double_t *pparam2 = (Double_t*)paramN->GetParameter();\r
1917   //\r
1918   AliKFParticle kfp1( *paramP, spdg[p1]  );\r
1919   AliKFParticle kfp2( *paramN, -1 *spdg[p2]  );\r
1920   AliKFParticle V0KF;\r
1921   (V0KF)+=kfp1;\r
1922   (V0KF)+=kfp2;\r
1923   kfparticle=V0KF;\r
1924   //\r
1925   // Pointing angle\r
1926   //\r
1927   Double_t  errPhi    = V0KF.GetErrPhi();\r
1928   Double_t  pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());\r
1929   if (pointAngle/errPhi>10) return 0;  \r
1930   //\r
1931   return ptype;  \r
1932 }\r
1933 \r
1934 //_____________________________________________________________________________\r
1935 Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)\r
1936 {\r
1937   //\r
1938   // Downscale randomly low pt V0\r
1939   //\r
1940   //return kFALSE;\r
1941   Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());\r
1942   Double_t scalempt= TMath::Min(maxPt,10.);\r
1943   Double_t downscaleF = gRandom->Rndm();\r
1944   downscaleF *= fLowPtV0DownscaligF;\r
1945   \r
1946   //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
1947   if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;\r
1948   return kFALSE;\r
1949 \r
1950   /*\r
1951     TH1F his1("his1","his1",100,0,10);\r
1952     TH1F his2("his2","his2",100,0,10);\r
1953     {for (Int_t i=0; i<10000; i++){\r
1954        Double_t rnd=gRandom->Exp(1);\r
1955        Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();\r
1956        his1->Fill(rnd); \r
1957        if (!isDownscaled) his2->Fill(rnd); \r
1958     }}\r
1959 \r
1960    */\r
1961 \r
1962 }\r
1963 \r
1964 \r
1965 \r
1966 //_____________________________________________________________________________\r
1967 Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
1968 {\r
1969  // Constrain TPC inner params constrained\r
1970  //\r
1971       if(!tpcInnerC) return kFALSE; \r
1972       if(!vtx) return kFALSE;\r
1973 \r
1974       Double_t dz[2],cov[3];\r
1975       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1976       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1977       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1978       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1979 \r
1980 \r
1981       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1982       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
1983       Double_t c[3]={covar[2],0.,covar[5]};\r
1984       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
1985       if (chi2C>kVeryBig) return kFALSE; \r
1986 \r
1987       if(!tpcInnerC->Update(p,c)) return kFALSE;\r
1988 \r
1989   return kTRUE;\r
1990 }\r
1991 \r
1992 //_____________________________________________________________________________\r
1993 Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
1994 {\r
1995  // Constrain TPC inner params constrained\r
1996  //\r
1997       if(!trackInnerC) return kFALSE; \r
1998       if(!vtx) return kFALSE;\r
1999 \r
2000       const Double_t kRadius  = 2.8; \r
2001       const Double_t kMaxStep = 1.0; \r
2002 \r
2003       Double_t dz[2],cov[3];\r
2004 \r
2005       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
2006       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
2007       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
2008 \r
2009       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
2010       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
2011 \r
2012       //\r
2013       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
2014       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
2015       Double_t c[3]={covar[2],0.,covar[5]};\r
2016       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
2017       if (chi2C>kVeryBig) return kFALSE; \r
2018 \r
2019       if(!trackInnerC->Update(p,c)) return kFALSE;\r
2020 \r
2021   return kTRUE;\r
2022 }\r
2023 \r
2024 \r
2025 //_____________________________________________________________________________\r
2026 TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack) \r
2027 {\r
2028   if(!particle) return NULL;\r
2029   if(!stack) return NULL;\r
2030 \r
2031   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
2032   TParticle* mother = NULL; \r
2033   mother = stack->Particle(motherLabel); \r
2034 \r
2035 return mother;\r
2036 }\r
2037 \r
2038 //_____________________________________________________________________________\r
2039 Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack) \r
2040 {\r
2041   Bool_t isFromConversion = kFALSE;\r
2042 \r
2043   if(stack) {\r
2044     TParticle* particle = stack->Particle(label);\r
2045 \r
2046     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2047     {\r
2048        Int_t mech = particle->GetUniqueID(); // production mechanism \r
2049        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2050 \r
2051        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
2052        TParticle* mother = stack->Particle(motherLabel); \r
2053        if(mother) {\r
2054           Int_t motherPdg = mother->GetPdgCode();\r
2055 \r
2056           if(!isPrim && mech==5 && motherPdg==kGamma) { \r
2057             isFromConversion=kTRUE; \r
2058           }\r
2059        }\r
2060     } \r
2061   } \r
2062 \r
2063   return isFromConversion;\r
2064 }\r
2065 \r
2066 //_____________________________________________________________________________\r
2067 Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack) \r
2068 {\r
2069   Bool_t isFromMaterial = kFALSE;\r
2070 \r
2071   if(stack) {\r
2072     TParticle* particle = stack->Particle(label);\r
2073 \r
2074     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2075     {\r
2076        Int_t mech = particle->GetUniqueID(); // production mechanism \r
2077        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2078 \r
2079        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
2080        TParticle* mother = stack->Particle(motherLabel); \r
2081        if(mother) {\r
2082           if(!isPrim && mech==13) { \r
2083             isFromMaterial=kTRUE; \r
2084           }\r
2085        }\r
2086      } \r
2087   } \r
2088 \r
2089   return isFromMaterial;\r
2090 }\r
2091 \r
2092 //_____________________________________________________________________________\r
2093 Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
2094 {\r
2095   Bool_t isFromStrangeness = kFALSE;\r
2096 \r
2097   if(stack) {\r
2098     TParticle* particle = stack->Particle(label);\r
2099 \r
2100     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
2101     {\r
2102        Int_t mech = particle->GetUniqueID(); // production mechanism \r
2103        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
2104 \r
2105        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
2106        TParticle* mother = stack->Particle(motherLabel); \r
2107        if(mother) {\r
2108           Int_t motherPdg = mother->GetPdgCode();\r
2109 \r
2110           // K+-, lambda, antilambda, K0s decays\r
2111           if(!isPrim && mech==4 && \r
2112               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
2113           {\r
2114             isFromStrangeness = kTRUE;\r
2115           } \r
2116        }\r
2117     } \r
2118   } \r
2119 \r
2120   return isFromStrangeness;\r
2121 }\r
2122 \r
2123 \r
2124 //_____________________________________________________________________________\r
2125 void AliAnalysisTaskFilteredTree::FinishTaskOutput() \r
2126 {\r
2127   //\r
2128   // Called one at the end \r
2129   // locally on working node\r
2130   //\r
2131 \r
2132   //// must be deleted to store trees\r
2133   //if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;\r
2134 \r
2135   //// open temporary file and copy trees to the ouptut container\r
2136 \r
2137   //TChain* chain = 0;\r
2138   ////\r
2139   //chain = new TChain("highPt");\r
2140   //if(chain) { \r
2141   //  chain->Add("jotwinow_Temp_Trees.root");\r
2142   //  fHighPtTree = chain->CopyTree("1");\r
2143   //  delete chain; chain=0; \r
2144   //}\r
2145   //if(fHighPtTree) fHighPtTree->Print();\r
2146 \r
2147   ////\r
2148   //chain = new TChain("V0s");\r
2149   //if(chain) { \r
2150   //  chain->Add("jotwinow_Temp_Trees.root");\r
2151   //  fV0Tree = chain->CopyTree("1");\r
2152   //  delete chain; chain=0; \r
2153   //}\r
2154   //if(fV0Tree) fV0Tree->Print();\r
2155 \r
2156   ////\r
2157   //chain = new TChain("dEdx");\r
2158   //if(chain) { \r
2159   //  chain->Add("jotwinow_Temp_Trees.root");\r
2160   //  fdEdxTree = chain->CopyTree("1");\r
2161   //  delete chain; chain=0; \r
2162   //}\r
2163   //if(fdEdxTree) fdEdxTree->Print();\r
2164 \r
2165   ////\r
2166   //chain = new TChain("Laser");\r
2167   //if(chain) { \r
2168   //  chain->Add("jotwinow_Temp_Trees.root");\r
2169   //  fLaserTree = chain->CopyTree("1");\r
2170   //  delete chain; chain=0; \r
2171   //}\r
2172   //if(fLaserTree) fLaserTree->Print();\r
2173 \r
2174   ////\r
2175   //chain = new TChain("MCEffTree");\r
2176   //if(chain) { \r
2177   //  chain->Add("jotwinow_Temp_Trees.root");\r
2178   //  fMCEffTree = chain->CopyTree("1");\r
2179   //  delete chain; chain=0; \r
2180   //}\r
2181   //if(fMCEffTree) fMCEffTree->Print();\r
2182 \r
2183   ////\r
2184   //chain = new TChain("CosmicPairs");\r
2185   //if(chain) { \r
2186   //  chain->Add("jotwinow_Temp_Trees.root");\r
2187   //  fCosmicPairsTree = chain->CopyTree("1");\r
2188   //  delete chain; chain=0; \r
2189   //}\r
2190   //if(fCosmicPairsTree) fCosmicPairsTree->Print();  \r
2191 \r
2192 \r
2193   //OpenFile(1);\r
2194 \r
2195   // Post output data.\r
2196   //PostData(1, fHighPtTree);\r
2197   //PostData(2, fV0Tree);\r
2198   //PostData(3, fdEdxTree);\r
2199   //PostData(4, fLaserTree);\r
2200   //PostData(5, fMCEffTree);\r
2201   //PostData(6, fCosmicPairsTree);\r
2202 }\r
2203 \r
2204 //_____________________________________________________________________________\r
2205 void AliAnalysisTaskFilteredTree::Terminate(Option_t *) \r
2206 {\r
2207   // Called one at the end \r
2208   /*\r
2209   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));\r
2210   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
2211   TChain* chain = new TChain("highPt");\r
2212   if(!chain) return;\r
2213   chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");\r
2214   TTree *tree = chain->CopyTree("1");\r
2215   if (chain) { delete chain; chain=0; }\r
2216   if(!tree) return;\r
2217   tree->Print();\r
2218   fOutputSummary = tree;\r
2219 \r
2220   if (!fOutputSummary) {\r
2221     Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1));\r
2222     return;\r
2223   }\r
2224   */\r
2225 }\r
2226 \r
2227 //_____________________________________________________________________________\r
2228 Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)\r
2229 {\r
2230   //\r
2231   // calculate mc event true track multiplicity\r
2232   //\r
2233   if(!mcEvent) return 0;\r
2234 \r
2235   AliStack* stack = 0;\r
2236   Int_t mult = 0;\r
2237 \r
2238   // MC particle stack\r
2239   stack = mcEvent->Stack();\r
2240   if (!stack) return 0;\r
2241 \r
2242   //\r
2243   //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());\r
2244   //\r
2245 \r
2246   Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);\r
2247   if(!isEventOK) return 0; \r
2248 \r
2249   Int_t nPart  = stack->GetNtrack();\r
2250   for (Int_t iMc = 0; iMc < nPart; ++iMc) \r
2251   {\r
2252      TParticle* particle = stack->Particle(iMc);\r
2253      if (!particle)\r
2254      continue;\r
2255 \r
2256      // only charged particles\r
2257      if(!particle->GetPDG()) continue;\r
2258      Double_t charge = particle->GetPDG()->Charge()/3.;\r
2259      if (TMath::Abs(charge) < 0.001)\r
2260      continue;\r
2261       \r
2262      // physical primary\r
2263      Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
2264      if(!prim) continue;\r
2265 \r
2266      // checked accepted without pt cut\r
2267      //if(accCuts->AcceptTrack(particle)) \r
2268      if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) \r
2269      {\r
2270        mult++;\r
2271      }\r
2272   }\r
2273 \r
2274 return mult;  \r
2275 }\r
2276 \r
2277 \r