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