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