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