PWGPP/AliAnalysisTaskFilteredTree.cxx

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