PWGPP/AliAnalysisTaskFilteredTree.cxx

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