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