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