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