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