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;
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 AliExternalTrackParam dummyexternaltrackparam;
1402         //if (!dummyexternaltrackparam) dummyexternaltrackparam=new AliExternalTrackParam();
1403         static AliTrackReference dummytrackreference;
1404         //if (!dummytrackreference) dummytrackreference=new AliTrackReference();
1405         static TParticle dummyparticle;
1406         //if (!dummyparticle) dummyparticle=new TParticle();
1407
1408         //assign the dummy objects if needed
1409         if (!track) {track=fDummyTrack;}
1410         if (!friendTrack) {friendTrack=fDummyFriendTrack;}
1411         if (!vtxESD) {vtxESD=&dummyvtxESD;}
1412         if (mcEvent)
1413         {
1414           if (!refTPCIn) {refTPCIn=&dummytrackreference;}
1415           if (!refITS) {refITS=&dummytrackreference;}
1416           if (!particle) {particle=&dummyparticle;}
1417           if (!particleMother) {particleMother=&dummyparticle;}
1418           if (!particleTPC) {particleTPC=&dummyparticle;}
1419           if (!particleMotherTPC) {particleMotherTPC=&dummyparticle;}
1420           if (!particleITS) {particleITS=&dummyparticle;}
1421           if (!particleMotherITS) {particleMotherITS=&dummyparticle;}
1422         }
1423         //the following are guaranteed to exist:
1424         //if (!tpcInnerC) {tpcInnerC=&dummyexternaltrackparam;}
1425         //if (!trackInnerC) {trackInnerC=&dummyexternaltrackparam;}
1426         //if (!trackInnerC2) {trackInnerC2=&dummyexternaltrackparam;}
1427         //if (!outerITSc) {outerITSc=&dummyexternaltrackparam;}
1428         //if (!trackInnerC3) {trackInnerC3=&dummyexternaltrackparam;}
1429         /////////////////////////
1430
1431         //Double_t vtxX=vtxESD->GetX();
1432         //Double_t vtxY=vtxESD->GetY();
1433         //Double_t vtxZ=vtxESD->GetZ();
1434
1435         //AliESDVertex* pvtxESD = (AliESDVertex*)vtxESD->Clone();
1436         //AliESDtrack* ptrack=(AliESDtrack*)track->Clone();
1437         //AliExternalTrackParam* ptpcInnerC = (AliExternalTrackParam*)tpcInnerC->Clone();
1438         //AliExternalTrackParam* ptrackInnerC = (AliExternalTrackParam*)trackInnerC->Clone();
1439         //AliExternalTrackParam* ptrackInnerC2 = (AliExternalTrackParam*)trackInnerC2->Clone();
1440         //AliExternalTrackParam* pouterITSc = (AliExternalTrackParam*)outerITSc->Clone();
1441         //AliExternalTrackParam* ptrackInnerC3 = (AliExternalTrackParam*)trackInnerC3->Clone();
1442         //AliESDVertex* pvtxESD = new AliESDVertex(*vtxESD);
1443         //AliESDtrack* ptrack= new AliESDtrack(*track);
1444         //AliExternalTrackParam* ptpcInnerC = new AliExternalTrackParam(*tpcInnerC);
1445         //AliExternalTrackParam* ptrackInnerC = new AliExternalTrackParam(*trackInnerC);
1446         //AliExternalTrackParam* ptrackInnerC2 = new AliExternalTrackParam(*trackInnerC2);
1447         //AliExternalTrackParam* pouterITSc =  new AliExternalTrackParam(*outerITSc);
1448         //AliExternalTrackParam* ptrackInnerC3 = new AliExternalTrackParam(*trackInnerC3);
1449
1450         Int_t ntracks = esdEvent->GetNumberOfTracks();
1451
1452         // Global event id calculation using orbitID, bunchCrossingID and periodID
1453         ULong64_t orbitID      = (ULong64_t)esdEvent->GetOrbitNumber();
1454         ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1455         ULong64_t periodID     = (ULong64_t)esdEvent->GetPeriodNumber();
1456         ULong64_t gid          = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1457
1458         // fill histograms
1459         FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0));
1460
1461         if(fTreeSRedirector && dumpToTree && fFillTree)
1462         {
1463
1464           //if (friendTrack)
1465           //{
1466           //  const AliTrackPointArray* array = friendTrack->GetTrackPointArray();
1467           //  if (!array) {printf("we have a friend, but the ponits are empty\n"); continue;}
1468           //  if (friendTrack==fDummyFriendTrack) printf("using the dummy friend track\n");
1469           //  cout<<array->GetX()[0]<<" "<<array->GetX()[2]<<endl;
1470           //}
1471           //else
1472           //{
1473           //  printf("no friend track\n");
1474           //}
1475
1476
1477           (*fTreeSRedirector)<<"highPt"<<
1478             "gid="<<gid<<
1479             "fileName.="<<&fileName<<                // name of the chunk file (hopefully full)
1480             "runNumber="<<runNumber<<                // runNumber
1481             "evtTimeStamp="<<evtTimeStamp<<          // time stamp of event (in seconds)
1482             "evtNumberInFile="<<evtNumberInFile<<    // event number
1483             "triggerClass="<<&triggerClass<<         // trigger class as a string
1484             "Bz="<<bz<<                              // solenoid magnetic field in the z direction (in kGaus)
1485             "vtxESD.="<<vtxESD<<                    // vertexer ESD tracks (can be biased by TPC pileup tracks)
1486             //"vtxESDx="<<vtxX<<
1487             //"vtxESDy="<<vtxY<<
1488             //"vtxESDz="<<vtxZ<<
1489             "IRtot="<<ir1<<                         // interaction record (trigger) counters - coutner 1
1490             "IRint2="<<ir2<<                        // interaction record (trigger) coutners - counter 2
1491             "mult="<<mult<<                         // multiplicity of tracks pointing to the primary vertex
1492             "ntracks="<<ntracks<<                   // number of the esd tracks (to take into account the pileup in the TPC)
1493             //                                           important variables for the pile-up studies
1494             "contTPC="<< contTPC<<                    // number of contributors to the TPC primary vertex candidate
1495             "contSPD="<< contSPD<<                    // number of contributors to the SPD primary vertex candidate
1496             "vertexPosTPC.="<<&vertexPosTPC<<          // TPC vertex position
1497             "vertexPosSPD.="<<&vertexPosSPD<<          // SPD vertex position
1498             "ntracksTPC="<<ntracksTPC<<               // total number of the TPC tracks which were refitted
1499             "ntracksITS="<<ntracksITS<<               // total number of the ITS tracks which were refitted
1500             //
1501             "esdTrack.="<<track<<                  // esdTrack as used in the physical analysis
1502             "friendTrack.="<<friendTrack<<      // esdFriendTrack associated to the esdTrack
1503             "extTPCInnerC.="<<tpcInnerC<<          // ???
1504             "extInnerParamC.="<<trackInnerC<<      // ???
1505             "extInnerParam.="<<trackInnerC2<<      // ???
1506             "extOuterITS.="<<outerITSc<<           // ???
1507             "extInnerParamRef.="<<trackInnerC3<<   // ???
1508             "chi2TPCInnerC="<<chi2(0,0)<<           // chi2   of tracks ???
1509             "chi2InnerC="<<chi2trackC(0,0)<<        // chi2s  of tracks TPCinner to the combined
1510             "chi2OuterITS="<<chi2OuterITS(0,0)<<    // chi2s  of tracks TPC at inner wall to the ITSout
1511             "centralityF="<<centralityF;
1512           if (mcEvent)
1513           {
1514             static AliTrackReference refDummy;
1515             if (!refITS) refITS = &refDummy;
1516             if (!refTRD) refTRD = &refDummy;
1517             if (!refTOF) refTOF = &refDummy;
1518             if (!refEMCAL) refEMCAL = &refDummy;
1519             if (!refPHOS) refPHOS = &refDummy;
1520             (*fTreeSRedirector)<<"highPt"<<
1521               "multMCTrueTracks="<<multMCTrueTracks<<   // mC track multiplicities
1522               "nrefITS="<<nrefITS<<              // number of track references in the ITS
1523               "nrefTPC="<<nrefTPC<<              // number of track references in the TPC
1524               "nrefTRD="<<nrefTRD<<              // number of track references in the TRD
1525               "nrefTOF="<<nrefTOF<<              // number of track references in the TOF
1526               "nrefEMCAL="<<nrefEMCAL<<              // number of track references in the TOF
1527               "nrefPHOS="<<nrefPHOS<<              // number of track references in the TOF
1528               "refTPCIn.="<<refTPCIn<<
1529               "refTPCOut.="<<refTPCOut<<
1530               "refITS.="<<refITS<<
1531               "refTRD.="<<refTRD<<
1532               "refTOF.="<<refTOF<<
1533               "refEMCAL.="<<refEMCAL<<
1534               "refPHOS.="<<refPHOS<<
1535               "particle.="<<particle<<
1536               "particleMother.="<<particleMother<<
1537               "mech="<<mech<<
1538               "isPrim="<<isPrim<<
1539               "isFromStrangess="<<isFromStrangess<<
1540               "isFromConversion="<<isFromConversion<<
1541               "isFromMaterial="<<isFromMaterial<<
1542               "particleTPC.="<<particleTPC<<
1543               "particleMotherTPC.="<<particleMotherTPC<<
1544               "mechTPC="<<mechTPC<<
1545               "isPrimTPC="<<isPrimTPC<<
1546               "isFromStrangessTPC="<<isFromStrangessTPC<<
1547               "isFromConversionTPC="<<isFromConversionTPC<<
1548               "isFromMaterialTPC="<<isFromMaterialTPC<<
1549               "particleITS.="<<particleITS<<
1550               "particleMotherITS.="<<particleMotherITS<<
1551               "mechITS="<<mechITS<<
1552               "isPrimITS="<<isPrimITS<<
1553               "isFromStrangessITS="<<isFromStrangessITS<<
1554               "isFromConversionITS="<<isFromConversionITS<<
1555               "isFromMaterialITS="<<isFromMaterialITS;
1556           }
1557           //finish writing the entry
1558           AliInfo("writing tree highPt");
1559           (*fTreeSRedirector)<<"highPt"<<"\n";
1560         }
1561         AliSysInfo::AddStamp("filteringTask",iTrack,numberOfTracks,numberOfFriendTracks,(friendTrack=fDummyFriendTrack)?0:1);
1562
1563         delete tpcInnerC;
1564         delete trackInnerC;
1565         delete trackInnerC2;
1566         delete outerITSc;
1567         delete trackInnerC3;
1568       }
1569   }
1570
1571 }
1572
1573
1574 //_____________________________________________________________________________
1575 void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1576 {
1577   //
1578   // Fill tree for efficiency studies MC only
1579   AliInfo("we start!");
1580
1581   if(!esdEvent) {
1582     AliDebug(AliLog::kError, "esdEvent not available");
1583     return;
1584   }
1585
1586   if(!mcEvent) {
1587     AliDebug(AliLog::kError, "mcEvent not available");
1588     return;
1589   }
1590
1591   // get selection cuts
1592   AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1593   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1594   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1595
1596   if(!evtCuts || !accCuts  || !esdTrackCuts) {
1597     AliDebug(AliLog::kError, "cuts not available");
1598     return;
1599   }
1600
1601   // trigger selection
1602   Bool_t isEventTriggered = kTRUE;
1603   AliPhysicsSelection *physicsSelection = NULL;
1604   AliTriggerAnalysis* triggerAnalysis = NULL;
1605
1606   //
1607   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1608   if (!inputHandler)
1609   {
1610     Printf("ERROR: Could not receive input handler");
1611     return;
1612   }
1613
1614   // get file name
1615   TTree *chain = (TChain*)GetInputData(0);
1616   if(!chain) {
1617     Printf("ERROR: Could not receive input chain");
1618     return;
1619   }
1620   TObjString fileName(chain->GetCurrentFile()->GetName());
1621
1622   // trigger
1623   if(evtCuts->IsTriggerRequired())
1624   {
1625     // always MB
1626     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1627
1628     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1629     if(!physicsSelection) return;
1630
1631     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1632       // set trigger (V0AND)
1633       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1634       if(!triggerAnalysis) return;
1635       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1636     }
1637   }
1638
1639   // centrality determination
1640   Float_t centralityF = -1;
1641   AliCentrality *esdCentrality = esdEvent->GetCentrality();
1642   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1643
1644   // use MC information
1645   AliHeader* header = 0;
1646   AliGenEventHeader* genHeader = 0;
1647   AliStack* stack = 0;
1648   Int_t mcStackSize=0;
1649   TArrayF vtxMC(3);
1650
1651   Int_t multMCTrueTracks = 0;
1652   //
1653   if(!mcEvent) {
1654     AliDebug(AliLog::kError, "mcEvent not available");
1655     return;
1656   }
1657   // get MC event header
1658   header = mcEvent->Header();
1659   if (!header) {
1660     AliDebug(AliLog::kError, "Header not available");
1661     return;
1662   }
1663   // MC particle stack
1664   stack = mcEvent->Stack();
1665   if (!stack) {
1666     AliDebug(AliLog::kError, "Stack not available");
1667     return;
1668   }
1669   mcStackSize=stack->GetNtrack();
1670
1671   // get MC vertex
1672   genHeader = header->GenEventHeader();
1673   if (!genHeader) {
1674     AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1675     return;
1676   }
1677   genHeader->PrimaryVertex(vtxMC);
1678
1679   // multipliticy of all MC primary tracks
1680   // in Zv, pt and eta ranges)
1681   multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1682
1683
1684   // get reconstructed vertex
1685   //const AliESDVertex* vtxESD = 0;
1686   AliESDVertex* vtxESD = 0;
1687   if(GetAnalysisMode() == kTPCAnalysisMode) {
1688     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1689   }
1690   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1691     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1692   }
1693   else {
1694     return;
1695   }
1696
1697   if(!vtxESD) return;
1698   //
1699   // Vertex info comparison and track multiplicity
1700   //
1701   AliESDVertex *vertexSPD =  (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
1702   AliESDVertex *vertexTPC =  (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
1703   Int_t contSPD = vertexSPD->GetNContributors();
1704   Int_t contTPC = vertexTPC->GetNContributors();
1705   TVectorD vertexPosTPC(3), vertexPosSPD(3);
1706   vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
1707   vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
1708   Int_t ntracksTPC=0;
1709   Int_t ntracksITS=0;
1710   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
1711     AliESDtrack *track = esdEvent->GetTrack(iTrack);
1712     if(!track) continue;
1713     if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
1714     if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
1715   }
1716
1717   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1718   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1719   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1720
1721   TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1722
1723   // check event cuts
1724   if(isEventOK && isEventTriggered)
1725   {
1726     //if(!stack) return;
1727
1728     //
1729     // MC info
1730     //
1731     TParticle *particle=NULL;
1732     TParticle *particleMother=NULL;
1733     Int_t mech=-1;
1734
1735     // reco event info
1736     Double_t vert[3] = {0};
1737     vert[0] = vtxESD->GetXv();
1738     vert[1] = vtxESD->GetYv();
1739     vert[2] = vtxESD->GetZv();
1740     Int_t mult = vtxESD->GetNContributors();
1741     Double_t bz = esdEvent->GetMagneticField();
1742     Double_t runNumber = esdEvent->GetRunNumber();
1743     Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1744     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1745     // loop over MC stack
1746     for (Int_t iMc = 0; iMc < mcStackSize; ++iMc)
1747     {
1748       particle = stack->Particle(iMc);
1749       if (!particle)
1750         continue;
1751
1752       // only charged particles
1753       if(!particle->GetPDG()) continue;
1754       Double_t charge = particle->GetPDG()->Charge()/3.;
1755       if (TMath::Abs(charge) < 0.001)
1756         continue;
1757
1758       // only primary particles
1759       Bool_t prim = stack->IsPhysicalPrimary(iMc);
1760       if(!prim) continue;
1761
1762       // downscale low-pT particles
1763       Double_t scalempt= TMath::Min(particle->Pt(),10.);
1764       Double_t downscaleF = gRandom->Rndm();
1765       downscaleF *= fLowPtTrackDownscaligF;
1766       if(TMath::Exp(2*scalempt)<downscaleF) continue;
1767
1768       // is particle in acceptance
1769       if(!accCuts->AcceptTrack(particle)) continue;
1770
1771       // check if particle reconstructed
1772       Bool_t isRec = kFALSE;
1773       Int_t trackIndex = -1;
1774       Int_t trackLoopIndex = -1;
1775       Int_t isESDtrackCut= 0;
1776       Int_t isAccCuts    = 0;
1777       Int_t nRec = 0;    // how many times reconstructed
1778       Int_t nFakes = 0;  // how many times reconstructed as a fake track
1779       AliESDtrack *recTrack = NULL;
1780
1781       for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1782       {
1783         AliESDtrack *track = esdEvent->GetTrack(iTrack);
1784         if(!track) continue;
1785         if(track->Charge()==0) continue;
1786         //
1787         Int_t label =  TMath::Abs(track->GetLabel());
1788         if (label >= mcStackSize) continue;
1789         if(label == iMc) {
1790           Bool_t isAcc=esdTrackCuts->AcceptTrack(track);
1791           if (isAcc) isESDtrackCut=1;
1792           if (accCuts->AcceptTrack(track)) isAccCuts=1;
1793           isRec = kTRUE;
1794           trackIndex = iTrack;
1795
1796           if (recTrack){
1797             if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track
1798             if (!isAcc) continue;
1799             trackLoopIndex = iTrack;
1800           }
1801           recTrack = esdEvent->GetTrack(trackIndex);
1802           nRec++;
1803           if(track->GetLabel()<0) nFakes++;
1804
1805           continue;
1806         }
1807       }
1808
1809       // Store information in the output tree
1810       if (trackLoopIndex>-1)  {
1811         recTrack = esdEvent->GetTrack(trackLoopIndex);
1812       } else if (trackIndex >-1) {
1813         recTrack = esdEvent->GetTrack(trackIndex);
1814       } else {
1815         recTrack = fDummyTrack;
1816       }
1817
1818       particleMother = GetMother(particle,stack);
1819       mech = particle->GetUniqueID();
1820
1821       //MC particle track length
1822       Double_t tpcTrackLength = 0.;
1823       AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
1824       if(mcParticle) {
1825         Int_t counter;
1826         tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
1827       }
1828
1829
1830       //
1831       if(fTreeSRedirector && fFillTree) {
1832         (*fTreeSRedirector)<<"MCEffTree"<<
1833           "fileName.="<<&fileName<<
1834           "triggerClass.="<<&triggerClass<<
1835           "runNumber="<<runNumber<<
1836           "evtTimeStamp="<<evtTimeStamp<<
1837           "evtNumberInFile="<<evtNumberInFile<<     //
1838           "Bz="<<bz<<                               // magnetic field
1839           "vtxESD.="<<vtxESD<<                      // vertex info
1840           //
1841           "mult="<<mult<<                           // primary vertex 9whatewe found) multiplicity
1842           "multMCTrueTracks="<<multMCTrueTracks<<   // mC track multiplicities
1843           //                                           important variables for the pile-up studies
1844           "contTPC="<< contTPC<<                    // number of contributors to the TPC primary vertex candidate
1845           "contSPD="<< contSPD<<                    // number of contributors to the SPD primary vertex candidate
1846           "vertexPosTPC.="<<&vertexPosTPC<<          // TPC vertex position
1847           "vertexPosSPD.="<<&vertexPosSPD<<          // SPD vertex position
1848           "ntracksTPC="<<ntracksTPC<<               // total number of the TPC tracks which were refitted
1849           "ntracksITS="<<ntracksITS<<               // total number of the ITS tracks which were refitted
1850           //
1851           //
1852           "isAcc0="<<isESDtrackCut<<                // track accepted by ESD track cuts
1853           "isAcc1="<<isAccCuts<<                    // track accepted by acceptance cuts flag
1854           "esdTrack.="<<recTrack<<                  // reconstructed track (only the longest from the loopers)
1855           "isRec="<<isRec<<                         // track was reconstructed
1856           "tpcTrackLength="<<tpcTrackLength<<       // track length in the TPC r projection
1857           "particle.="<<particle<<                  // particle properties
1858           "particleMother.="<<particleMother<<      // particle mother
1859           "mech="<<mech<<                           // production mechanizm
1860           "nRec="<<nRec<<                           // how many times reconstruted
1861           "nFakes="<<nFakes<<                       // how many times reconstructed as a fake track
1862           "\n";
1863       }
1864
1865       //if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
1866     }
1867   }
1868
1869 }
1870
1871 //_____________________________________________________________________________
1872 Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {
1873   //
1874   // check if particle is Z > 1
1875   //
1876   if (track->GetTPCNcls() < 60) return kFALSE;
1877   Double_t mom = track->GetInnerParam()->GetP();
1878   if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
1879   Float_t dca[2], bCov[3];
1880   track->GetImpactParameters(dca,bCov);
1881   //
1882
1883   Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
1884
1885   if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
1886
1887   return kFALSE;
1888 }
1889
1890 //_____________________________________________________________________________
1891 void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
1892 {
1893   //
1894   // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates
1895   //
1896   if(!esdEvent) {
1897     AliDebug(AliLog::kError, "esdEvent not available");
1898     return;
1899   }
1900
1901   // get selection cuts
1902   AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1903   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1904   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1905
1906   if(!evtCuts || !accCuts  || !esdTrackCuts) {
1907     AliDebug(AliLog::kError, "cuts not available");
1908     return;
1909   }
1910
1911   // trigger selection
1912   Bool_t isEventTriggered = kTRUE;
1913   AliPhysicsSelection *physicsSelection = NULL;
1914   AliTriggerAnalysis* triggerAnalysis = NULL;
1915
1916   //
1917   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1918   if (!inputHandler)
1919   {
1920     Printf("ERROR: Could not receive input handler");
1921     return;
1922   }
1923
1924   // get file name
1925   TTree *chain = (TChain*)GetInputData(0);
1926   if(!chain) {
1927     Printf("ERROR: Could not receive input chain");
1928     return;
1929   }
1930   TObjString fileName(chain->GetCurrentFile()->GetName());
1931
1932   // trigger
1933   if(evtCuts->IsTriggerRequired())
1934   {
1935     // always MB
1936     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1937
1938     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1939     if(!physicsSelection) return;
1940     //SetPhysicsTriggerSelection(physicsSelection);
1941
1942     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1943       // set trigger (V0AND)
1944       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1945       if(!triggerAnalysis) return;
1946       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1947     }
1948   }
1949
1950   // centrality determination
1951   Float_t centralityF = -1;
1952   AliCentrality *esdCentrality = esdEvent->GetCentrality();
1953   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1954
1955
1956   // get reconstructed vertex
1957   //const AliESDVertex* vtxESD = 0;
1958   AliESDVertex* vtxESD = 0;
1959   if(GetAnalysisMode() == kTPCAnalysisMode) {
1960     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1961   }
1962   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1963     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1964   }
1965   else {
1966     return;
1967   }
1968
1969   if(!vtxESD) return;
1970
1971   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1972   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1973   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1974
1975   // check event cuts
1976   if(isEventOK && isEventTriggered) {
1977     //
1978     // Dump the pt downscaled V0 into the tree
1979     //
1980     Int_t ntracks = esdEvent->GetNumberOfTracks();
1981     Int_t nV0s = esdEvent->GetNumberOfV0s();
1982     Int_t run = esdEvent->GetRunNumber();
1983     Int_t time= esdEvent->GetTimeStamp();
1984     Int_t evNr=esdEvent->GetEventNumberInFile();
1985     Double_t bz = esdEvent->GetMagneticField();
1986
1987     // Global event id calculation using orbitID, bunchCrossingID and periodID
1988     ULong64_t orbitID      = (ULong64_t)esdEvent->GetOrbitNumber();
1989     ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1990     ULong64_t periodID     = (ULong64_t)esdEvent->GetPeriodNumber();
1991     ULong64_t gid          = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1992
1993
1994     for (Int_t iv0=0; iv0<nV0s; iv0++){
1995
1996       AliESDv0 * v0 = esdEvent->GetV0(iv0);
1997       if (!v0) continue;
1998       AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
1999       AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
2000       if (!track0) continue;
2001       if (!track1) continue;
2002       AliESDfriendTrack* friendTrack0=NULL;
2003       AliESDfriendTrack* friendTrack1=NULL;
2004       if (esdFriend)
2005       {
2006         if (!esdFriend->TestSkipBit())
2007         {
2008           friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
2009           friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
2010         }
2011       }
2012       //Bool_t newFriendTrack0=kFALSE;
2013       //Bool_t newFriendTrack1=kFALSE;
2014       //if (!friendTrack0) {friendTrack0=new AliESDfriendTrack(); newFriendTrack0=kTRUE;}
2015       //if (!friendTrack1) {friendTrack1=new AliESDfriendTrack(); newFriendTrack1=kTRUE;}
2016       if (!friendTrack0) {friendTrack0=fDummyFriendTrack;}
2017       if (!friendTrack1) {friendTrack1=fDummyFriendTrack;}
2018       if (track0->GetSign()<0) {
2019         track1 = esdEvent->GetTrack(v0->GetIndex(0));
2020         track0 = esdEvent->GetTrack(v0->GetIndex(1));
2021       }
2022       //
2023       Bool_t isDownscaled = IsV0Downscaled(v0);
2024       if (isDownscaled) continue;
2025       AliKFParticle kfparticle; //
2026       Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
2027       if (type==0) continue;
2028       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2029
2030       if(!fFillTree) return;
2031       if(!fTreeSRedirector) return;
2032       (*fTreeSRedirector)<<"V0s"<<
2033         "gid="<<gid<<
2034         "isDownscaled="<<isDownscaled<<
2035         "triggerClass="<<&triggerClass<<      //  trigger
2036         "Bz="<<bz<<
2037         "fileName.="<<&fileName<<
2038         "runNumber="<<run<<
2039         "evtTimeStamp="<<time<<
2040         "evtNumberInFile="<<evNr<<
2041         "type="<<type<<
2042         "ntracks="<<ntracks<<
2043         "v0.="<<v0<<
2044         "kf.="<<&kfparticle<<
2045         "track0.="<<track0<<
2046         "track1.="<<track1<<
2047         "friendTrack0.="<<friendTrack0<<
2048         "friendTrack1.="<<friendTrack1<<
2049         "centralityF="<<centralityF<<
2050         "\n";
2051       //if (newFriendTrack0) {delete friendTrack0;}
2052       //if (newFriendTrack1) {delete friendTrack1;}
2053     }
2054   }
2055 }
2056
2057 //_____________________________________________________________________________
2058 void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
2059 {
2060   //
2061   // Select real events with large TPC dEdx signal
2062   //
2063   if(!esdEvent) {
2064     AliDebug(AliLog::kError, "esdEvent not available");
2065     return;
2066   }
2067
2068   // get selection cuts
2069   AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
2070   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
2071   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
2072
2073   if(!evtCuts || !accCuts  || !esdTrackCuts) {
2074     AliDebug(AliLog::kError, "cuts not available");
2075     return;
2076   }
2077
2078   // get file name
2079   TTree *chain = (TChain*)GetInputData(0);
2080   if(!chain) {
2081     Printf("ERROR: Could not receive input chain");
2082     return;
2083   }
2084   TObjString fileName(chain->GetCurrentFile()->GetName());
2085
2086   // trigger
2087   Bool_t isEventTriggered = kTRUE;
2088   AliPhysicsSelection *physicsSelection = NULL;
2089   AliTriggerAnalysis* triggerAnalysis = NULL;
2090
2091   //
2092   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
2093   if (!inputHandler)
2094   {
2095     Printf("ERROR: Could not receive input handler");
2096     return;
2097   }
2098
2099   if(evtCuts->IsTriggerRequired())
2100   {
2101     // always MB
2102     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
2103
2104     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
2105     if(!physicsSelection) return;
2106
2107     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
2108       // set trigger (V0AND)
2109       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
2110       if(!triggerAnalysis) return;
2111       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
2112     }
2113   }
2114
2115   // get reconstructed vertex
2116   AliESDVertex* vtxESD = 0;
2117   if(GetAnalysisMode() == kTPCAnalysisMode) {
2118     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
2119   }
2120   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
2121     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
2122   }
2123   else {
2124     return;
2125   }
2126   if(!vtxESD) return;
2127
2128   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
2129   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
2130   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
2131
2132
2133   // check event cuts
2134   if(isEventOK && isEventTriggered)
2135   {
2136     Double_t vert[3] = {0};
2137     vert[0] = vtxESD->GetXv();
2138     vert[1] = vtxESD->GetYv();
2139     vert[2] = vtxESD->GetZv();
2140     Int_t mult = vtxESD->GetNContributors();
2141     Double_t bz = esdEvent->GetMagneticField();
2142     Double_t runNumber = esdEvent->GetRunNumber();
2143     Double_t evtTimeStamp = esdEvent->GetTimeStamp();
2144     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
2145
2146     // Global event id calculation using orbitID, bunchCrossingID and periodID
2147     ULong64_t orbitID      = (ULong64_t)esdEvent->GetOrbitNumber();
2148     ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
2149     ULong64_t periodID     = (ULong64_t)esdEvent->GetPeriodNumber();
2150     ULong64_t gid          = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
2151
2152     // large dEdx
2153     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
2154     {
2155       AliESDtrack *track = esdEvent->GetTrack(iTrack);
2156       if(!track) continue;
2157       AliESDfriendTrack* friendTrack=NULL;
2158       if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
2159       //Bool_t newFriendTrack=kFALSE;
2160       //if (!friendTrack) {friendTrack=new AliESDfriendTrack(); newFriendTrack=kTRUE;}
2161       if (!friendTrack) {friendTrack=fDummyFriendTrack;}
2162       if(track->Charge()==0) continue;
2163       if(!esdTrackCuts->AcceptTrack(track)) continue;
2164       if(!accCuts->AcceptTrack(track)) continue;
2165
2166       if(!IsHighDeDxParticle(track)) continue;
2167       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2168
2169       if(!fFillTree) return;
2170       if(!fTreeSRedirector) return;
2171       (*fTreeSRedirector)<<"dEdx"<<
2172         "gid="<<gid<<
2173         "fileName.="<<&fileName<<
2174         "runNumber="<<runNumber<<
2175         "evtTimeStamp="<<evtTimeStamp<<
2176         "evtNumberInFile="<<evtNumberInFile<<
2177         "triggerClass="<<&triggerClass<<      //  trigger
2178         "Bz="<<bz<<
2179         "vtxESD.="<<vtxESD<<
2180         "mult="<<mult<<
2181         "esdTrack.="<<track<<
2182         "friendTrack.="<<friendTrack<<
2183         "\n";
2184       //if (newFriendTrack) delete friendTrack;
2185     }
2186   }
2187 }
2188
2189 //_____________________________________________________________________________
2190 Int_t   AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
2191 {
2192   //
2193   // Create KF particle in case the V0 fullfill selection criteria
2194   //
2195   // Selection criteria
2196   //  0. algorithm cut
2197   //  1. track cut
2198   //  3. chi2 cut
2199   //  4. rough mass cut
2200   //  5. Normalized pointing angle cut
2201   //
2202   const Double_t cutMass=0.2;
2203   const Double_t kSigmaDCACut=3;
2204   //
2205   // 0.) algo cut - accept only on the fly
2206   //
2207   if (v0->GetOnFlyStatus() ==kFALSE) return 0;
2208   //
2209   // 1.) track cut
2210   //
2211   AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
2212   AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
2213   /*
2214      TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
2215      TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
2216      TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
2217      */
2218   if (TMath::Abs(track0->GetTgl())>1) return 0;
2219   if (TMath::Abs(track1->GetTgl())>1) return 0;
2220   if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
2221   if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
2222   //if ((track0->GetITSclusters(0))<2) return 0;
2223   //if ((track1->GetITSclusters(0))<2) return 0;
2224   Float_t pos0[2]={0}, cov0[3]={0};
2225   Float_t pos1[2]={0}, cov1[3]={0};
2226   track0->GetImpactParameters(pos0,cov0);
2227   track0->GetImpactParameters(pos1,cov1);
2228   //
2229   if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
2230   if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
2231   //
2232   //
2233   // 3.) Chi2 cut
2234   //
2235   Double_t chi2KF = v0->GetKFInfo(2,2,2);
2236   if (chi2KF>25) return 0;
2237   //
2238   // 4.) Rough mass cut - 0.200 GeV
2239   //
2240   static Double_t masses[2]={-1};
2241   if (masses[0]<0){
2242     masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
2243     masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
2244   }
2245   Double_t mass00=  v0->GetEffMass(0,0);
2246   Double_t mass22=  v0->GetEffMass(2,2);
2247   Double_t mass42=  v0->GetEffMass(4,2);
2248   Double_t mass24=  v0->GetEffMass(2,4);
2249   Bool_t massOK=kFALSE;
2250   Int_t type=0;
2251   Int_t ptype=0;
2252   Double_t dmass=1;
2253   Int_t p1=0, p2=0;
2254   if (TMath::Abs(mass00-0)<cutMass) {
2255     massOK=kTRUE; type+=1;
2256     if (TMath::Abs(mass00-0)<dmass) {
2257       ptype=1;
2258       dmass=TMath::Abs(mass00-0);
2259       p1=0; p2=0;
2260     }
2261   }
2262   if (TMath::Abs(mass24-masses[1])<cutMass) {
2263     massOK=kTRUE; type+=2;
2264     if (TMath::Abs(mass24-masses[1])<dmass){
2265       dmass = TMath::Abs(mass24-masses[1]);
2266       ptype=2;
2267       p1=2; p2=4;
2268     }
2269   }
2270   if (TMath::Abs(mass42-masses[1])<cutMass) {
2271     massOK=kTRUE; type+=4;
2272     if (TMath::Abs(mass42-masses[1])<dmass){
2273       dmass = TMath::Abs(mass42-masses[1]);
2274       ptype=4;
2275       p1=4; p2=2;
2276     }
2277   }
2278   if (TMath::Abs(mass22-masses[0])<cutMass) {
2279     massOK=kTRUE; type+=8;
2280     if (TMath::Abs(mass22-masses[0])<dmass){
2281       dmass = TMath::Abs(mass22-masses[0]);
2282       ptype=8;
2283       p1=2; p2=2;
2284     }
2285   }
2286   if (type==0) return 0;
2287   //
2288   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
2289   const AliExternalTrackParam *paramP = v0->GetParamP();
2290   const AliExternalTrackParam *paramN = v0->GetParamN();
2291   if (paramP->GetSign()<0){
2292     paramP=v0->GetParamP();
2293     paramN=v0->GetParamN();
2294   }
2295   //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
2296   //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
2297   //
2298   AliKFParticle kfp1( *paramP, spdg[p1]  );
2299   AliKFParticle kfp2( *paramN, -1 *spdg[p2]  );
2300   AliKFParticle V0KF;
2301   (V0KF)+=kfp1;
2302   (V0KF)+=kfp2;
2303   kfparticle=V0KF;
2304   //
2305   // Pointing angle
2306   //
2307   Double_t  errPhi    = V0KF.GetErrPhi();
2308   Double_t  pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
2309   if (pointAngle/errPhi>10) return 0;
2310   //
2311   return ptype;
2312 }
2313
2314 //_____________________________________________________________________________
2315 Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)
2316 {
2317   //
2318   // Downscale randomly low pt V0
2319   //
2320   //return kFALSE;
2321   Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
2322   Double_t scalempt= TMath::Min(maxPt,10.);
2323   Double_t downscaleF = gRandom->Rndm();
2324   downscaleF *= fLowPtV0DownscaligF;
2325   //
2326   // Special treatment of the gamma conversion pt spectra is softer -
2327   Double_t mass00=  v0->GetEffMass(0,0);
2328   const Double_t cutMass=0.2;
2329   if (TMath::Abs(mass00-0)<cutMass){
2330     downscaleF/=10.;  // 10 times smaller downscaling for the gamma concersion candidate
2331   }
2332   //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
2333   if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;
2334   return kFALSE;
2335
2336   /*
2337      TH1F his1("his1","his1",100,0,10);
2338      TH1F his2("his2","his2",100,0,10);
2339      {for (Int_t i=0; i<10000; i++){
2340      Double_t rnd=gRandom->Exp(1);
2341      Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
2342      his1->Fill(rnd);
2343      if (!isDownscaled) his2->Fill(rnd);
2344      }}
2345
2346 */
2347
2348 }
2349
2350
2351
2352 //_____________________________________________________________________________
2353 Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
2354 {
2355   // Constrain TPC inner params constrained
2356   //
2357   if(!tpcInnerC) return kFALSE;
2358   if(!vtx) return kFALSE;
2359
2360   Double_t dz[2],cov[3];
2361   //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2362   //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2363   //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2364   if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
2365
2366
2367   Double_t covar[6]; vtx->GetCovMatrix(covar);
2368   Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
2369   Double_t c[3]={covar[2],0.,covar[5]};
2370   Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
2371   if (chi2C>kVeryBig) return kFALSE;
2372
2373   if(!tpcInnerC->Update(p,c)) return kFALSE;
2374
2375   return kTRUE;
2376 }
2377
2378 //_____________________________________________________________________________
2379 Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
2380 {
2381   // Constrain TPC inner params constrained
2382   //
2383   if(!trackInnerC) return kFALSE;
2384   if(!vtx) return kFALSE;
2385
2386   const Double_t kRadius  = 2.8;
2387   const Double_t kMaxStep = 1.0;
2388
2389   Double_t dz[2],cov[3];
2390
2391   //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2392   //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2393   //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2394
2395   if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
2396   if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
2397
2398   //
2399   Double_t covar[6]; vtx->GetCovMatrix(covar);
2400   Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
2401   Double_t c[3]={covar[2],0.,covar[5]};
2402   Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
2403   if (chi2C>kVeryBig) return kFALSE;
2404
2405   if(!trackInnerC->Update(p,c)) return kFALSE;
2406
2407   return kTRUE;
2408 }
2409
2410
2411 //_____________________________________________________________________________
2412 TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack)
2413 {
2414   if(!particle) return NULL;
2415   if(!stack) return NULL;
2416
2417   Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2418   TParticle* mother = NULL;
2419   Int_t mcStackSize=stack->GetNtrack();
2420   if (motherLabel>=mcStackSize) return NULL;
2421   mother = stack->Particle(motherLabel);
2422
2423   return mother;
2424 }
2425
2426 //_____________________________________________________________________________
2427 Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack)
2428 {
2429   Bool_t isFromConversion = kFALSE;
2430
2431   if(stack) {
2432     Int_t mcStackSize=stack->GetNtrack();
2433     if (label>=mcStackSize) return kFALSE;
2434     TParticle* particle = stack->Particle(label);
2435     if (!particle) return kFALSE;
2436
2437     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2438     {
2439       Int_t mech = particle->GetUniqueID(); // production mechanism
2440       Bool_t isPrim = stack->IsPhysicalPrimary(label);
2441
2442       Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2443       if (motherLabel>=mcStackSize) return kFALSE;
2444       TParticle* mother = stack->Particle(motherLabel);
2445       if(mother) {
2446         Int_t motherPdg = mother->GetPdgCode();
2447
2448         if(!isPrim && mech==5 && motherPdg==kGamma) {
2449           isFromConversion=kTRUE;
2450         }
2451       }
2452     }
2453   }
2454
2455   return isFromConversion;
2456 }
2457
2458 //_____________________________________________________________________________
2459 Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack)
2460 {
2461   Bool_t isFromMaterial = kFALSE;
2462
2463   if(stack) {
2464     Int_t mcStackSize=stack->GetNtrack();
2465     if (label>=mcStackSize) return kFALSE;
2466     TParticle* particle = stack->Particle(label);
2467     if (!particle) return kFALSE;
2468
2469     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2470     {
2471       Int_t mech = particle->GetUniqueID(); // production mechanism
2472       Bool_t isPrim = stack->IsPhysicalPrimary(label);
2473
2474       Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2475       if (motherLabel>=mcStackSize) return kFALSE;
2476       TParticle* mother = stack->Particle(motherLabel);
2477       if(mother) {
2478         if(!isPrim && mech==13) {
2479           isFromMaterial=kTRUE;
2480         }
2481       }
2482     }
2483   }
2484
2485   return isFromMaterial;
2486 }
2487
2488 //_____________________________________________________________________________
2489 Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack)
2490 {
2491   Bool_t isFromStrangeness = kFALSE;
2492
2493   if(stack) {
2494     Int_t mcStackSize=stack->GetNtrack();
2495     if (label>=mcStackSize) return kFALSE;
2496     TParticle* particle = stack->Particle(label);
2497     if (!particle) return kFALSE;
2498
2499     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2500     {
2501       Int_t mech = particle->GetUniqueID(); // production mechanism
2502       Bool_t isPrim = stack->IsPhysicalPrimary(label);
2503
2504       Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2505       if (motherLabel>=mcStackSize) return kFALSE;
2506       TParticle* mother = stack->Particle(motherLabel);
2507       if(mother) {
2508         Int_t motherPdg = mother->GetPdgCode();
2509
2510         // K+-, lambda, antilambda, K0s decays
2511         if(!isPrim && mech==4 &&
2512             (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
2513         {
2514           isFromStrangeness = kTRUE;
2515         }
2516       }
2517     }
2518   }
2519
2520   return isFromStrangeness;
2521 }
2522
2523
2524 //_____________________________________________________________________________
2525 void AliAnalysisTaskFilteredTree::FinishTaskOutput()
2526 {
2527   //
2528   // Called one at the end
2529   // locally on working node
2530   //
2531
2532   //// must be deleted to store trees
2533   //if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;
2534
2535   //// open temporary file and copy trees to the ouptut container
2536
2537   //TChain* chain = 0;
2538   ////
2539   //chain = new TChain("highPt");
2540   //if(chain) {
2541   //  chain->Add("jotwinow_Temp_Trees.root");
2542   //  fHighPtTree = chain->CopyTree("1");
2543   //  delete chain; chain=0;
2544   //}
2545   //if(fHighPtTree) fHighPtTree->Print();
2546
2547   ////
2548   //chain = new TChain("V0s");
2549   //if(chain) {
2550   //  chain->Add("jotwinow_Temp_Trees.root");
2551   //  fV0Tree = chain->CopyTree("1");
2552   //  delete chain; chain=0;
2553   //}
2554   //if(fV0Tree) fV0Tree->Print();
2555
2556   ////
2557   //chain = new TChain("dEdx");
2558   //if(chain) {
2559   //  chain->Add("jotwinow_Temp_Trees.root");
2560   //  fdEdxTree = chain->CopyTree("1");
2561   //  delete chain; chain=0;
2562   //}
2563   //if(fdEdxTree) fdEdxTree->Print();
2564
2565   ////
2566   //chain = new TChain("Laser");
2567   //if(chain) {
2568   //  chain->Add("jotwinow_Temp_Trees.root");
2569   //  fLaserTree = chain->CopyTree("1");
2570   //  delete chain; chain=0;
2571   //}
2572   //if(fLaserTree) fLaserTree->Print();
2573
2574   ////
2575   //chain = new TChain("MCEffTree");
2576   //if(chain) {
2577   //  chain->Add("jotwinow_Temp_Trees.root");
2578   //  fMCEffTree = chain->CopyTree("1");
2579   //  delete chain; chain=0;
2580   //}
2581   //if(fMCEffTree) fMCEffTree->Print();
2582
2583   ////
2584   //chain = new TChain("CosmicPairs");
2585   //if(chain) {
2586   //  chain->Add("jotwinow_Temp_Trees.root");
2587   //  fCosmicPairsTree = chain->CopyTree("1");
2588   //  delete chain; chain=0;
2589   //}
2590   //if(fCosmicPairsTree) fCosmicPairsTree->Print();
2591
2592   Bool_t deleteTrees=kTRUE;
2593   if ((AliAnalysisManager::GetAnalysisManager()))
2594   {
2595     if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
2596         AliAnalysisManager::kProofAnalysis)
2597       deleteTrees=kFALSE;
2598   }
2599   if (deleteTrees) delete fTreeSRedirector;
2600   fTreeSRedirector=NULL;
2601
2602   //OpenFile(1);
2603
2604   // Post output data.
2605   //PostData(1, fHighPtTree);
2606   //PostData(2, fV0Tree);
2607   //PostData(3, fdEdxTree);
2608   //PostData(4, fLaserTree);
2609   //PostData(5, fMCEffTree);
2610   //PostData(6, fCosmicPairsTree);
2611 }
2612
2613 //_____________________________________________________________________________
2614 void AliAnalysisTaskFilteredTree::Terminate(Option_t *)
2615 {
2616   // Called one at the end
2617   /*
2618      fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));
2619      if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;
2620      TChain* chain = new TChain("highPt");
2621      if(!chain) return;
2622      chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");
2623      TTree *tree = chain->CopyTree("1");
2624      if (chain) { delete chain; chain=0; }
2625      if(!tree) return;
2626      tree->Print();
2627      fOutputSummary = tree;
2628
2629      if (!fOutputSummary) {
2630      Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1));
2631      return;
2632      }
2633      */
2634
2635 }
2636
2637 //_____________________________________________________________________________
2638 Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
2639 {
2640   //
2641   // calculate mc event true track multiplicity
2642   //
2643   if(!mcEvent) return 0;
2644
2645   AliStack* stack = 0;
2646   Int_t mult = 0;
2647
2648   // MC particle stack
2649   stack = mcEvent->Stack();
2650   if (!stack) return 0;
2651
2652   //
2653   //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
2654   //
2655
2656   Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
2657   if(!isEventOK) return 0;
2658
2659   Int_t nPart  = stack->GetNtrack();
2660   for (Int_t iMc = 0; iMc < nPart; ++iMc)
2661   {
2662     TParticle* particle = stack->Particle(iMc);
2663     if (!particle)
2664       continue;
2665
2666     // only charged particles
2667     if(!particle->GetPDG()) continue;
2668     Double_t charge = particle->GetPDG()->Charge()/3.;
2669     if (TMath::Abs(charge) < 0.001)
2670       continue;
2671
2672     // physical primary
2673     Bool_t prim = stack->IsPhysicalPrimary(iMc);
2674     if(!prim) continue;
2675
2676     // checked accepted without pt cut
2677     //if(accCuts->AcceptTrack(particle))
2678     if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
2679     {
2680       mult++;
2681     }
2682   }
2683
2684   return mult;
2685 }
2686
2687 //_____________________________________________________________________________
2688 void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, const Double_t centralityF, const Double_t chi2TPCInnerC)
2689 {
2690   //
2691   // Fill pT relative resolution histograms for
2692   // TPC only, TPC only constrained to vertex and TPC+ITS tracking
2693   //
2694   if(!ptrack) return;
2695   if(!ptpcInnerC) return;
2696
2697   const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();
2698   if(!innerParam) return;
2699
2700   Float_t dxy, dz;
2701   ptrack->GetImpactParameters(dxy,dz);
2702
2703   // TPC+ITS primary tracks
2704   if( abs(ptrack->Eta())<0.8 &&
2705       ptrack->GetTPCClusterInfo(3,1)>120 &&
2706       ptrack->IsOn(0x40) &&
2707       ptrack->GetTPCclusters(0)>0.0 &&
2708       ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2709       //abs(innerParam->GetX())>0.0 &&
2710       //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2711       //abs(innerParam->GetTgl())<0.85 &&
2712       ptrack->IsOn(0x0004) &&
2713       ptrack->GetNcls(0)>0 &&
2714       ptrack->GetITSchi2()>0 &&
2715       sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 &&
2716       sqrt(chi2TPCInnerC)<6 &&
2717       (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) &&
2718       abs(dz)<2.0 &&
2719       abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )
2720   {
2721     fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2722     fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2723     fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2724   }
2725
2726   // TPC primary tracks
2727   // and TPC constrained primary tracks
2728
2729   AliExternalTrackParam *ptpcInner  = (AliExternalTrackParam *) ptrack->GetTPCInnerParam();
2730   if(!ptpcInner) return;
2731
2732
2733   Float_t dxyTPC, dzTPC;
2734   ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);
2735
2736   if( abs(ptrack->Eta())<0.8 &&
2737       ptrack->GetTPCClusterInfo(3,1)>120 &&
2738       ptrack->IsOn(0x40)&&
2739       ptrack->GetTPCclusters(0)>0.0 &&
2740       ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2741       //abs(innerParam->GetX())>0.0 &&
2742       //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2743       //abs(innerParam->GetTgl())<0.85 &&
2744       abs(dzTPC)<3.2 &&
2745       abs(dxyTPC)<2.4 )
2746   {
2747     // TPC only
2748     fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2749     fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2750     fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2751
2752     // TPC constrained to vertex
2753     fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2754     fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2755     fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2756   }
2757 }