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