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