3e4cce3c3ca8f3b84c907e72e4a6cf193d30e3b5
[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=NULL;
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 AliESDtrack* dummytrack=NULL;
1402         if (!dummytrack) dummytrack=new AliESDtrack();
1403         static AliExternalTrackParam* dummyexternaltrackparam=NULL;
1404         if (!dummyexternaltrackparam) dummyexternaltrackparam=new AliExternalTrackParam();
1405         static AliTrackReference* dummytrackreference=NULL;
1406         if (!dummytrackreference) dummytrackreference=new AliTrackReference();
1407         static TParticle* dummyparticle=NULL;
1408         if (!dummyparticle) dummyparticle=new TParticle();
1409
1410         //assign the dummy objects if needed
1411         if (!track) {track=dummytrack;}
1412         if (!friendTrack) {friendTrack=fDummyFriendTrack;}
1413         if (mcEvent)
1414         {
1415           if (!refTPCIn) {refTPCIn=dummytrackreference;}
1416           if (!refITS) {refITS=dummytrackreference;}
1417           if (!particle) {particle=dummyparticle;}
1418           if (!particleMother) {particleMother=dummyparticle;}
1419           if (!particleTPC) {particleTPC=dummyparticle;}
1420           if (!particleMotherTPC) {particleMotherTPC=dummyparticle;}
1421           if (!particleITS) {particleITS=dummyparticle;}
1422           if (!particleMotherITS) {particleMotherITS=dummyparticle;}
1423         }
1424         //the following are guaranteed to exist:
1425         //if (!tpcInnerC) {tpcInnerC=dummyexternaltrackparam;}
1426         //if (!trackInnerC) {trackInnerC=dummyexternaltrackparam;}
1427         //if (!trackInnerC2) {trackInnerC2=dummyexternaltrackparam;}
1428         //if (!outerITSc) {outerITSc=dummyexternaltrackparam;}
1429         //if (!trackInnerC3) {trackInnerC3=dummyexternaltrackparam;}
1430         /////////////////////////
1431
1432         //Double_t vtxX=vtxESD->GetX();
1433         //Double_t vtxY=vtxESD->GetY();
1434         //Double_t vtxZ=vtxESD->GetZ();
1435
1436         //AliESDVertex* pvtxESD = (AliESDVertex*)vtxESD->Clone();
1437         //AliESDtrack* ptrack=(AliESDtrack*)track->Clone();
1438         //AliExternalTrackParam* ptpcInnerC = (AliExternalTrackParam*)tpcInnerC->Clone();
1439         //AliExternalTrackParam* ptrackInnerC = (AliExternalTrackParam*)trackInnerC->Clone();
1440         //AliExternalTrackParam* ptrackInnerC2 = (AliExternalTrackParam*)trackInnerC2->Clone();
1441         //AliExternalTrackParam* pouterITSc = (AliExternalTrackParam*)outerITSc->Clone();
1442         //AliExternalTrackParam* ptrackInnerC3 = (AliExternalTrackParam*)trackInnerC3->Clone();
1443         //AliESDVertex* pvtxESD = new AliESDVertex(*vtxESD);
1444         //AliESDtrack* ptrack= new AliESDtrack(*track);
1445         //AliExternalTrackParam* ptpcInnerC = new AliExternalTrackParam(*tpcInnerC);
1446         //AliExternalTrackParam* ptrackInnerC = new AliExternalTrackParam(*trackInnerC);
1447         //AliExternalTrackParam* ptrackInnerC2 = new AliExternalTrackParam(*trackInnerC2);
1448         //AliExternalTrackParam* pouterITSc =  new AliExternalTrackParam(*outerITSc);
1449         //AliExternalTrackParam* ptrackInnerC3 = new AliExternalTrackParam(*trackInnerC3);
1450
1451         Int_t ntracks = esdEvent->GetNumberOfTracks();
1452
1453         // Global event id calculation using orbitID, bunchCrossingID and periodID 
1454         ULong64_t orbitID      = (ULong64_t)esdEvent->GetOrbitNumber();
1455         ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1456         ULong64_t periodID     = (ULong64_t)esdEvent->GetPeriodNumber();
1457         ULong64_t gid          = ((periodID << 36) | (orbitID << 12) | bunchCrossID); 
1458
1459         // fill histograms
1460         FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0));
1461
1462         if(fTreeSRedirector && dumpToTree && fFillTree) 
1463         {
1464
1465           //if (friendTrack)
1466           //{
1467           //  const AliTrackPointArray* array = friendTrack->GetTrackPointArray();
1468           //  if (!array) {printf("we have a friend, but the ponits are empty\n"); continue;}
1469           //  if (friendTrack==fDummyFriendTrack) printf("using the dummy friend track\n");
1470           //  cout<<array->GetX()[0]<<" "<<array->GetX()[2]<<endl;
1471           //}
1472           //else
1473           //{
1474           //  printf("no friend track\n");
1475           //}
1476
1477
1478           (*fTreeSRedirector)<<"highPt"<<
1479             "gid="<<gid<<
1480             "fileName.="<<&fileName<<                // name of the chunk file (hopefully full)
1481             "runNumber="<<runNumber<<                // runNumber
1482             "evtTimeStamp="<<evtTimeStamp<<          // time stamp of event (in seconds)
1483             "evtNumberInFile="<<evtNumberInFile<<    // event number
1484             "triggerClass="<<&triggerClass<<         // trigger class as a string
1485             "Bz="<<bz<<                              // solenoid magnetic field in the z direction (in kGaus)
1486             "vtxESD.="<<vtxESD<<                    // vertexer ESD tracks (can be biased by TPC pileup tracks)
1487             //"vtxESDx="<<vtxX<<
1488             //"vtxESDy="<<vtxY<<
1489             //"vtxESDz="<<vtxZ<<
1490             "IRtot="<<ir1<<                         // interaction record (trigger) counters - coutner 1
1491             "IRint2="<<ir2<<                        // interaction record (trigger) coutners - counter 2
1492             "mult="<<mult<<                         // multiplicity of tracks pointing to the primary vertex
1493             "ntracks="<<ntracks<<                   // number of the esd tracks (to take into account the pileup in the TPC)
1494             //                                           important variables for the pile-up studies
1495             "contTPC="<< contTPC<<                    // number of contributors to the TPC primary vertex candidate
1496             "contSPD="<< contSPD<<                    // number of contributors to the SPD primary vertex candidate
1497             "vertexPosTPC.="<<&vertexPosTPC<<          // TPC vertex position
1498             "vertexPosSPD.="<<&vertexPosSPD<<          // SPD vertex position
1499             "ntracksTPC="<<ntracksTPC<<               // total number of the TPC tracks which were refitted
1500             "ntracksITS="<<ntracksITS<<               // total number of the ITS tracks which were refitted
1501             //
1502             "esdTrack.="<<track<<                  // esdTrack as used in the physical analysis
1503             "friendTrack.="<<friendTrack<<      // esdFriendTrack associated to the esdTrack
1504             "extTPCInnerC.="<<tpcInnerC<<          // ??? 
1505             "extInnerParamC.="<<trackInnerC<<      // ???
1506             "extInnerParam.="<<trackInnerC2<<      // ???
1507             "extOuterITS.="<<outerITSc<<           // ???
1508             "extInnerParamRef.="<<trackInnerC3<<   // ???
1509             "chi2TPCInnerC="<<chi2(0,0)<<           // chi2   of tracks ???
1510             "chi2InnerC="<<chi2trackC(0,0)<<        // chi2s  of tracks TPCinner to the combined
1511             "chi2OuterITS="<<chi2OuterITS(0,0)<<    // chi2s  of tracks TPC at inner wall to the ITSout
1512             "centralityF="<<centralityF;
1513           if (mcEvent)
1514           {
1515             AliTrackReference refDummy;
1516             if (!refITS) refITS = &refDummy;
1517             if (!refTRD) refTRD = &refDummy;
1518             if (!refTOF) refTOF = &refDummy;
1519             if (!refEMCAL) refEMCAL = &refDummy;
1520             if (!refPHOS) refPHOS = &refDummy;
1521             (*fTreeSRedirector)<<"highPt"<<     
1522               "multMCTrueTracks="<<multMCTrueTracks<<   // mC track multiplicities
1523               "nrefITS="<<nrefITS<<              // number of track references in the ITS
1524               "nrefTPC="<<nrefTPC<<              // number of track references in the TPC
1525               "nrefTRD="<<nrefTRD<<              // number of track references in the TRD
1526               "nrefTOF="<<nrefTOF<<              // number of track references in the TOF
1527               "nrefEMCAL="<<nrefEMCAL<<              // number of track references in the TOF
1528               "nrefPHOS="<<nrefPHOS<<              // number of track references in the TOF
1529               "refTPCIn.="<<refTPCIn<<
1530               "refTPCOut.="<<refTPCOut<<
1531               "refITS.="<<refITS<<          
1532               "refTRD.="<<refTRD<<          
1533               "refTOF.="<<refTOF<<          
1534               "refEMCAL.="<<refEMCAL<<      
1535               "refPHOS.="<<refPHOS<<        
1536               "particle.="<<particle<<
1537               "particleMother.="<<particleMother<<
1538               "mech="<<mech<<
1539               "isPrim="<<isPrim<<
1540               "isFromStrangess="<<isFromStrangess<<
1541               "isFromConversion="<<isFromConversion<<
1542               "isFromMaterial="<<isFromMaterial<<
1543               "particleTPC.="<<particleTPC<<
1544               "particleMotherTPC.="<<particleMotherTPC<<
1545               "mechTPC="<<mechTPC<<
1546               "isPrimTPC="<<isPrimTPC<<
1547               "isFromStrangessTPC="<<isFromStrangessTPC<<
1548               "isFromConversionTPC="<<isFromConversionTPC<<
1549               "isFromMaterialTPC="<<isFromMaterialTPC<<
1550               "particleITS.="<<particleITS<<
1551               "particleMotherITS.="<<particleMotherITS<<
1552               "mechITS="<<mechITS<<
1553               "isPrimITS="<<isPrimITS<<
1554               "isFromStrangessITS="<<isFromStrangessITS<<
1555               "isFromConversionITS="<<isFromConversionITS<<
1556               "isFromMaterialITS="<<isFromMaterialITS;
1557           }
1558           //finish writing the entry
1559           AliInfo("writing tree highPt");
1560           (*fTreeSRedirector)<<"highPt"<<"\n";
1561         }
1562         //AliSysInfo::AddStamp("a",iTrack,numberOfTracks,numberOfFriendTracks);
1563
1564         delete tpcInnerC;
1565         delete trackInnerC;
1566         delete trackInnerC2;
1567         delete outerITSc;
1568         delete trackInnerC3;
1569       }
1570   }
1571
1572 }
1573
1574
1575 //_____________________________________________________________________________
1576 void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1577 {
1578   //
1579   // Fill tree for efficiency studies MC only
1580   AliInfo("we start!");
1581
1582   if(!esdEvent) {
1583     AliDebug(AliLog::kError, "esdEvent not available");
1584     return;
1585   }
1586
1587   if(!mcEvent) {
1588     AliDebug(AliLog::kError, "mcEvent not available");
1589     return;
1590   }
1591
1592   // get selection cuts
1593   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); 
1594   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); 
1595   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); 
1596
1597   if(!evtCuts || !accCuts  || !esdTrackCuts) {
1598     AliDebug(AliLog::kError, "cuts not available");
1599     return;
1600   }
1601
1602   // trigger selection
1603   Bool_t isEventTriggered = kTRUE;
1604   AliPhysicsSelection *physicsSelection = NULL;
1605   AliTriggerAnalysis* triggerAnalysis = NULL;
1606
1607   // 
1608   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1609   if (!inputHandler)
1610   {
1611     Printf("ERROR: Could not receive input handler");
1612     return;
1613   }
1614
1615   // get file name
1616   TTree *chain = (TChain*)GetInputData(0);
1617   if(!chain) { 
1618     Printf("ERROR: Could not receive input chain");
1619     return;
1620   }
1621   TObjString fileName(chain->GetCurrentFile()->GetName());
1622
1623   // trigger
1624   if(evtCuts->IsTriggerRequired())  
1625   {
1626     // always MB
1627     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1628
1629     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1630     if(!physicsSelection) return;
1631
1632     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1633       // set trigger (V0AND)
1634       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1635       if(!triggerAnalysis) return;
1636       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1637     }
1638   }
1639
1640   // centrality determination
1641   Float_t centralityF = -1;
1642   AliCentrality *esdCentrality = esdEvent->GetCentrality();
1643   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1644
1645   // use MC information
1646   AliHeader* header = 0;
1647   AliGenEventHeader* genHeader = 0;
1648   AliStack* stack = 0;
1649   Int_t mcStackSize=0;
1650   TArrayF vtxMC(3);
1651
1652   Int_t multMCTrueTracks = 0;
1653   //
1654   if(!mcEvent) {
1655     AliDebug(AliLog::kError, "mcEvent not available");
1656     return;
1657   }
1658   // get MC event header
1659   header = mcEvent->Header();
1660   if (!header) {
1661     AliDebug(AliLog::kError, "Header not available");
1662     return;
1663   }
1664   // MC particle stack
1665   stack = mcEvent->Stack();
1666   if (!stack) {
1667     AliDebug(AliLog::kError, "Stack not available");
1668     return;
1669   }
1670   mcStackSize=stack->GetNtrack();
1671
1672   // get MC vertex
1673   genHeader = header->GenEventHeader();
1674   if (!genHeader) {
1675     AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1676     return;
1677   }
1678   genHeader->PrimaryVertex(vtxMC);
1679
1680   // multipliticy of all MC primary tracks
1681   // in Zv, pt and eta ranges)
1682   multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1683
1684
1685   // get reconstructed vertex  
1686   //const AliESDVertex* vtxESD = 0; 
1687   AliESDVertex* vtxESD = 0; 
1688   if(GetAnalysisMode() == kTPCAnalysisMode) {
1689     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1690   }
1691   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1692     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1693   }
1694   else {
1695     return;
1696   }
1697
1698   if(!vtxESD) return;   
1699   //
1700   // Vertex info comparison and track multiplicity
1701   //
1702   AliESDVertex *vertexSPD =  (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
1703   AliESDVertex *vertexTPC =  (AliESDVertex *)esdEvent->GetPrimaryVertexTPC(); 
1704   Int_t contSPD = vertexSPD->GetNContributors();
1705   Int_t contTPC = vertexTPC->GetNContributors();        
1706   TVectorD vertexPosTPC(3), vertexPosSPD(3);
1707   vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
1708   vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
1709   Int_t ntracksTPC=0;
1710   Int_t ntracksITS=0;
1711   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
1712     AliESDtrack *track = esdEvent->GetTrack(iTrack);    
1713     if(!track) continue;
1714     if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
1715     if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
1716   }
1717
1718   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); 
1719   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1720   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1721
1722   TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1723
1724   // check event cuts
1725   if(isEventOK && isEventTriggered)
1726   {
1727     //if(!stack) return;
1728
1729     //
1730     // MC info
1731     //
1732     TParticle *particle=NULL;
1733     TParticle *particleMother=NULL;
1734     Int_t mech=-1;
1735
1736     // reco event info
1737     Double_t vert[3] = {0}; 
1738     vert[0] = vtxESD->GetXv();
1739     vert[1] = vtxESD->GetYv();
1740     vert[2] = vtxESD->GetZv();
1741     Int_t mult = vtxESD->GetNContributors();
1742     Double_t bz = esdEvent->GetMagneticField();
1743     Double_t runNumber = esdEvent->GetRunNumber();
1744     Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1745     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1746     // loop over MC stack
1747     for (Int_t iMc = 0; iMc < mcStackSize; ++iMc) 
1748     {
1749       particle = stack->Particle(iMc);
1750       if (!particle)
1751         continue;
1752
1753       // only charged particles
1754       if(!particle->GetPDG()) continue;
1755       Double_t charge = particle->GetPDG()->Charge()/3.;
1756       if (TMath::Abs(charge) < 0.001)
1757         continue;
1758
1759       // only primary particles
1760       Bool_t prim = stack->IsPhysicalPrimary(iMc);
1761       if(!prim) continue;
1762
1763       // downscale low-pT particles
1764       Double_t scalempt= TMath::Min(particle->Pt(),10.);
1765       Double_t downscaleF = gRandom->Rndm();
1766       downscaleF *= fLowPtTrackDownscaligF;
1767       if(TMath::Exp(2*scalempt)<downscaleF) continue;
1768
1769       // is particle in acceptance
1770       if(!accCuts->AcceptTrack(particle)) continue;
1771
1772       // check if particle reconstructed
1773       Bool_t isRec = kFALSE;
1774       Int_t trackIndex = -1;  
1775       Int_t trackLoopIndex = -1;  
1776       Int_t isESDtrackCut= 0;
1777       Int_t isAccCuts    = 0;
1778       Int_t nRec = 0;    // how many times reconstructed 
1779       Int_t nFakes = 0;  // how many times reconstructed as a fake track
1780       AliESDtrack *recTrack = NULL; 
1781
1782       for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1783       {
1784         AliESDtrack *track = esdEvent->GetTrack(iTrack);
1785         if(!track) continue;
1786         if(track->Charge()==0) continue;
1787         //
1788         Int_t label =  TMath::Abs(track->GetLabel());
1789         if (label >= mcStackSize) continue;
1790         if(label == iMc) {        
1791           Bool_t isAcc=esdTrackCuts->AcceptTrack(track);
1792           if (isAcc) isESDtrackCut=1;
1793           if (accCuts->AcceptTrack(track)) isAccCuts=1;
1794           isRec = kTRUE;
1795           trackIndex = iTrack;
1796
1797           if (recTrack){
1798             if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track
1799             if (!isAcc) continue;
1800             trackLoopIndex = iTrack;
1801           }
1802           recTrack = esdEvent->GetTrack(trackIndex); 
1803           nRec++;
1804           if(track->GetLabel()<0) nFakes++;
1805
1806           continue;
1807         }        
1808       }
1809
1810       // Store information in the output tree
1811       if (trackLoopIndex>-1)  { 
1812         recTrack = esdEvent->GetTrack(trackLoopIndex); 
1813       } else if (trackIndex >-1) {
1814         recTrack = esdEvent->GetTrack(trackIndex); 
1815       } else {
1816         recTrack = fDummyTrack; 
1817       } 
1818
1819       particleMother = GetMother(particle,stack);
1820       mech = particle->GetUniqueID();
1821
1822       //MC particle track length
1823       Double_t tpcTrackLength = 0.;
1824       AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
1825       if(mcParticle) {
1826         Int_t counter;
1827         tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
1828       } 
1829
1830
1831       //
1832       if(fTreeSRedirector && fFillTree) {
1833         (*fTreeSRedirector)<<"MCEffTree"<<
1834           "fileName.="<<&fileName<<
1835           "triggerClass.="<<&triggerClass<<
1836           "runNumber="<<runNumber<<
1837           "evtTimeStamp="<<evtTimeStamp<<
1838           "evtNumberInFile="<<evtNumberInFile<<     // 
1839           "Bz="<<bz<<                               // magnetic field
1840           "vtxESD.="<<vtxESD<<                      // vertex info
1841           //
1842           "mult="<<mult<<                           // primary vertex 9whatewe found) multiplicity
1843           "multMCTrueTracks="<<multMCTrueTracks<<   // mC track multiplicities
1844           //                                           important variables for the pile-up studies
1845           "contTPC="<< contTPC<<                    // number of contributors to the TPC primary vertex candidate
1846           "contSPD="<< contSPD<<                    // number of contributors to the SPD primary vertex candidate
1847           "vertexPosTPC.="<<&vertexPosTPC<<          // TPC vertex position
1848           "vertexPosSPD.="<<&vertexPosSPD<<          // SPD vertex position
1849           "ntracksTPC="<<ntracksTPC<<               // total number of the TPC tracks which were refitted
1850           "ntracksITS="<<ntracksITS<<               // total number of the ITS tracks which were refitted
1851           //
1852           //
1853           "isAcc0="<<isESDtrackCut<<                // track accepted by ESD track cuts
1854           "isAcc1="<<isAccCuts<<                    // track accepted by acceptance cuts flag
1855           "esdTrack.="<<recTrack<<                  // reconstructed track (only the longest from the loopers)
1856           "isRec="<<isRec<<                         // track was reconstructed
1857           "tpcTrackLength="<<tpcTrackLength<<       // track length in the TPC r projection
1858           "particle.="<<particle<<                  // particle properties
1859           "particleMother.="<<particleMother<<      // particle mother
1860           "mech="<<mech<<                           // production mechanizm
1861           "nRec="<<nRec<<                           // how many times reconstruted
1862           "nFakes="<<nFakes<<                       // how many times reconstructed as a fake track
1863           "\n";
1864       }
1865
1866       //if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
1867     }
1868   }
1869
1870 }
1871
1872 //_____________________________________________________________________________
1873 Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {
1874   //
1875   // check if particle is Z > 1 
1876   //
1877   if (track->GetTPCNcls() < 60) return kFALSE;
1878   Double_t mom = track->GetInnerParam()->GetP();
1879   if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
1880   Float_t dca[2], bCov[3];
1881   track->GetImpactParameters(dca,bCov);
1882   //
1883
1884   Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
1885
1886   if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
1887
1888   return kFALSE;
1889 }
1890
1891 //_____________________________________________________________________________
1892 void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
1893 {
1894   //
1895   // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates
1896   //
1897   if(!esdEvent) {
1898     AliDebug(AliLog::kError, "esdEvent not available");
1899     return;
1900   }
1901
1902   // get selection cuts
1903   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); 
1904   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); 
1905   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); 
1906
1907   if(!evtCuts || !accCuts  || !esdTrackCuts) {
1908     AliDebug(AliLog::kError, "cuts not available");
1909     return;
1910   }
1911
1912   // trigger selection
1913   Bool_t isEventTriggered = kTRUE;
1914   AliPhysicsSelection *physicsSelection = NULL;
1915   AliTriggerAnalysis* triggerAnalysis = NULL;
1916
1917   // 
1918   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1919   if (!inputHandler)
1920   {
1921     Printf("ERROR: Could not receive input handler");
1922     return;
1923   }
1924
1925   // get file name
1926   TTree *chain = (TChain*)GetInputData(0);
1927   if(!chain) { 
1928     Printf("ERROR: Could not receive input chain");
1929     return;
1930   }
1931   TObjString fileName(chain->GetCurrentFile()->GetName());
1932
1933   // trigger
1934   if(evtCuts->IsTriggerRequired())  
1935   {
1936     // always MB
1937     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1938
1939     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1940     if(!physicsSelection) return;
1941     //SetPhysicsTriggerSelection(physicsSelection);
1942
1943     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1944       // set trigger (V0AND)
1945       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1946       if(!triggerAnalysis) return;
1947       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1948     }
1949   }
1950
1951   // centrality determination
1952   Float_t centralityF = -1;
1953   AliCentrality *esdCentrality = esdEvent->GetCentrality();
1954   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1955
1956
1957   // get reconstructed vertex  
1958   //const AliESDVertex* vtxESD = 0; 
1959   AliESDVertex* vtxESD = 0; 
1960   if(GetAnalysisMode() == kTPCAnalysisMode) {
1961     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1962   }
1963   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1964     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1965   }
1966   else {
1967     return;
1968   }
1969
1970   if(!vtxESD) return;
1971
1972   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); 
1973   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1974   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1975
1976   // check event cuts
1977   if(isEventOK && isEventTriggered) {
1978     //
1979     // Dump the pt downscaled V0 into the tree
1980     // 
1981     Int_t ntracks = esdEvent->GetNumberOfTracks();
1982     Int_t nV0s = esdEvent->GetNumberOfV0s();
1983     Int_t run = esdEvent->GetRunNumber();
1984     Int_t time= esdEvent->GetTimeStamp();
1985     Int_t evNr=esdEvent->GetEventNumberInFile();
1986     Double_t bz = esdEvent->GetMagneticField();
1987
1988     // Global event id calculation using orbitID, bunchCrossingID and periodID
1989     ULong64_t orbitID      = (ULong64_t)esdEvent->GetOrbitNumber();
1990     ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1991     ULong64_t periodID     = (ULong64_t)esdEvent->GetPeriodNumber();
1992     ULong64_t gid          = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1993
1994
1995     for (Int_t iv0=0; iv0<nV0s; iv0++){
1996  
1997       AliESDv0 * v0 = esdEvent->GetV0(iv0);
1998       if (!v0) continue;
1999       AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
2000       AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
2001       if (!track0) continue;
2002       if (!track1) continue;
2003       AliESDfriendTrack* friendTrack0=NULL;
2004       AliESDfriendTrack* friendTrack1=NULL;
2005       if (esdFriend) 
2006       {
2007         if (!esdFriend->TestSkipBit())
2008         {
2009           friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
2010           friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
2011         }
2012       }
2013       //Bool_t newFriendTrack0=kFALSE;
2014       //Bool_t newFriendTrack1=kFALSE;
2015       //if (!friendTrack0) {friendTrack0=new AliESDfriendTrack(); newFriendTrack0=kTRUE;}
2016       //if (!friendTrack1) {friendTrack1=new AliESDfriendTrack(); newFriendTrack1=kTRUE;}
2017       if (!friendTrack0) {friendTrack0=fDummyFriendTrack;}
2018       if (!friendTrack1) {friendTrack1=fDummyFriendTrack;}
2019       if (track0->GetSign()<0) {
2020         track1 = esdEvent->GetTrack(v0->GetIndex(0));
2021         track0 = esdEvent->GetTrack(v0->GetIndex(1));
2022       }
2023       //
2024       Bool_t isDownscaled = IsV0Downscaled(v0);
2025       if (isDownscaled) continue;
2026       AliKFParticle kfparticle; //
2027       Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
2028       if (type==0) continue;   
2029       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2030
2031       if(!fFillTree) return;
2032       if(!fTreeSRedirector) return;
2033       (*fTreeSRedirector)<<"V0s"<<
2034         "gid="<<gid<<
2035         "isDownscaled="<<isDownscaled<<
2036         "triggerClass="<<&triggerClass<<      //  trigger
2037         "Bz="<<bz<<
2038         "fileName.="<<&fileName<<
2039         "runNumber="<<run<<
2040         "evtTimeStamp="<<time<<
2041         "evtNumberInFile="<<evNr<<
2042         "type="<<type<<
2043         "ntracks="<<ntracks<<
2044         "v0.="<<v0<<
2045         "kf.="<<&kfparticle<<
2046         "track0.="<<track0<<
2047         "track1.="<<track1<<
2048         "friendTrack0.="<<friendTrack0<<
2049         "friendTrack1.="<<friendTrack1<<
2050         "centralityF="<<centralityF<<
2051         "\n";
2052       //if (newFriendTrack0) {delete friendTrack0;}
2053       //if (newFriendTrack1) {delete friendTrack1;}
2054     }
2055   }
2056 }
2057
2058 //_____________________________________________________________________________
2059 void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
2060 {
2061   //
2062   // Select real events with large TPC dEdx signal
2063   //
2064   if(!esdEvent) {
2065     AliDebug(AliLog::kError, "esdEvent not available");
2066     return;
2067   }
2068
2069   // get selection cuts
2070   AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); 
2071   AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); 
2072   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); 
2073
2074   if(!evtCuts || !accCuts  || !esdTrackCuts) {
2075     AliDebug(AliLog::kError, "cuts not available");
2076     return;
2077   }
2078
2079   // get file name
2080   TTree *chain = (TChain*)GetInputData(0);
2081   if(!chain) { 
2082     Printf("ERROR: Could not receive input chain");
2083     return;
2084   }
2085   TObjString fileName(chain->GetCurrentFile()->GetName());
2086
2087   // trigger
2088   Bool_t isEventTriggered = kTRUE;
2089   AliPhysicsSelection *physicsSelection = NULL;
2090   AliTriggerAnalysis* triggerAnalysis = NULL;
2091
2092   // 
2093   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
2094   if (!inputHandler)
2095   {
2096     Printf("ERROR: Could not receive input handler");
2097     return;
2098   }
2099
2100   if(evtCuts->IsTriggerRequired())  
2101   {
2102     // always MB
2103     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
2104
2105     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
2106     if(!physicsSelection) return;
2107
2108     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
2109       // set trigger (V0AND)
2110       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
2111       if(!triggerAnalysis) return;
2112       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
2113     }
2114   }
2115
2116   // get reconstructed vertex  
2117   AliESDVertex* vtxESD = 0; 
2118   if(GetAnalysisMode() == kTPCAnalysisMode) {
2119     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
2120   }
2121   else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
2122     vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
2123   }
2124   else {
2125     return;
2126   }
2127   if(!vtxESD) return;
2128
2129   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); 
2130   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
2131   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
2132
2133
2134   // check event cuts
2135   if(isEventOK && isEventTriggered)
2136   {
2137     Double_t vert[3] = {0}; 
2138     vert[0] = vtxESD->GetXv();
2139     vert[1] = vtxESD->GetYv();
2140     vert[2] = vtxESD->GetZv();
2141     Int_t mult = vtxESD->GetNContributors();
2142     Double_t bz = esdEvent->GetMagneticField();
2143     Double_t runNumber = esdEvent->GetRunNumber();
2144     Double_t evtTimeStamp = esdEvent->GetTimeStamp();
2145     Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
2146    
2147     // Global event id calculation using orbitID, bunchCrossingID and periodID 
2148     ULong64_t orbitID      = (ULong64_t)esdEvent->GetOrbitNumber();
2149     ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
2150     ULong64_t periodID     = (ULong64_t)esdEvent->GetPeriodNumber();
2151     ULong64_t gid          = ((periodID << 36) | (orbitID << 12) | bunchCrossID); 
2152     
2153     // large dEdx
2154     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
2155     {
2156       AliESDtrack *track = esdEvent->GetTrack(iTrack);
2157       if(!track) continue;
2158       AliESDfriendTrack* friendTrack=NULL;
2159       if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
2160       //Bool_t newFriendTrack=kFALSE;
2161       //if (!friendTrack) {friendTrack=new AliESDfriendTrack(); newFriendTrack=kTRUE;}
2162       if (!friendTrack) {friendTrack=fDummyFriendTrack;}
2163       if(track->Charge()==0) continue;
2164       if(!esdTrackCuts->AcceptTrack(track)) continue;
2165       if(!accCuts->AcceptTrack(track)) continue;
2166
2167       if(!IsHighDeDxParticle(track)) continue;
2168       TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2169
2170       if(!fFillTree) return;
2171       if(!fTreeSRedirector) return;
2172       (*fTreeSRedirector)<<"dEdx"<<
2173         "gid="<<gid<<
2174         "fileName.="<<&fileName<<
2175         "runNumber="<<runNumber<<
2176         "evtTimeStamp="<<evtTimeStamp<<
2177         "evtNumberInFile="<<evtNumberInFile<<
2178         "triggerClass="<<&triggerClass<<      //  trigger
2179         "Bz="<<bz<<
2180         "vtxESD.="<<vtxESD<<
2181         "mult="<<mult<<
2182         "esdTrack.="<<track<<
2183         "friendTrack.="<<friendTrack<<
2184         "\n";
2185       //if (newFriendTrack) delete friendTrack;
2186     }
2187   }
2188 }
2189
2190 //_____________________________________________________________________________
2191 Int_t   AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
2192 {
2193   //
2194   // Create KF particle in case the V0 fullfill selection criteria
2195   //
2196   // Selection criteria
2197   //  0. algorithm cut
2198   //  1. track cut
2199   //  3. chi2 cut
2200   //  4. rough mass cut
2201   //  5. Normalized pointing angle cut
2202   //
2203   const Double_t cutMass=0.2;
2204   const Double_t kSigmaDCACut=3;
2205   //
2206   // 0.) algo cut - accept only on the fly
2207   //
2208   if (v0->GetOnFlyStatus() ==kFALSE) return 0;     
2209   //
2210   // 1.) track cut
2211   // 
2212   AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
2213   AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
2214   /*
2215      TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
2216      TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
2217      TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
2218      */  
2219   if (TMath::Abs(track0->GetTgl())>1) return 0;
2220   if (TMath::Abs(track1->GetTgl())>1) return 0;
2221   if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
2222   if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
2223   //if ((track0->GetITSclusters(0))<2) return 0;
2224   //if ((track1->GetITSclusters(0))<2) return 0; 
2225   Float_t pos0[2]={0}, cov0[3]={0};
2226   Float_t pos1[2]={0}, cov1[3]={0};
2227   track0->GetImpactParameters(pos0,cov0);
2228   track0->GetImpactParameters(pos1,cov1);
2229   //
2230   if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
2231   if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
2232   // 
2233   //
2234   // 3.) Chi2 cut
2235   //
2236   Double_t chi2KF = v0->GetKFInfo(2,2,2);
2237   if (chi2KF>25) return 0;
2238   //
2239   // 4.) Rough mass cut - 0.200 GeV
2240   //
2241   static Double_t masses[2]={-1};
2242   if (masses[0]<0){
2243     masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
2244     masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
2245   }
2246   Double_t mass00=  v0->GetEffMass(0,0);
2247   Double_t mass22=  v0->GetEffMass(2,2);
2248   Double_t mass42=  v0->GetEffMass(4,2);
2249   Double_t mass24=  v0->GetEffMass(2,4);
2250   Bool_t massOK=kFALSE;
2251   Int_t type=0;
2252   Int_t ptype=0;
2253   Double_t dmass=1;
2254   Int_t p1=0, p2=0;
2255   if (TMath::Abs(mass00-0)<cutMass) {
2256     massOK=kTRUE; type+=1; 
2257     if (TMath::Abs(mass00-0)<dmass) {
2258       ptype=1;
2259       dmass=TMath::Abs(mass00-0);      
2260       p1=0; p2=0;
2261     } 
2262   }
2263   if (TMath::Abs(mass24-masses[1])<cutMass) {
2264     massOK=kTRUE; type+=2; 
2265     if (TMath::Abs(mass24-masses[1])<dmass){
2266       dmass = TMath::Abs(mass24-masses[1]);
2267       ptype=2;
2268       p1=2; p2=4;
2269     }
2270   }
2271   if (TMath::Abs(mass42-masses[1])<cutMass) {
2272     massOK=kTRUE; type+=4;
2273     if (TMath::Abs(mass42-masses[1])<dmass){
2274       dmass = TMath::Abs(mass42-masses[1]);
2275       ptype=4;
2276       p1=4; p2=2;
2277     }
2278   }
2279   if (TMath::Abs(mass22-masses[0])<cutMass) {
2280     massOK=kTRUE; type+=8;
2281     if (TMath::Abs(mass22-masses[0])<dmass){
2282       dmass = TMath::Abs(mass22-masses[0]);
2283       ptype=8;
2284       p1=2; p2=2;
2285     }
2286   }
2287   if (type==0) return 0;
2288   //
2289   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
2290   const AliExternalTrackParam *paramP = v0->GetParamP();
2291   const AliExternalTrackParam *paramN = v0->GetParamN();
2292   if (paramP->GetSign()<0){
2293     paramP=v0->GetParamP();
2294     paramN=v0->GetParamN();
2295   }
2296   //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
2297   //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
2298   //
2299   AliKFParticle kfp1( *paramP, spdg[p1]  );
2300   AliKFParticle kfp2( *paramN, -1 *spdg[p2]  );
2301   AliKFParticle V0KF;
2302   (V0KF)+=kfp1;
2303   (V0KF)+=kfp2;
2304   kfparticle=V0KF;
2305   //
2306   // Pointing angle
2307   //
2308   Double_t  errPhi    = V0KF.GetErrPhi();
2309   Double_t  pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
2310   if (pointAngle/errPhi>10) return 0;  
2311   //
2312   return ptype;  
2313 }
2314
2315 //_____________________________________________________________________________
2316 Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)
2317 {
2318   //
2319   // Downscale randomly low pt V0
2320   //
2321   //return kFALSE;
2322   Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
2323   Double_t scalempt= TMath::Min(maxPt,10.);
2324   Double_t downscaleF = gRandom->Rndm();
2325   downscaleF *= fLowPtV0DownscaligF;
2326   //
2327   // Special treatment of the gamma conversion pt spectra is softer - 
2328   Double_t mass00=  v0->GetEffMass(0,0);
2329   const Double_t cutMass=0.2;
2330   if (TMath::Abs(mass00-0)<cutMass){
2331     downscaleF/=10.;  // 10 times smaller downscaling for the gamma concersion candidate
2332   }
2333   //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
2334   if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;
2335   return kFALSE;
2336
2337   /*
2338      TH1F his1("his1","his1",100,0,10);
2339      TH1F his2("his2","his2",100,0,10);
2340      {for (Int_t i=0; i<10000; i++){
2341      Double_t rnd=gRandom->Exp(1);
2342      Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
2343      his1->Fill(rnd); 
2344      if (!isDownscaled) his2->Fill(rnd); 
2345      }}
2346
2347 */
2348
2349 }
2350
2351
2352
2353 //_____________________________________________________________________________
2354 Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
2355 {
2356   // Constrain TPC inner params constrained
2357   //
2358   if(!tpcInnerC) return kFALSE; 
2359   if(!vtx) return kFALSE;
2360
2361   Double_t dz[2],cov[3];
2362   //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2363   //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; 
2364   //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; 
2365   if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; 
2366
2367
2368   Double_t covar[6]; vtx->GetCovMatrix(covar);
2369   Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
2370   Double_t c[3]={covar[2],0.,covar[5]};
2371   Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
2372   if (chi2C>kVeryBig) return kFALSE; 
2373
2374   if(!tpcInnerC->Update(p,c)) return kFALSE;
2375
2376   return kTRUE;
2377 }
2378
2379 //_____________________________________________________________________________
2380 Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
2381 {
2382   // Constrain TPC inner params constrained
2383   //
2384   if(!trackInnerC) return kFALSE; 
2385   if(!vtx) return kFALSE;
2386
2387   const Double_t kRadius  = 2.8; 
2388   const Double_t kMaxStep = 1.0; 
2389
2390   Double_t dz[2],cov[3];
2391
2392   //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2393   //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; 
2394   //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; 
2395
2396   if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
2397   if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; 
2398
2399   //
2400   Double_t covar[6]; vtx->GetCovMatrix(covar);
2401   Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
2402   Double_t c[3]={covar[2],0.,covar[5]};
2403   Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
2404   if (chi2C>kVeryBig) return kFALSE; 
2405
2406   if(!trackInnerC->Update(p,c)) return kFALSE;
2407
2408   return kTRUE;
2409 }
2410
2411
2412 //_____________________________________________________________________________
2413 TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack) 
2414 {
2415   if(!particle) return NULL;
2416   if(!stack) return NULL;
2417
2418   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  
2419   TParticle* mother = NULL; 
2420   Int_t mcStackSize=stack->GetNtrack();
2421   if (motherLabel>=mcStackSize) return NULL;
2422   mother = stack->Particle(motherLabel); 
2423
2424   return mother;
2425 }
2426
2427 //_____________________________________________________________________________
2428 Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack) 
2429 {
2430   Bool_t isFromConversion = kFALSE;
2431
2432   if(stack) {
2433     Int_t mcStackSize=stack->GetNtrack();
2434     if (label>=mcStackSize) return kFALSE;
2435     TParticle* particle = stack->Particle(label);
2436     if (!particle) return kFALSE;
2437
2438     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) 
2439     {
2440       Int_t mech = particle->GetUniqueID(); // production mechanism 
2441       Bool_t isPrim = stack->IsPhysicalPrimary(label);
2442
2443       Int_t motherLabel = TMath::Abs(particle->GetMother(0));  
2444       if (motherLabel>=mcStackSize) return kFALSE;
2445       TParticle* mother = stack->Particle(motherLabel); 
2446       if(mother) {
2447         Int_t motherPdg = mother->GetPdgCode();
2448
2449         if(!isPrim && mech==5 && motherPdg==kGamma) { 
2450           isFromConversion=kTRUE; 
2451         }
2452       }
2453     } 
2454   } 
2455
2456   return isFromConversion;
2457 }
2458
2459 //_____________________________________________________________________________
2460 Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack) 
2461 {
2462   Bool_t isFromMaterial = kFALSE;
2463
2464   if(stack) {
2465     Int_t mcStackSize=stack->GetNtrack();
2466     if (label>=mcStackSize) return kFALSE;
2467     TParticle* particle = stack->Particle(label);
2468     if (!particle) return kFALSE;
2469
2470     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) 
2471     {
2472       Int_t mech = particle->GetUniqueID(); // production mechanism 
2473       Bool_t isPrim = stack->IsPhysicalPrimary(label);
2474
2475       Int_t motherLabel = TMath::Abs(particle->GetMother(0));  
2476       if (motherLabel>=mcStackSize) return kFALSE;
2477       TParticle* mother = stack->Particle(motherLabel); 
2478       if(mother) {
2479         if(!isPrim && mech==13) { 
2480           isFromMaterial=kTRUE; 
2481         }
2482       }
2483     } 
2484   } 
2485
2486   return isFromMaterial;
2487 }
2488
2489 //_____________________________________________________________________________
2490 Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack) 
2491 {
2492   Bool_t isFromStrangeness = kFALSE;
2493
2494   if(stack) {
2495     Int_t mcStackSize=stack->GetNtrack();
2496     if (label>=mcStackSize) return kFALSE;
2497     TParticle* particle = stack->Particle(label);
2498     if (!particle) return kFALSE;
2499
2500     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) 
2501     {
2502       Int_t mech = particle->GetUniqueID(); // production mechanism 
2503       Bool_t isPrim = stack->IsPhysicalPrimary(label);
2504
2505       Int_t motherLabel = TMath::Abs(particle->GetMother(0));  
2506       if (motherLabel>=mcStackSize) return kFALSE;
2507       TParticle* mother = stack->Particle(motherLabel); 
2508       if(mother) {
2509         Int_t motherPdg = mother->GetPdgCode();
2510
2511         // K+-, lambda, antilambda, K0s decays
2512         if(!isPrim && mech==4 && 
2513             (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
2514         {
2515           isFromStrangeness = kTRUE;
2516         } 
2517       }
2518     } 
2519   } 
2520
2521   return isFromStrangeness;
2522 }
2523
2524
2525 //_____________________________________________________________________________
2526 void AliAnalysisTaskFilteredTree::FinishTaskOutput() 
2527 {
2528   //
2529   // Called one at the end 
2530   // locally on working node
2531   //
2532
2533   //// must be deleted to store trees
2534   //if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;
2535
2536   //// open temporary file and copy trees to the ouptut container
2537
2538   //TChain* chain = 0;
2539   ////
2540   //chain = new TChain("highPt");
2541   //if(chain) { 
2542   //  chain->Add("jotwinow_Temp_Trees.root");
2543   //  fHighPtTree = chain->CopyTree("1");
2544   //  delete chain; chain=0; 
2545   //}
2546   //if(fHighPtTree) fHighPtTree->Print();
2547
2548   ////
2549   //chain = new TChain("V0s");
2550   //if(chain) { 
2551   //  chain->Add("jotwinow_Temp_Trees.root");
2552   //  fV0Tree = chain->CopyTree("1");
2553   //  delete chain; chain=0; 
2554   //}
2555   //if(fV0Tree) fV0Tree->Print();
2556
2557   ////
2558   //chain = new TChain("dEdx");
2559   //if(chain) { 
2560   //  chain->Add("jotwinow_Temp_Trees.root");
2561   //  fdEdxTree = chain->CopyTree("1");
2562   //  delete chain; chain=0; 
2563   //}
2564   //if(fdEdxTree) fdEdxTree->Print();
2565
2566   ////
2567   //chain = new TChain("Laser");
2568   //if(chain) { 
2569   //  chain->Add("jotwinow_Temp_Trees.root");
2570   //  fLaserTree = chain->CopyTree("1");
2571   //  delete chain; chain=0; 
2572   //}
2573   //if(fLaserTree) fLaserTree->Print();
2574
2575   ////
2576   //chain = new TChain("MCEffTree");
2577   //if(chain) { 
2578   //  chain->Add("jotwinow_Temp_Trees.root");
2579   //  fMCEffTree = chain->CopyTree("1");
2580   //  delete chain; chain=0; 
2581   //}
2582   //if(fMCEffTree) fMCEffTree->Print();
2583
2584   ////
2585   //chain = new TChain("CosmicPairs");
2586   //if(chain) { 
2587   //  chain->Add("jotwinow_Temp_Trees.root");
2588   //  fCosmicPairsTree = chain->CopyTree("1");
2589   //  delete chain; chain=0; 
2590   //}
2591   //if(fCosmicPairsTree) fCosmicPairsTree->Print();  
2592
2593   Bool_t deleteTrees=kTRUE;
2594   if ((AliAnalysisManager::GetAnalysisManager()))
2595   {
2596     if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == 
2597         AliAnalysisManager::kProofAnalysis)
2598       deleteTrees=kFALSE;
2599   }
2600   if (deleteTrees) delete fTreeSRedirector;
2601   fTreeSRedirector=NULL;
2602
2603   //OpenFile(1);
2604
2605   // Post output data.
2606   //PostData(1, fHighPtTree);
2607   //PostData(2, fV0Tree);
2608   //PostData(3, fdEdxTree);
2609   //PostData(4, fLaserTree);
2610   //PostData(5, fMCEffTree);
2611   //PostData(6, fCosmicPairsTree);
2612 }
2613
2614 //_____________________________________________________________________________
2615 void AliAnalysisTaskFilteredTree::Terminate(Option_t *) 
2616 {
2617   // Called one at the end 
2618   /*
2619      fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));
2620      if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;
2621      TChain* chain = new TChain("highPt");
2622      if(!chain) return;
2623      chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");
2624      TTree *tree = chain->CopyTree("1");
2625      if (chain) { delete chain; chain=0; }
2626      if(!tree) return;
2627      tree->Print();
2628      fOutputSummary = tree;
2629
2630      if (!fOutputSummary) {
2631      Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1));
2632      return;
2633      }
2634      */
2635
2636 }
2637
2638 //_____________________________________________________________________________
2639 Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
2640 {
2641   //
2642   // calculate mc event true track multiplicity
2643   //
2644   if(!mcEvent) return 0;
2645
2646   AliStack* stack = 0;
2647   Int_t mult = 0;
2648
2649   // MC particle stack
2650   stack = mcEvent->Stack();
2651   if (!stack) return 0;
2652
2653   //
2654   //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
2655   //
2656
2657   Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
2658   if(!isEventOK) return 0; 
2659
2660   Int_t nPart  = stack->GetNtrack();
2661   for (Int_t iMc = 0; iMc < nPart; ++iMc) 
2662   {
2663     TParticle* particle = stack->Particle(iMc);
2664     if (!particle)
2665       continue;
2666
2667     // only charged particles
2668     if(!particle->GetPDG()) continue;
2669     Double_t charge = particle->GetPDG()->Charge()/3.;
2670     if (TMath::Abs(charge) < 0.001)
2671       continue;
2672
2673     // physical primary
2674     Bool_t prim = stack->IsPhysicalPrimary(iMc);
2675     if(!prim) continue;
2676
2677     // checked accepted without pt cut
2678     //if(accCuts->AcceptTrack(particle)) 
2679     if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) 
2680     {
2681       mult++;
2682     }
2683   }
2684
2685   return mult;  
2686 }
2687
2688 //_____________________________________________________________________________
2689 void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, const Double_t centralityF, const Double_t chi2TPCInnerC) 
2690 {
2691   //
2692   // Fill pT relative resolution histograms for 
2693   // TPC only, TPC only constrained to vertex and TPC+ITS tracking
2694   //
2695   if(!ptrack) return;    
2696   if(!ptpcInnerC) return;    
2697
2698   const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();
2699   if(!innerParam) return;
2700
2701   Float_t dxy, dz;
2702   ptrack->GetImpactParameters(dxy,dz);
2703
2704   // TPC+ITS primary tracks 
2705   if( abs(ptrack->Eta())<0.8 && 
2706       ptrack->GetTPCClusterInfo(3,1)>120 && 
2707       ptrack->IsOn(0x40) && 
2708       ptrack->GetTPCclusters(0)>0.0 &&  
2709       ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 && 
2710       //abs(innerParam->GetX())>0.0 && 
2711       //abs(innerParam->GetY()/innerParam->GetX())<0.14 && 
2712       //abs(innerParam->GetTgl())<0.85 && 
2713       ptrack->IsOn(0x0004) && 
2714       ptrack->GetNcls(0)>0 &&
2715       ptrack->GetITSchi2()>0 && 
2716       sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 &&
2717       sqrt(chi2TPCInnerC)<6 &&
2718       (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) &&
2719       abs(dz)<2.0 && 
2720       abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )
2721   {
2722     fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2723     fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2724     fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2725   }
2726
2727   // TPC primary tracks 
2728   // and TPC constrained primary tracks 
2729
2730   AliExternalTrackParam *ptpcInner  = (AliExternalTrackParam *) ptrack->GetTPCInnerParam(); 
2731   if(!ptpcInner) return;
2732
2733
2734   Float_t dxyTPC, dzTPC;
2735   ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);
2736
2737   if( abs(ptrack->Eta())<0.8 && 
2738       ptrack->GetTPCClusterInfo(3,1)>120 && 
2739       ptrack->IsOn(0x40)&& 
2740       ptrack->GetTPCclusters(0)>0.0 &&  
2741       ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 && 
2742       //abs(innerParam->GetX())>0.0 && 
2743       //abs(innerParam->GetY()/innerParam->GetX())<0.14 && 
2744       //abs(innerParam->GetTgl())<0.85 && 
2745       abs(dzTPC)<3.2 && 
2746       abs(dxyTPC)<2.4 )
2747   {
2748     // TPC only
2749     fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2750     fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2751     fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2752
2753     // TPC constrained to vertex 
2754     fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2755     fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2756     fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2757   }
2758 }