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