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