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