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