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