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