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