modification in task to filter events
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / AlidNdPtTrackDumpTask.cxx
CommitLineData
a26e43aa 1/**************************************************************************\r
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3* *\r
4* Author: The ALICE Off-line Project. *\r
5* Contributors are mentioned in the code where appropriate. *\r
6* *\r
7* Permission to use, copy, modify and distribute this software and its *\r
8* documentation strictly for non-commercial purposes is hereby granted *\r
9* without fee, provided that the above copyright notice appears in all *\r
10* copies and that both the copyright notice and this permission notice *\r
11* appear in the supporting documentation. The authors make no claims *\r
12* about the suitability of this software for any purpose. It is *\r
13* provided "as is" without express or implied warranty. *\r
14**************************************************************************/\r
15\r
16#include "iostream"\r
17\r
18#include <TPDGCode.h>\r
91326969 19#include <TDatabasePDG.h>\r
a26e43aa 20\r
21#include "TChain.h"\r
22#include "TTreeStream.h"\r
23#include "TTree.h"\r
24#include "TH1F.h"\r
25#include "TCanvas.h"\r
26#include "TList.h"\r
99496190 27#include "TObjArray.h"\r
a26e43aa 28#include "TFile.h"\r
29#include "TMatrixD.h"\r
91326969 30#include "TRandom3.h"\r
a26e43aa 31\r
32#include "AliHeader.h" \r
33#include "AliGenEventHeader.h" \r
34#include "AliInputEventHandler.h" \r
35#include "AliStack.h" \r
36#include "AliTrackReference.h" \r
37\r
38#include "AliPhysicsSelection.h"\r
39#include "AliAnalysisTask.h"\r
40#include "AliAnalysisManager.h"\r
41#include "AliESDEvent.h"\r
42#include "AliESDfriend.h"\r
43#include "AliMCEvent.h"\r
44#include "AliESDInputHandler.h"\r
45#include "AliESDVertex.h"\r
46#include "AliTracker.h"\r
47#include "AliGeomManager.h"\r
48\r
49#include "AliCentrality.h"\r
50#include "AliESDVZERO.h"\r
51#include "AliMultiplicity.h"\r
52\r
53#include "AliESDtrackCuts.h"\r
54#include "AliMCEventHandler.h"\r
bdd49ee6 55#include "AlidNdPt.h"\r
56#include "AlidNdPtEventCuts.h"\r
57#include "AlidNdPtAcceptanceCuts.h"\r
a26e43aa 58\r
bdd49ee6 59#include "AlidNdPtTrackDumpTask.h"\r
91326969 60#include "AliKFParticle.h"\r
61#include "AliESDv0.h"\r
a26e43aa 62\r
63using namespace std;\r
64\r
65ClassImp(AlidNdPtTrackDumpTask)\r
66\r
67//_____________________________________________________________________________\r
68AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name) \r
69 : AliAnalysisTaskSE(name)\r
70 , fESD(0)\r
71 , fMC(0)\r
72 , fESDfriend(0)\r
73 , fOutput(0)\r
74 , fPitList(0)\r
75 , fUseMCInfo(kFALSE)\r
76 , fdNdPtEventCuts(0)\r
77 , fdNdPtAcceptanceCuts(0)\r
78 , fdNdPtRecAcceptanceCuts(0)\r
79 , fEsdTrackCuts(0)\r
80 , fTrigger(AliTriggerAnalysis::kMB1) \r
81 , fAnalysisMode(AlidNdPtHelper::kTPC) \r
a26e43aa 82 , fTreeSRedirector(0)\r
83 , fCentralityEstimator(0)\r
91326969 84 , fLowPtTrackDownscaligF(0)\r
85 , fLowPtV0DownscaligF(0)\r
298fa346 86 , fProcessAll(kFALSE)\r
a26e43aa 87{\r
88 // Constructor\r
89\r
90 // Define input and output slots here\r
298fa346 91 DefineOutput(1, TList::Class());\r
a26e43aa 92\r
93}\r
94\r
95//_____________________________________________________________________________\r
96AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()\r
97{\r
98 if(fOutput) delete fOutput; fOutput =0; \r
a26e43aa 99 if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector =0; \r
100\r
101 if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL; \r
102 if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;\r
103 if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL; \r
104 if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;\r
105}\r
106\r
107//____________________________________________________________________________\r
108Bool_t AlidNdPtTrackDumpTask::Notify()\r
109{\r
110 static Int_t count = 0;\r
111 count++;\r
112 TTree *chain = (TChain*)GetInputData(0);\r
113 if(chain)\r
114 {\r
115 Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());\r
116 }\r
117\r
118return kTRUE;\r
119}\r
120\r
121//_____________________________________________________________________________\r
122void AlidNdPtTrackDumpTask::UserCreateOutputObjects()\r
123{\r
124 // Create histograms\r
125 // Called once\r
126 fOutput = new TList;\r
127 fOutput->SetOwner();\r
128\r
129 //\r
99496190 130 // create temporary file for output tree\r
131 fTreeSRedirector = new TTreeSRedirector("jotwinow_Temp_Trees.root");\r
a26e43aa 132\r
298fa346 133 PostData(1, fOutput);\r
a26e43aa 134}\r
135\r
136//_____________________________________________________________________________\r
137void AlidNdPtTrackDumpTask::UserExec(Option_t *) \r
138{\r
139 //\r
140 // Called for each event\r
141 //\r
142\r
143 // ESD event\r
144 fESD = (AliESDEvent*) (InputEvent());\r
145 if (!fESD) {\r
146 Printf("ERROR: ESD event not available");\r
147 return;\r
148 }\r
149\r
150 // MC event\r
151 if(fUseMCInfo) {\r
152 fMC = MCEvent();\r
153 if (!fMC) {\r
154 Printf("ERROR: MC event not available");\r
155 return;\r
156 }\r
157 }\r
158\r
159 fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
160 if(!fESDfriend) {\r
161 Printf("ERROR: ESD friends not available");\r
162 }\r
163\r
164 //\r
298fa346 165 if(fProcessAll) { \r
166 ProcessAll(fESD,fMC,fESDfriend);\r
167 ProcessV0(fESD,fMC,fESDfriend);\r
168 }\r
169 else {\r
170 Process(fESD,fMC,fESDfriend);\r
171 ProcessV0(fESD,fMC,fESDfriend);\r
172 }\r
a26e43aa 173\r
174 // Post output data.\r
298fa346 175 PostData(1, fOutput);\r
a26e43aa 176}\r
177\r
178//_____________________________________________________________________________\r
298fa346 179void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
180{\r
181 //\r
182 // Process real and/or simulated events\r
183 // Only filtering\r
184 //\r
185 if(!esdEvent) {\r
186 AliDebug(AliLog::kError, "esdEvent not available");\r
187 return;\r
188 }\r
189\r
190 // get selection cuts\r
191 AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
192 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
193 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
194\r
195 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
196 AliDebug(AliLog::kError, "cuts not available");\r
197 return;\r
198 }\r
199\r
200 // trigger selection\r
201 Bool_t isEventTriggered = kTRUE;\r
202 AliPhysicsSelection *physicsSelection = NULL;\r
203 AliTriggerAnalysis* triggerAnalysis = NULL;\r
204\r
205 // \r
206 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
207 if (!inputHandler)\r
208 {\r
209 Printf("ERROR: Could not receive input handler");\r
210 return;\r
211 }\r
212 \r
213 // get file name\r
214 TTree *chain = (TChain*)GetInputData(0);\r
215 if(!chain) { \r
216 Printf("ERROR: Could not receive input chain");\r
217 return;\r
218 }\r
219 TObjString fileName(chain->GetCurrentFile()->GetName());\r
220\r
221 // trigger\r
222 if(evtCuts->IsTriggerRequired()) \r
223 {\r
224 // always MB\r
225 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
226\r
227 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
228 if(!physicsSelection) return;\r
229 //SetPhysicsTriggerSelection(physicsSelection);\r
230\r
231 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
232 // set trigger (V0AND)\r
233 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
234 if(!triggerAnalysis) return;\r
235 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
236 }\r
237 }\r
238\r
239 // centrality determination\r
240 Float_t centralityF = -1;\r
241 AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
242 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
243\r
244 // use MC information\r
245 AliHeader* header = 0;\r
246 AliGenEventHeader* genHeader = 0;\r
247 AliStack* stack = 0;\r
248 TArrayF vtxMC(3);\r
249\r
250 Int_t multMCTrueTracks = 0;\r
251 if(IsUseMCInfo())\r
252 {\r
253 //\r
254 if(!mcEvent) {\r
255 AliDebug(AliLog::kError, "mcEvent not available");\r
256 return;\r
257 }\r
258 // get MC event header\r
259 header = mcEvent->Header();\r
260 if (!header) {\r
261 AliDebug(AliLog::kError, "Header not available");\r
262 return;\r
263 }\r
264 // MC particle stack\r
265 stack = mcEvent->Stack();\r
266 if (!stack) {\r
267 AliDebug(AliLog::kError, "Stack not available");\r
268 return;\r
269 }\r
270\r
271 // get MC vertex\r
272 genHeader = header->GenEventHeader();\r
273 if (!genHeader) {\r
274 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
275 return;\r
276 }\r
277 genHeader->PrimaryVertex(vtxMC);\r
278\r
279 // multipliticy of all MC primary tracks\r
280 // in Zv, pt and eta ranges)\r
281 multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
282\r
283 } // end bUseMC\r
284\r
285 // laser events \r
286 if(esdEvent)\r
287 {\r
288 const AliESDHeader* esdHeader = esdEvent->GetHeader();\r
289 if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) \r
290 {\r
291 Int_t countLaserTracks = 0;\r
292 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
293 {\r
294 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
295 if(!track) continue;\r
296\r
297 if(track->GetTPCInnerParam()) countLaserTracks++;\r
298 }\r
299 \r
300 if(countLaserTracks > 100) { \r
301\r
302 Double_t runNumber = esdEvent->GetRunNumber();\r
303 Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
304 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
305\r
306 if(!fTreeSRedirector) return;\r
307 (*fTreeSRedirector)<<"Laser"<<\r
308 "fileName.="<<&fileName<<\r
309 "runNumber="<<runNumber<<\r
310 "evtTimeStamp="<<evtTimeStamp<<\r
311 "evtNumberInFile="<<evtNumberInFile<<\r
312 "multTPCtracks="<<countLaserTracks<<\r
313 "\n";\r
314 }\r
315 }\r
316 }\r
317\r
318\r
319 // get reconstructed vertex \r
320 //const AliESDVertex* vtxESD = 0; \r
321 const AliESDVertex* vtxESD = 0; \r
322 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
323 vtxESD = esdEvent->GetPrimaryVertexTPC();\r
324 }\r
325 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
326 vtxESD = esdEvent->GetPrimaryVertexTracks();\r
327 }\r
328 else {\r
329 return;\r
330 }\r
331\r
332 if(!vtxESD) return;\r
333\r
334 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
335 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
336 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
337\r
338\r
339 // check event cuts\r
340 if(isEventOK && isEventTriggered)\r
341 {\r
342 //\r
343 Double_t vert[3] = {0}; \r
344 vert[0] = vtxESD->GetXv();\r
345 vert[1] = vtxESD->GetYv();\r
346 vert[2] = vtxESD->GetZv();\r
347 Int_t mult = vtxESD->GetNContributors();\r
348 Double_t bz = esdEvent->GetMagneticField();\r
349 Double_t runNumber = esdEvent->GetRunNumber();\r
350 Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
351 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
352\r
353 // high pT tracks\r
354 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
355 {\r
356 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
357 if(!track) continue;\r
358 if(track->Charge()==0) continue;\r
359 if(!esdTrackCuts->AcceptTrack(track)) continue;\r
360 if(!accCuts->AcceptTrack(track)) continue;\r
361 \r
362 // downscale low-pT tracks\r
363 Double_t scalempt= TMath::Min(track->Pt(),10.);\r
364 Double_t downscaleF = gRandom->Rndm();\r
365 downscaleF *= fLowPtTrackDownscaligF;\r
366 if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
367 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
368\r
369 // Dump to the tree \r
370 // vertex\r
371 // TPC-ITS tracks\r
372 //\r
373 if(!fTreeSRedirector) return;\r
374 (*fTreeSRedirector)<<"dNdPtTree"<<\r
375 "fileName.="<<&fileName<<\r
376 "runNumber="<<runNumber<<\r
377 "evtTimeStamp="<<evtTimeStamp<<\r
378 "evtNumberInFile="<<evtNumberInFile<<\r
379 "Bz="<<bz<<\r
380 "vertX="<<vert[0]<<\r
381 "vertY="<<vert[1]<<\r
382 "vertZ="<<vert[2]<<\r
383 "mult="<<mult<<\r
384 "esdTrack.="<<track<<\r
385 "centralityF="<<centralityF<<\r
386 "\n";\r
387 }\r
388\r
389 // high dEdx\r
390 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
391 {\r
392 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
393 if(!track) continue;\r
394 if(track->Charge()==0) continue;\r
395 if(!esdTrackCuts->AcceptTrack(track)) continue;\r
396 if(!accCuts->AcceptTrack(track)) continue;\r
397\r
398 if(!IsHighDeDxParticle(track)) continue;\r
399 \r
400 if(!fTreeSRedirector) return;\r
401 (*fTreeSRedirector)<<"dEdx"<<\r
402 "fileName.="<<&fileName<<\r
403 "runNumber="<<runNumber<<\r
404 "evtTimeStamp="<<evtTimeStamp<<\r
405 "evtNumberInFile="<<evtNumberInFile<<\r
406 "Bz="<<bz<<\r
407 "vertX="<<vert[0]<<\r
408 "vertY="<<vert[1]<<\r
409 "vertZ="<<vert[2]<<\r
410 "mult="<<mult<<\r
411 "esdTrack.="<<track<<\r
412 "\n";\r
413 }\r
414 }\r
415 \r
416 PostData(1, fOutput);\r
417}\r
418\r
419\r
420\r
421\r
422\r
423\r
424\r
425//_____________________________________________________________________________\r
426void AlidNdPtTrackDumpTask::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
a26e43aa 427{\r
428 //\r
429 // Process real and/or simulated events\r
430 //\r
431 if(!esdEvent) {\r
432 AliDebug(AliLog::kError, "esdEvent not available");\r
433 return;\r
434 }\r
435\r
436 // get selection cuts\r
437 AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
438 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
439 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
440\r
441 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
442 AliDebug(AliLog::kError, "cuts not available");\r
443 return;\r
444 }\r
445\r
446 // trigger selection\r
447 Bool_t isEventTriggered = kTRUE;\r
448 AliPhysicsSelection *physicsSelection = NULL;\r
449 AliTriggerAnalysis* triggerAnalysis = NULL;\r
450\r
451 // \r
452 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
453 if (!inputHandler)\r
454 {\r
455 Printf("ERROR: Could not receive input handler");\r
456 return;\r
457 }\r
bdd49ee6 458 \r
459 // get file name\r
460 TTree *chain = (TChain*)GetInputData(0);\r
461 if(!chain) { \r
462 Printf("ERROR: Could not receive input chain");\r
463 return;\r
464 }\r
465 TObjString fileName(chain->GetCurrentFile()->GetName());\r
a26e43aa 466\r
bdd49ee6 467 // trigger\r
a26e43aa 468 if(evtCuts->IsTriggerRequired()) \r
469 {\r
470 // always MB\r
471 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
472\r
473 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
474 if(!physicsSelection) return;\r
475 //SetPhysicsTriggerSelection(physicsSelection);\r
476\r
477 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
478 // set trigger (V0AND)\r
479 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
480 if(!triggerAnalysis) return;\r
481 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
482 }\r
483 }\r
484\r
485 // centrality determination\r
486 Float_t centralityF = -1;\r
487 AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
488 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
489\r
490 // use MC information\r
491 AliHeader* header = 0;\r
492 AliGenEventHeader* genHeader = 0;\r
493 AliStack* stack = 0;\r
494 TArrayF vtxMC(3);\r
495\r
496 Int_t multMCTrueTracks = 0;\r
497 if(IsUseMCInfo())\r
498 {\r
499 //\r
500 if(!mcEvent) {\r
501 AliDebug(AliLog::kError, "mcEvent not available");\r
502 return;\r
503 }\r
504 // get MC event header\r
505 header = mcEvent->Header();\r
506 if (!header) {\r
507 AliDebug(AliLog::kError, "Header not available");\r
508 return;\r
509 }\r
510 // MC particle stack\r
511 stack = mcEvent->Stack();\r
512 if (!stack) {\r
513 AliDebug(AliLog::kError, "Stack not available");\r
514 return;\r
515 }\r
516\r
517 // get MC vertex\r
518 genHeader = header->GenEventHeader();\r
519 if (!genHeader) {\r
520 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
521 return;\r
522 }\r
523 genHeader->PrimaryVertex(vtxMC);\r
524\r
525 // multipliticy of all MC primary tracks\r
526 // in Zv, pt and eta ranges)\r
527 multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
528\r
529 } // end bUseMC\r
530\r
99496190 531 // laser events \r
532 if(esdEvent)\r
533 {\r
534 const AliESDHeader* esdHeader = esdEvent->GetHeader();\r
535 if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) \r
536 {\r
537 Int_t countLaserTracks = 0;\r
538 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
539 {\r
540 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
541 if(!track) continue;\r
542\r
543 if(track->GetTPCInnerParam()) countLaserTracks++;\r
544 }\r
545 \r
546 if(countLaserTracks > 100) { \r
547\r
548 Double_t runNumber = esdEvent->GetRunNumber();\r
549 Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
550 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
551\r
552 if(!fTreeSRedirector) return;\r
553 (*fTreeSRedirector)<<"Laser"<<\r
554 "fileName.="<<&fileName<<\r
555 "runNumber="<<runNumber<<\r
556 "evtTimeStamp="<<evtTimeStamp<<\r
557 "evtNumberInFile="<<evtNumberInFile<<\r
558 "multTPCtracks="<<countLaserTracks<<\r
559 "\n";\r
560 }\r
561 }\r
562 }\r
563\r
564\r
a26e43aa 565 // get reconstructed vertex \r
566 //const AliESDVertex* vtxESD = 0; \r
567 const AliESDVertex* vtxESD = 0; \r
568 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
569 vtxESD = esdEvent->GetPrimaryVertexTPC();\r
570 }\r
571 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
572 vtxESD = esdEvent->GetPrimaryVertexTracks();\r
573 }\r
574 else {\r
575 return;\r
576 }\r
577\r
578 if(!vtxESD) return;\r
579\r
580 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
581 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
582 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
583\r
99496190 584\r
a26e43aa 585 // check event cuts\r
586 if(isEventOK && isEventTriggered)\r
587 {\r
99496190 588 //\r
589 Double_t vert[3] = {0}; \r
590 vert[0] = vtxESD->GetXv();\r
591 vert[1] = vtxESD->GetYv();\r
592 vert[2] = vtxESD->GetZv();\r
593 Int_t mult = vtxESD->GetNContributors();\r
594 Double_t bz = esdEvent->GetMagneticField();\r
595 Double_t runNumber = esdEvent->GetRunNumber();\r
596 Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
597 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
598\r
599 // high pT tracks\r
a26e43aa 600 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
601 {\r
602 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
603 if(!track) continue;\r
604 if(track->Charge()==0) continue;\r
605 if(!esdTrackCuts->AcceptTrack(track)) continue;\r
606 if(!accCuts->AcceptTrack(track)) continue;\r
99496190 607 \r
08f8415e 608 // downscale low-pT tracks\r
91326969 609 Double_t scalempt= TMath::Min(track->Pt(),10.);\r
99496190 610 Double_t downscaleF = gRandom->Rndm();\r
611 downscaleF *= fLowPtTrackDownscaligF;\r
612 if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
613 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
08f8415e 614\r
a26e43aa 615 // Dump to the tree \r
616 // vertex\r
617 // TPC constrained tracks\r
618 // InnerParams constrained tracks\r
619 // TPC-ITS tracks\r
620 // ITSout-InnerParams tracks\r
621 // chi2 distance between TPC constrained and TPC-ITS tracks\r
622 // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
623 // chi2 distance between ITSout and InnerParams tracks\r
624 // MC information\r
625 \r
626 Double_t x[3]; track->GetXYZ(x);\r
627 Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
628 Bool_t isOK = kFALSE;\r
629\r
630 //\r
631 // Constrain TPC-only tracks (TPCinner) to vertex\r
632 //\r
633 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
634 if (!tpcInner) continue;\r
635 // transform to the track reference frame \r
636 isOK = tpcInner->Rotate(track->GetAlpha());\r
637 isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
638 if(!isOK) continue;\r
639\r
640 // clone TPCinner has to be deleted\r
458d14c4 641 AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
a26e43aa 642 if (!tpcInnerC) continue;\r
643 \r
644 // constrain TPCinner \r
645 //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());\r
646 isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
647\r
648 // transform to the track reference frame \r
649 isOK = tpcInnerC->Rotate(track->GetAlpha());\r
650 isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
651\r
652 if(!isOK) {\r
653 if(tpcInnerC) delete tpcInnerC;\r
654 continue;\r
655 }\r
656\r
657\r
658 //\r
659 // Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
660 // to vertex\r
661 //\r
662 // clone track InnerParams has to be deleted\r
458d14c4 663 AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));\r
a26e43aa 664 if (!trackInnerC) continue;\r
665 \r
666 // constrain track InnerParams \r
667 isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
668\r
669 // transform to the track reference frame \r
670 isOK = trackInnerC->Rotate(track->GetAlpha());\r
671 isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
672\r
673 if(!isOK) {\r
674 if(trackInnerC) delete trackInnerC;\r
675 continue;\r
676 }\r
677 \r
678 //\r
679 // calculate chi2 between vi and vj vectors\r
680 // with covi and covj covariance matrices\r
681 // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
682 //\r
683 TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
684 TMatrixD delta(1,5), deltatrackC(1,5);\r
685 TMatrixD covarM(5,5), covarMtrackC(5,5);\r
686\r
687 for (Int_t ipar=0; ipar<5; ipar++) {\r
688 deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
689 delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
690\r
691 deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
692 deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
693\r
694 for (Int_t jpar=0; jpar<5; jpar++) {\r
695 Int_t index=track->GetIndex(ipar,jpar);\r
696 covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
697 covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
698 }\r
699 }\r
700 // chi2 distance TPC constrained and TPC+ITS\r
701 TMatrixD covarMInv = covarM.Invert();\r
702 TMatrixD mat2 = covarMInv*deltaT;\r
703 TMatrixD chi2 = delta*mat2; \r
704 //chi2.Print();\r
705\r
706 // chi2 distance TPC refitted constrained and TPC+ITS\r
707 TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
708 TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
709 TMatrixD chi2trackC = deltatrackC*mat2trackC; \r
710 //chi2trackC.Print();\r
711\r
712\r
713 //\r
714 // Propagate ITSout to TPC inner wall \r
715 // and calculate chi2 distance to track (InnerParams)\r
716 //\r
717 const Double_t kTPCRadius=85; \r
718 const Double_t kStep=3; \r
719\r
720 // clone track InnerParams has to be deleted\r
458d14c4 721 AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
a26e43aa 722 if (!trackInnerC2) continue;\r
723 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
724 {\r
725 if(trackInnerC2) delete trackInnerC2;\r
726 continue;\r
727 }\r
728\r
729 AliExternalTrackParam *outerITSc = new AliExternalTrackParam();\r
730 if(!outerITSc) continue;\r
731\r
732 TMatrixD chi2OuterITS(1,1);\r
733 chi2OuterITS(0,0) = 0;\r
734\r
735 if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
736 {\r
737 // propagate ITSout to TPC inner wall\r
738 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
739\r
740 if(friendTrack) \r
741 {\r
458d14c4 742 if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) ) \r
a26e43aa 743 {\r
744 if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
745 {\r
746 // transform ITSout to the track InnerParams reference frame \r
747 isOK = kFALSE;\r
748 isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
749 isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
750\r
751 if(!isOK) {\r
752 if(outerITSc) delete outerITSc;\r
753 if(trackInnerC2) delete trackInnerC2;\r
754 continue;\r
755 }\r
756 \r
757 //\r
758 // calculate chi2 between outerITS and innerParams\r
759 // cov without z-coordinate at the moment\r
760 //\r
761 TMatrixD deltaTouterITS(4,1);\r
762 TMatrixD deltaouterITS(1,4);\r
763 TMatrixD covarMouterITS(4,4);\r
764\r
765 Int_t kipar = 0;\r
766 Int_t kjpar = 0;\r
767 for (Int_t ipar=0; ipar<5; ipar++) {\r
768 if(ipar!=1) {\r
769 deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
770 deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
771 }\r
772\r
773 kjpar=0;\r
774 for (Int_t jpar=0; jpar<5; jpar++) {\r
775 Int_t index=outerITSc->GetIndex(ipar,jpar);\r
776 if(ipar !=1 || jpar!=1) {\r
777 covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
778 }\r
779 if(jpar!=1) kjpar++;\r
780 }\r
781 if(ipar!=1) kipar++;\r
782 }\r
783\r
784 // chi2 distance ITSout and InnerParams\r
785 TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
786 TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
787 chi2OuterITS = deltaouterITS*mat2outerITS; \r
788 //chi2OuterITS.Print();\r
789 } \r
790 }\r
791 }\r
792 }\r
793\r
794 //\r
795 // MC info\r
796 //\r
797 TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
798 TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
799 Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
800 Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
801 Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
802 Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
803 Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
804\r
805 AliTrackReference *refTPCIn = NULL;\r
806 AliTrackReference *refITS = NULL;\r
458d14c4 807 AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
a26e43aa 808 if(!trackInnerC3) continue;\r
809\r
810 if(IsUseMCInfo()) \r
811 {\r
812 if(!stack) return;\r
813\r
814 //\r
815 // global track\r
816 //\r
817 Int_t label = TMath::Abs(track->GetLabel()); \r
818 particle = stack->Particle(label);\r
819 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
820 {\r
821 particleMother = GetMother(particle,stack);\r
822 mech = particle->GetUniqueID();\r
823 isPrim = stack->IsPhysicalPrimary(label);\r
824 isFromStrangess = IsFromStrangeness(label,stack);\r
825 isFromConversion = IsFromConversion(label,stack);\r
826 isFromMaterial = IsFromMaterial(label,stack);\r
827 }\r
828\r
829 //\r
830 // TPC track\r
831 //\r
832 Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
833 particleTPC = stack->Particle(labelTPC);\r
834 if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
835 {\r
836 particleMotherTPC = GetMother(particleTPC,stack);\r
837 mechTPC = particleTPC->GetUniqueID();\r
838 isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
839 isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);\r
840 isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
841 isFromMaterialTPC = IsFromMaterial(labelTPC,stack);\r
842 }\r
843\r
844 //\r
845 // store first track reference\r
846 // for TPC track\r
847 //\r
848 TParticle *part=0;\r
849 TClonesArray *trefs=0;\r
850 Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
851\r
852 if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
853 {\r
854 Int_t nTrackRef = trefs->GetEntries();\r
855 //printf("nTrackRef %d \n",nTrackRef);\r
856\r
857 Int_t countITS = 0;\r
858 for (Int_t iref = 0; iref < nTrackRef; iref++) \r
859 {\r
860 AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
861 //printf("ref %p \n",ref);\r
862 //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
863 //printf("AliTrackReference::kTPC %d \n",AliTrackReference::kTPC);\r
864\r
865\r
866 // Ref. in the middle ITS \r
867 if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
868 {\r
869 if(!refITS && countITS==2) {\r
870 refITS = ref;\r
871 //printf("refITS %p \n",refITS);\r
872 }\r
873 countITS++;\r
874 }\r
875\r
876 // TPC\r
877 if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
878 {\r
879 if(!refTPCIn) {\r
880 refTPCIn = ref;\r
881 //printf("refTPCIn %p \n",refTPCIn);\r
882 //break;\r
883 }\r
884 }\r
885 }\r
886\r
887 // transform inner params to TrackRef\r
888 // reference frame\r
889 if(refTPCIn && trackInnerC3) \r
890 {\r
891 Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
892 isOK = kFALSE;\r
893 isOK = trackInnerC3->Rotate(kRefPhi);\r
894 isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
895\r
896 if(!isOK){\r
897 if(trackInnerC3) delete trackInnerC3;\r
898 if(refTPCIn) delete refTPCIn;\r
899 }\r
900 }\r
901 }\r
902\r
903 //\r
904 // ITS track\r
905 //\r
906 Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
907 particleITS = stack->Particle(labelITS);\r
908 if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
909 {\r
910 particleMotherITS = GetMother(particleITS,stack);\r
911 mechITS = particleITS->GetUniqueID();\r
912 isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
913 isFromStrangessITS = IsFromStrangeness(labelITS,stack);\r
914 isFromConversionITS = IsFromConversion(labelITS,stack);\r
915 isFromMaterialITS = IsFromMaterial(labelITS,stack);\r
916 }\r
917 }\r
918\r
a26e43aa 919 //\r
920 if(!fTreeSRedirector) return;\r
921 (*fTreeSRedirector)<<"dNdPtTree"<<\r
bdd49ee6 922 "fileName.="<<&fileName<<\r
a26e43aa 923 "runNumber="<<runNumber<<\r
08f8415e 924 "evtTimeStamp="<<evtTimeStamp<<\r
bdd49ee6 925 "evtNumberInFile="<<evtNumberInFile<<\r
a26e43aa 926 "Bz="<<bz<<\r
08f8415e 927 "vertX="<<vert[0]<<\r
928 "vertY="<<vert[1]<<\r
929 "vertZ="<<vert[2]<<\r
a26e43aa 930 "mult="<<mult<<\r
931 "esdTrack.="<<track<<\r
932 "extTPCInnerC.="<<tpcInnerC<<\r
933 "extInnerParamC.="<<trackInnerC<<\r
934 "extInnerParam.="<<trackInnerC2<<\r
935 "extOuterITS.="<<outerITSc<<\r
936 "extInnerParamRef.="<<trackInnerC3<<\r
937 "refTPCIn.="<<refTPCIn<<\r
938 "refITS.="<<refITS<<\r
939 "chi2TPCInnerC="<<chi2(0,0)<<\r
940 "chi2InnerC="<<chi2trackC(0,0)<<\r
941 "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
942 "centralityF="<<centralityF<<\r
943 "particle.="<<particle<<\r
944 "particleMother.="<<particleMother<<\r
945 "mech="<<mech<<\r
946 "isPrim="<<isPrim<<\r
947 "isFromStrangess="<<isFromStrangess<<\r
948 "isFromConversion="<<isFromConversion<<\r
949 "isFromMaterial="<<isFromMaterial<<\r
950 "particleTPC.="<<particleTPC<<\r
951 "particleMotherTPC.="<<particleMotherTPC<<\r
952 "mechTPC="<<mechTPC<<\r
953 "isPrimTPC="<<isPrimTPC<<\r
954 "isFromStrangessTPC="<<isFromStrangessTPC<<\r
955 "isFromConversionTPC="<<isFromConversionTPC<<\r
956 "isFromMaterialTPC="<<isFromMaterialTPC<<\r
957 "particleITS.="<<particleITS<<\r
958 "particleMotherITS.="<<particleMotherITS<<\r
959 "mechITS="<<mechITS<<\r
960 "isPrimITS="<<isPrimITS<<\r
961 "isFromStrangessITS="<<isFromStrangessITS<<\r
962 "isFromConversionITS="<<isFromConversionITS<<\r
963 "isFromMaterialITS="<<isFromMaterialITS<<\r
964 "\n";\r
965\r
966 if(tpcInnerC) delete tpcInnerC;\r
967 if(trackInnerC) delete trackInnerC;\r
968 if(trackInnerC2) delete trackInnerC2;\r
969 if(outerITSc) delete outerITSc;\r
970\r
971 if(trackInnerC3) delete trackInnerC3;\r
972 }\r
99496190 973\r
974 // high dEdx\r
975 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
976 {\r
977 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
978 if(!track) continue;\r
979 if(track->Charge()==0) continue;\r
980 if(!esdTrackCuts->AcceptTrack(track)) continue;\r
981 if(!accCuts->AcceptTrack(track)) continue;\r
982\r
983 if(!IsHighDeDxParticle(track)) continue;\r
984 \r
985 if(!fTreeSRedirector) return;\r
986 (*fTreeSRedirector)<<"dEdx"<<\r
987 "fileName.="<<&fileName<<\r
988 "runNumber="<<runNumber<<\r
989 "evtTimeStamp="<<evtTimeStamp<<\r
990 "evtNumberInFile="<<evtNumberInFile<<\r
991 "Bz="<<bz<<\r
992 "vertX="<<vert[0]<<\r
993 "vertY="<<vert[1]<<\r
994 "vertZ="<<vert[2]<<\r
995 "mult="<<mult<<\r
996 "esdTrack.="<<track<<\r
997 "\n";\r
998 }\r
a26e43aa 999 }\r
99496190 1000 \r
298fa346 1001 PostData(1, fOutput);\r
99496190 1002}\r
1003\r
1004//_____________________________________________________________________________\r
1005Bool_t AlidNdPtTrackDumpTask::IsHighDeDxParticle(AliESDtrack * track) {\r
1006 //\r
1007 // check if particle is Z > 1 \r
1008 //\r
1009 if (track->GetTPCNcls() < 60) return kFALSE;\r
1010 Double_t mom = track->GetInnerParam()->GetP();\r
1011 if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization\r
1012 Float_t dca[2], bCov[3];\r
1013 track->GetImpactParameters(dca,bCov);\r
1014 //\r
a26e43aa 1015\r
99496190 1016 Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);\r
1017\r
1018 if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;\r
1019\r
1020 return kFALSE;\r
a26e43aa 1021}\r
1022\r
91326969 1023//_____________________________________________________________________________\r
1024void AlidNdPtTrackDumpTask::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
1025{\r
1026 //\r
1027 // Process real and/or simulated events\r
1028 //\r
1029 if(!esdEvent) {\r
1030 AliDebug(AliLog::kError, "esdEvent not available");\r
1031 return;\r
1032 }\r
1033\r
1034 // get selection cuts\r
1035 AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
1036 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1037 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1038\r
1039 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
1040 AliDebug(AliLog::kError, "cuts not available");\r
1041 return;\r
1042 }\r
1043\r
91326969 1044 // trigger selection\r
1045 Bool_t isEventTriggered = kTRUE;\r
1046 AliPhysicsSelection *physicsSelection = NULL;\r
1047 AliTriggerAnalysis* triggerAnalysis = NULL;\r
1048\r
1049 // \r
1050 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
1051 if (!inputHandler)\r
1052 {\r
1053 Printf("ERROR: Could not receive input handler");\r
1054 return;\r
1055 }\r
1056 \r
1057 // get file name\r
1058 TTree *chain = (TChain*)GetInputData(0);\r
1059 if(!chain) { \r
1060 Printf("ERROR: Could not receive input chain");\r
1061 return;\r
1062 }\r
1063 TObjString fileName(chain->GetCurrentFile()->GetName());\r
1064\r
1065 // trigger\r
1066 if(evtCuts->IsTriggerRequired()) \r
1067 {\r
1068 // always MB\r
1069 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
1070\r
1071 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
1072 if(!physicsSelection) return;\r
1073 //SetPhysicsTriggerSelection(physicsSelection);\r
1074\r
1075 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
1076 // set trigger (V0AND)\r
1077 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
1078 if(!triggerAnalysis) return;\r
1079 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
1080 }\r
1081 }\r
1082\r
1083 // centrality determination\r
1084 Float_t centralityF = -1;\r
1085 AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
1086 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
1087\r
1088\r
1089 // get reconstructed vertex \r
1090 //const AliESDVertex* vtxESD = 0; \r
1091 const AliESDVertex* vtxESD = 0; \r
1092 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
1093 vtxESD = esdEvent->GetPrimaryVertexTPC();\r
1094 }\r
1095 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
1096 vtxESD = esdEvent->GetPrimaryVertexTracks();\r
1097 }\r
1098 else {\r
1099 return;\r
1100 }\r
1101\r
1102 if(!vtxESD) return;\r
1103\r
1104 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
1105 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
1106 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
1107\r
1108 // check event cuts\r
1109 if(isEventOK && isEventTriggered) {\r
1110 //\r
1111 // Dump the pt downscaled V0 into the tree\r
1112 // \r
91326969 1113 Int_t ntracks = esdEvent->GetNumberOfTracks();\r
1114 Int_t nV0s = esdEvent->GetNumberOfV0s();\r
1115 Int_t run = esdEvent->GetRunNumber();\r
1116 Int_t time= esdEvent->GetTimeStamp();\r
1117 Int_t evNr=esdEvent->GetEventNumberInFile();\r
1118 \r
99496190 1119\r
1120\r
1121\r
91326969 1122 for (Int_t iv0=0; iv0<nV0s; iv0++){\r
1123 AliESDv0 * v0 = esdEvent->GetV0(iv0);\r
1124 if (!v0) continue;\r
1125 AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));\r
1126 AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));\r
1127 if (!track0) continue;\r
1128 if (!track1) continue;\r
1129 if (track0->GetSign()<0) {\r
1130 track1 = esdEvent->GetTrack(v0->GetIndex(0));\r
1131 track0 = esdEvent->GetTrack(v0->GetIndex(1));\r
1132 }\r
1133 //\r
1134 Bool_t isDownscaled = IsV0Downscaled(v0);\r
1135 if (isDownscaled) continue;\r
1136 AliKFParticle kfparticle; //\r
1137 Int_t type=GetKFParticle(v0,esdEvent,kfparticle);\r
1138 if (type==0) continue; \r
1139\r
1140 if(!fTreeSRedirector) return;\r
1141 (*fTreeSRedirector)<<"V0s"<<\r
1142 "isDownscaled="<<isDownscaled<<\r
99496190 1143 "fileName.="<<&fileName<<\r
1144 "runNumber="<<run<<\r
1145 "evtTimeStamp="<<time<<\r
1146 "evtNumberInFile="<<evNr<<\r
91326969 1147 "type="<<type<<\r
1148 "ntracks="<<ntracks<<\r
1149 "v0.="<<v0<<\r
1150 "kf.="<<&kfparticle<<\r
1151 "track0.="<<track0<<\r
1152 "track1.="<<track1<<\r
1153 "centralityF="<<centralityF<<\r
1154 "\n";\r
1155 }\r
1156 }\r
298fa346 1157 PostData(1, fOutput);\r
91326969 1158}\r
1159\r
99496190 1160\r
91326969 1161//_____________________________________________________________________________\r
1162Int_t AlidNdPtTrackDumpTask::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)\r
1163{\r
1164 //\r
1165 // Create KF particle in case the V0 fullfill selection criteria\r
1166 //\r
1167 // Selection criteria\r
1168 // 0. algorithm cut\r
1169 // 1. track cut\r
1170 // 3. chi2 cut\r
1171 // 4. rough mass cut\r
1172 // 5. Normalized pointing angle cut\r
1173 //\r
1174 const Double_t cutMass=0.2;\r
1175 const Double_t kSigmaDCACut=3;\r
1176 //\r
1177 // 0.) algo cut - accept only on the fly\r
1178 //\r
1179 if (v0->GetOnFlyStatus() ==kFALSE) return 0; \r
1180 //\r
1181 // 1.) track cut\r
1182 // \r
1183 AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));\r
1184 AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));\r
1185 /*\r
1186 TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";\r
1187 TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";\r
1188 TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";\r
1189 */ \r
1190 if (TMath::Abs(track0->GetTgl())>1) return 0;\r
1191 if (TMath::Abs(track1->GetTgl())>1) return 0;\r
1192 if ((track0->GetTPCClusterInfo(2,1))<100) return 0;\r
1193 if ((track1->GetTPCClusterInfo(2,1))<100) return 0;\r
1194 //if ((track0->GetITSclusters(0))<2) return 0;\r
1195 //if ((track1->GetITSclusters(0))<2) return 0; \r
1196 Float_t pos0[2]={0}, cov0[3]={0};\r
1197 Float_t pos1[2]={0}, cov1[3]={0};\r
1198 track0->GetImpactParameters(pos0,cov0);\r
1199 track0->GetImpactParameters(pos1,cov1);\r
1200 //\r
1201 if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;\r
1202 if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;\r
1203 // \r
1204 //\r
1205 // 3.) Chi2 cut\r
1206 //\r
1207 Double_t chi2KF = v0->GetKFInfo(2,2,2);\r
1208 if (chi2KF>25) return 0;\r
1209 //\r
1210 // 4.) Rough mass cut - 0.200 GeV\r
1211 //\r
1212 static Double_t masses[2]={-1};\r
1213 if (masses[0]<0){\r
1214 masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();\r
1215 masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();\r
1216 }\r
1217 Double_t mass00= v0->GetEffMass(0,0);\r
1218 Double_t mass22= v0->GetEffMass(2,2);\r
1219 Double_t mass42= v0->GetEffMass(4,2);\r
1220 Double_t mass24= v0->GetEffMass(2,4);\r
1221 Bool_t massOK=kFALSE;\r
1222 Int_t type=0;\r
1223 Int_t ptype=0;\r
1224 Double_t dmass=1;\r
1225 Int_t p1=0, p2=0;\r
1226 if (TMath::Abs(mass00-0)<cutMass) {\r
1227 massOK=kTRUE; type+=1; \r
1228 if (TMath::Abs(mass00-0)<dmass) {\r
1229 ptype=1;\r
1230 dmass=TMath::Abs(mass00-0); \r
1231 p1=0; p2=0;\r
1232 } \r
1233 }\r
1234 if (TMath::Abs(mass24-masses[1])<cutMass) {\r
1235 massOK=kTRUE; type+=2; \r
1236 if (TMath::Abs(mass24-masses[1])<dmass){\r
1237 dmass = TMath::Abs(mass24-masses[1]);\r
1238 ptype=2;\r
1239 p1=2; p2=4;\r
1240 }\r
1241 }\r
1242 if (TMath::Abs(mass42-masses[1])<cutMass) {\r
1243 massOK=kTRUE; type+=4;\r
1244 if (TMath::Abs(mass42-masses[1])<dmass){\r
1245 dmass = TMath::Abs(mass42-masses[1]);\r
1246 ptype=4;\r
1247 p1=4; p2=2;\r
1248 }\r
1249 }\r
1250 if (TMath::Abs(mass22-masses[0])<cutMass) {\r
1251 massOK=kTRUE; type+=8;\r
1252 if (TMath::Abs(mass22-masses[0])<dmass){\r
1253 dmass = TMath::Abs(mass22-masses[0]);\r
1254 ptype=8;\r
1255 p1=2; p2=2;\r
1256 }\r
1257 }\r
1258 if (type==0) return 0;\r
1259 //\r
1260 const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};\r
1261 const AliExternalTrackParam *paramP = v0->GetParamP();\r
1262 const AliExternalTrackParam *paramN = v0->GetParamN();\r
1263 if (paramP->GetSign()<0){\r
1264 paramP=v0->GetParamP();\r
1265 paramN=v0->GetParamN();\r
1266 }\r
1267 //Double_t *pparam1 = (Double_t*)paramP->GetParameter();\r
1268 //Double_t *pparam2 = (Double_t*)paramN->GetParameter();\r
1269 //\r
1270 AliKFParticle kfp1( *paramP, spdg[p1] );\r
1271 AliKFParticle kfp2( *paramN, -1 *spdg[p2] );\r
1272 AliKFParticle V0KF;\r
1273 (V0KF)+=kfp1;\r
1274 (V0KF)+=kfp2;\r
1275 kfparticle=V0KF;\r
1276 //\r
1277 // Pointing angle\r
1278 //\r
1279 Double_t errPhi = V0KF.GetErrPhi();\r
1280 Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());\r
1281 if (pointAngle/errPhi>10) return 0; \r
1282 //\r
1283 return ptype; \r
1284}\r
1285\r
1286//_____________________________________________________________________________\r
1287Bool_t AlidNdPtTrackDumpTask::IsV0Downscaled(AliESDv0 *const v0)\r
1288{\r
1289 //\r
1290 // Downscale randomly low pt V0\r
1291 //\r
1292 //return kFALSE;\r
1293 Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());\r
1294 Double_t scalempt= TMath::Min(maxPt,10.);\r
99496190 1295 Double_t downscaleF = gRandom->Rndm();\r
1296 downscaleF *= fLowPtV0DownscaligF;\r
1297 \r
1298 //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
1299 if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;\r
91326969 1300 return kFALSE;\r
99496190 1301\r
91326969 1302 /*\r
91326969 1303 TH1F his1("his1","his1",100,0,10);\r
1304 TH1F his2("his2","his2",100,0,10);\r
1305 {for (Int_t i=0; i<10000; i++){\r
1306 Double_t rnd=gRandom->Exp(1);\r
1307 Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();\r
1308 his1->Fill(rnd); \r
1309 if (!isDownscaled) his2->Fill(rnd); \r
1310 }}\r
1311\r
1312 */\r
1313\r
1314}\r
1315\r
1316\r
1317\r
1318\r
1319\r
1320\r
1321\r
1322\r
a26e43aa 1323\r
1324//_____________________________________________________________________________\r
1325Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
1326{\r
1327 // Constrain TPC inner params constrained\r
1328 //\r
1329 if(!tpcInnerC) return kFALSE; \r
1330 if(!vtx) return kFALSE;\r
1331\r
1332 Double_t dz[2],cov[3];\r
1333 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1334 //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1335 //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1336 if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1337\r
1338\r
1339 Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1340 Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
1341 Double_t c[3]={covar[2],0.,covar[5]};\r
1342 Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
1343 if (chi2C>kVeryBig) return kFALSE; \r
1344\r
1345 if(!tpcInnerC->Update(p,c)) return kFALSE;\r
1346\r
1347 return kTRUE;\r
1348}\r
1349\r
1350//_____________________________________________________________________________\r
1351Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
1352{\r
1353 // Constrain TPC inner params constrained\r
1354 //\r
1355 if(!trackInnerC) return kFALSE; \r
1356 if(!vtx) return kFALSE;\r
1357\r
1358 const Double_t kRadius = 2.8; \r
1359 const Double_t kMaxStep = 1.0; \r
1360\r
1361 Double_t dz[2],cov[3];\r
1362\r
1363 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1364 //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1365 //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1366\r
1367 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
1368 if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1369\r
1370 //\r
1371 Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1372 Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
1373 Double_t c[3]={covar[2],0.,covar[5]};\r
1374 Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
1375 if (chi2C>kVeryBig) return kFALSE; \r
1376\r
1377 if(!trackInnerC->Update(p,c)) return kFALSE;\r
1378\r
1379 return kTRUE;\r
1380}\r
1381\r
1382\r
1383//_____________________________________________________________________________\r
1384TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
1385{\r
1386 if(!particle) return NULL;\r
1387 if(!stack) return NULL;\r
1388\r
1389 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
1390 TParticle* mother = NULL; \r
1391 mother = stack->Particle(motherLabel); \r
1392\r
1393return mother;\r
1394}\r
1395\r
1396//_____________________________________________________________________________\r
1397Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
1398{\r
1399 Bool_t isFromConversion = kFALSE;\r
1400\r
1401 if(stack) {\r
1402 TParticle* particle = stack->Particle(label);\r
1403\r
1404 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1405 {\r
1406 Int_t mech = particle->GetUniqueID(); // production mechanism \r
1407 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1408\r
1409 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
1410 TParticle* mother = stack->Particle(motherLabel); \r
1411 if(mother) {\r
1412 Int_t motherPdg = mother->GetPdgCode();\r
1413\r
1414 if(!isPrim && mech==5 && motherPdg==kGamma) { \r
1415 isFromConversion=kTRUE; \r
1416 }\r
1417 }\r
1418 } \r
1419 } \r
1420\r
1421 return isFromConversion;\r
1422}\r
1423\r
1424//_____________________________________________________________________________\r
1425Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
1426{\r
1427 Bool_t isFromMaterial = kFALSE;\r
1428\r
1429 if(stack) {\r
1430 TParticle* particle = stack->Particle(label);\r
1431\r
1432 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1433 {\r
1434 Int_t mech = particle->GetUniqueID(); // production mechanism \r
1435 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1436\r
1437 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
1438 TParticle* mother = stack->Particle(motherLabel); \r
1439 if(mother) {\r
1440 if(!isPrim && mech==13) { \r
1441 isFromMaterial=kTRUE; \r
1442 }\r
1443 }\r
1444 } \r
1445 } \r
1446\r
1447 return isFromMaterial;\r
1448}\r
1449\r
1450//_____________________________________________________________________________\r
1451Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
1452{\r
1453 Bool_t isFromStrangeness = kFALSE;\r
1454\r
1455 if(stack) {\r
1456 TParticle* particle = stack->Particle(label);\r
1457\r
1458 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1459 {\r
1460 Int_t mech = particle->GetUniqueID(); // production mechanism \r
1461 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1462\r
1463 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
1464 TParticle* mother = stack->Particle(motherLabel); \r
1465 if(mother) {\r
1466 Int_t motherPdg = mother->GetPdgCode();\r
1467\r
1468 // K+-, lambda, antilambda, K0s decays\r
1469 if(!isPrim && mech==4 && \r
1470 (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
1471 {\r
1472 isFromStrangeness = kTRUE;\r
1473 } \r
1474 }\r
1475 } \r
1476 } \r
1477\r
1478 return isFromStrangeness;\r
1479}\r
1480\r
1481\r
1482//_____________________________________________________________________________\r
1483void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
1484{\r
1485 //\r
1486 // Called one at the end \r
1487 // locally on working node\r
1488 //\r
1489\r
99496190 1490 // must be deleted to store trees\r
1491 if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0;\r
1492\r
1493 // open temporary file and copy trees to the ouptut container\r
1494\r
298fa346 1495 TChain* chain = 0;\r
1496 TTree* tree1 = 0;\r
1497 TTree* tree2 = 0;\r
1498 TTree* tree3 = 0;\r
1499 TTree* tree4 = 0;\r
99496190 1500 //\r
1501 chain = new TChain("dNdPtTree");\r
298fa346 1502 if(chain) { \r
1503 chain->Add("jotwinow_Temp_Trees.root");\r
1504 tree1 = chain->CopyTree("1");\r
1505 delete chain; chain=0; \r
1506 }\r
99496190 1507 if(tree1) tree1->Print();\r
298fa346 1508\r
99496190 1509 //\r
1510 chain = new TChain("V0s");\r
298fa346 1511 if(chain) { \r
1512 chain->Add("jotwinow_Temp_Trees.root");\r
1513 tree2 = chain->CopyTree("1");\r
1514 delete chain; chain=0; \r
1515 }\r
99496190 1516 if(tree2) tree2->Print();\r
298fa346 1517\r
99496190 1518 //\r
1519 chain = new TChain("dEdx");\r
298fa346 1520 if(chain) { \r
1521 chain->Add("jotwinow_Temp_Trees.root");\r
1522 tree3 = chain->CopyTree("1");\r
1523 delete chain; chain=0; \r
1524 }\r
99496190 1525 if(tree3) tree3->Print();\r
298fa346 1526\r
99496190 1527 //\r
1528 chain = new TChain("Laser");\r
298fa346 1529 if(chain) { \r
1530 chain->Add("jotwinow_Temp_Trees.root");\r
1531 tree4 = chain->CopyTree("1");\r
1532 delete chain; chain=0; \r
1533 }\r
99496190 1534 if(tree4) tree4->Print();\r
1535\r
298fa346 1536\r
1537 OpenFile(1);\r
99496190 1538\r
1539 if(tree1) fOutput->Add(tree1);\r
1540 if(tree2) fOutput->Add(tree2);\r
1541 if(tree3) fOutput->Add(tree3);\r
1542 if(tree4) fOutput->Add(tree4);\r
1543 \r
a26e43aa 1544 // Post output data.\r
298fa346 1545 PostData(1, fOutput);\r
a26e43aa 1546}\r
1547\r
1548//_____________________________________________________________________________\r
1549void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
1550{\r
1551 // Called one at the end \r
99496190 1552 /*\r
1553 fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));\r
a26e43aa 1554 if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
a26e43aa 1555 TChain* chain = new TChain("dNdPtTree");\r
1556 if(!chain) return;\r
91326969 1557 chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");\r
a26e43aa 1558 TTree *tree = chain->CopyTree("1");\r
1559 if (chain) { delete chain; chain=0; }\r
1560 if(!tree) return;\r
1561 tree->Print();\r
1562 fOutputSummary = tree;\r
1563\r
1564 if (!fOutputSummary) {\r
99496190 1565 Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable %p \n", GetOutputData(1));\r
a26e43aa 1566 return;\r
1567 }\r
99496190 1568 */\r
1569\r
298fa346 1570 PostData(1, fOutput);\r
a26e43aa 1571\r
a26e43aa 1572}\r