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