PWGPP/AliAnalysisTaskFilteredTree.cxx

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