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