PWGPP/AliAnalysisTaskFilteredTree.cxx

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