- classes for pT spectra charged hadrons analysis 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
19\r
20#include "TChain.h"\r
21#include "TTreeStream.h"\r
22#include "TTree.h"\r
23#include "TH1F.h"\r
24#include "TCanvas.h"\r
25#include "TList.h"\r
26#include "TFile.h"\r
27#include "TMatrixD.h"\r
08f8415e 28#include "TRandom.h"\r
a26e43aa 29\r
30#include "AliHeader.h" \r
31#include "AliGenEventHeader.h" \r
32#include "AliInputEventHandler.h" \r
33#include "AliStack.h" \r
34#include "AliTrackReference.h" \r
35\r
36#include "AliPhysicsSelection.h"\r
37#include "AliAnalysisTask.h"\r
38#include "AliAnalysisManager.h"\r
39#include "AliESDEvent.h"\r
40#include "AliESDfriend.h"\r
41#include "AliMCEvent.h"\r
42#include "AliESDInputHandler.h"\r
43#include "AliESDVertex.h"\r
44#include "AliTracker.h"\r
45#include "AliGeomManager.h"\r
46\r
47#include "AliCentrality.h"\r
48#include "AliESDVZERO.h"\r
49#include "AliMultiplicity.h"\r
50\r
51#include "AliESDtrackCuts.h"\r
52#include "AliMCEventHandler.h"\r
bdd49ee6 53#include "AlidNdPt.h"\r
54#include "AlidNdPtEventCuts.h"\r
55#include "AlidNdPtAcceptanceCuts.h"\r
a26e43aa 56\r
bdd49ee6 57#include "AlidNdPtTrackDumpTask.h"\r
a26e43aa 58\r
59using namespace std;\r
60\r
61ClassImp(AlidNdPtTrackDumpTask)\r
62\r
63//_____________________________________________________________________________\r
64AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name) \r
65 : AliAnalysisTaskSE(name)\r
66 , fESD(0)\r
67 , fMC(0)\r
68 , fESDfriend(0)\r
69 , fOutput(0)\r
70 , fPitList(0)\r
71 , fUseMCInfo(kFALSE)\r
72 , fdNdPtEventCuts(0)\r
73 , fdNdPtAcceptanceCuts(0)\r
74 , fdNdPtRecAcceptanceCuts(0)\r
75 , fEsdTrackCuts(0)\r
76 , fTrigger(AliTriggerAnalysis::kMB1) \r
77 , fAnalysisMode(AlidNdPtHelper::kTPC) \r
78 , fOutputSummary(0)\r
79 , fTreeSRedirector(0)\r
80 , fCentralityEstimator(0)\r
81{\r
82 // Constructor\r
83\r
84 // Define input and output slots here\r
85 DefineOutput(0, TTree::Class());\r
86 //DefineOutput(1, TList::Class());\r
87\r
88}\r
89\r
90//_____________________________________________________________________________\r
91AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()\r
92{\r
93 if(fOutput) delete fOutput; fOutput =0; \r
d55b2a83 94 //if(fOutputSummary) delete fOutputSummary; fOutputSummary =0; \r
a26e43aa 95 if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector =0; \r
96\r
97 if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL; \r
98 if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;\r
99 if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL; \r
100 if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;\r
101}\r
102\r
103//____________________________________________________________________________\r
104Bool_t AlidNdPtTrackDumpTask::Notify()\r
105{\r
106 static Int_t count = 0;\r
107 count++;\r
108 TTree *chain = (TChain*)GetInputData(0);\r
109 if(chain)\r
110 {\r
111 Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());\r
112 }\r
113\r
114return kTRUE;\r
115}\r
116\r
117//_____________________________________________________________________________\r
118void AlidNdPtTrackDumpTask::UserCreateOutputObjects()\r
119{\r
120 // Create histograms\r
121 // Called once\r
122 fOutput = new TList;\r
123 fOutput->SetOwner();\r
124\r
125 //\r
126 // create output tree\r
127 //\r
128 fTreeSRedirector = new TTreeSRedirector("dNdPtOutliersAnalysisPbPb.root");\r
129\r
130 PostData(0, fOutputSummary);\r
131 //PostData(1, fOutput);\r
132}\r
133\r
134//_____________________________________________________________________________\r
135void AlidNdPtTrackDumpTask::UserExec(Option_t *) \r
136{\r
137 //\r
138 // Called for each event\r
139 //\r
140\r
141 // ESD event\r
142 fESD = (AliESDEvent*) (InputEvent());\r
143 if (!fESD) {\r
144 Printf("ERROR: ESD event not available");\r
145 return;\r
146 }\r
147\r
148 // MC event\r
149 if(fUseMCInfo) {\r
150 fMC = MCEvent();\r
151 if (!fMC) {\r
152 Printf("ERROR: MC event not available");\r
153 return;\r
154 }\r
155 }\r
156\r
157 fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
158 if(!fESDfriend) {\r
159 Printf("ERROR: ESD friends not available");\r
160 }\r
161\r
162 //\r
163 Process(fESD,fMC,fESDfriend);\r
164\r
165 // Post output data.\r
166 PostData(0, fOutputSummary);\r
167 //PostData(1, fOutput);\r
168}\r
169\r
170//_____________________________________________________________________________\r
171void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
172{\r
173 //\r
174 // Process real and/or simulated events\r
175 //\r
176 if(!esdEvent) {\r
177 AliDebug(AliLog::kError, "esdEvent not available");\r
178 return;\r
179 }\r
180\r
181 // get selection cuts\r
182 AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
183 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
184 AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
185\r
186 if(!evtCuts || !accCuts || !esdTrackCuts) {\r
187 AliDebug(AliLog::kError, "cuts not available");\r
188 return;\r
189 }\r
190\r
191 // trigger selection\r
192 Bool_t isEventTriggered = kTRUE;\r
193 AliPhysicsSelection *physicsSelection = NULL;\r
194 AliTriggerAnalysis* triggerAnalysis = NULL;\r
195\r
196 // \r
197 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
198 if (!inputHandler)\r
199 {\r
200 Printf("ERROR: Could not receive input handler");\r
201 return;\r
202 }\r
bdd49ee6 203 \r
204 // get file name\r
205 TTree *chain = (TChain*)GetInputData(0);\r
206 if(!chain) { \r
207 Printf("ERROR: Could not receive input chain");\r
208 return;\r
209 }\r
210 TObjString fileName(chain->GetCurrentFile()->GetName());\r
a26e43aa 211\r
bdd49ee6 212 // trigger\r
a26e43aa 213 if(evtCuts->IsTriggerRequired()) \r
214 {\r
215 // always MB\r
216 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
217\r
218 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
219 if(!physicsSelection) return;\r
220 //SetPhysicsTriggerSelection(physicsSelection);\r
221\r
222 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
223 // set trigger (V0AND)\r
224 triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
225 if(!triggerAnalysis) return;\r
226 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
227 }\r
228 }\r
229\r
230 // centrality determination\r
231 Float_t centralityF = -1;\r
232 AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
233 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
234\r
235 // use MC information\r
236 AliHeader* header = 0;\r
237 AliGenEventHeader* genHeader = 0;\r
238 AliStack* stack = 0;\r
239 TArrayF vtxMC(3);\r
240\r
241 Int_t multMCTrueTracks = 0;\r
242 if(IsUseMCInfo())\r
243 {\r
244 //\r
245 if(!mcEvent) {\r
246 AliDebug(AliLog::kError, "mcEvent not available");\r
247 return;\r
248 }\r
249 // get MC event header\r
250 header = mcEvent->Header();\r
251 if (!header) {\r
252 AliDebug(AliLog::kError, "Header not available");\r
253 return;\r
254 }\r
255 // MC particle stack\r
256 stack = mcEvent->Stack();\r
257 if (!stack) {\r
258 AliDebug(AliLog::kError, "Stack not available");\r
259 return;\r
260 }\r
261\r
262 // get MC vertex\r
263 genHeader = header->GenEventHeader();\r
264 if (!genHeader) {\r
265 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
266 return;\r
267 }\r
268 genHeader->PrimaryVertex(vtxMC);\r
269\r
270 // multipliticy of all MC primary tracks\r
271 // in Zv, pt and eta ranges)\r
272 multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
273\r
274 } // end bUseMC\r
275\r
276 // get reconstructed vertex \r
277 //const AliESDVertex* vtxESD = 0; \r
278 const AliESDVertex* vtxESD = 0; \r
279 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
280 vtxESD = esdEvent->GetPrimaryVertexTPC();\r
281 }\r
282 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
283 vtxESD = esdEvent->GetPrimaryVertexTracks();\r
284 }\r
285 else {\r
286 return;\r
287 }\r
288\r
289 if(!vtxESD) return;\r
290\r
291 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
292 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
293 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
294\r
295 // check event cuts\r
296 if(isEventOK && isEventTriggered)\r
297 {\r
08f8415e 298 TRandom random;\r
a26e43aa 299\r
300 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
301 {\r
302 AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
303 if(!track) continue;\r
304 if(track->Charge()==0) continue;\r
305 if(!esdTrackCuts->AcceptTrack(track)) continue;\r
306 if(!accCuts->AcceptTrack(track)) continue;\r
307\r
08f8415e 308 // downscale low-pT tracks\r
309 if(TMath::Exp(2*track->Pt())<1000*random.Rndm()) continue;\r
310\r
a26e43aa 311 // Dump to the tree \r
312 // vertex\r
313 // TPC constrained tracks\r
314 // InnerParams constrained tracks\r
315 // TPC-ITS tracks\r
316 // ITSout-InnerParams tracks\r
317 // chi2 distance between TPC constrained and TPC-ITS tracks\r
318 // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
319 // chi2 distance between ITSout and InnerParams tracks\r
320 // MC information\r
321 \r
322 Double_t x[3]; track->GetXYZ(x);\r
323 Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
324 Bool_t isOK = kFALSE;\r
325\r
326 //\r
327 // Constrain TPC-only tracks (TPCinner) to vertex\r
328 //\r
329 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
330 if (!tpcInner) continue;\r
331 // transform to the track reference frame \r
332 isOK = tpcInner->Rotate(track->GetAlpha());\r
333 isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
334 if(!isOK) continue;\r
335\r
336 // clone TPCinner has to be deleted\r
458d14c4 337 AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
a26e43aa 338 if (!tpcInnerC) continue;\r
339 \r
340 // constrain TPCinner \r
341 //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());\r
342 isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
343\r
344 // transform to the track reference frame \r
345 isOK = tpcInnerC->Rotate(track->GetAlpha());\r
346 isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
347\r
348 if(!isOK) {\r
349 if(tpcInnerC) delete tpcInnerC;\r
350 continue;\r
351 }\r
352\r
353\r
354 //\r
355 // Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
356 // to vertex\r
357 //\r
358 // clone track InnerParams has to be deleted\r
458d14c4 359 AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));\r
a26e43aa 360 if (!trackInnerC) continue;\r
361 \r
362 // constrain track InnerParams \r
363 isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
364\r
365 // transform to the track reference frame \r
366 isOK = trackInnerC->Rotate(track->GetAlpha());\r
367 isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
368\r
369 if(!isOK) {\r
370 if(trackInnerC) delete trackInnerC;\r
371 continue;\r
372 }\r
373 \r
374 //\r
375 // calculate chi2 between vi and vj vectors\r
376 // with covi and covj covariance matrices\r
377 // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
378 //\r
379 TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
380 TMatrixD delta(1,5), deltatrackC(1,5);\r
381 TMatrixD covarM(5,5), covarMtrackC(5,5);\r
382\r
383 for (Int_t ipar=0; ipar<5; ipar++) {\r
384 deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
385 delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
386\r
387 deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
388 deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
389\r
390 for (Int_t jpar=0; jpar<5; jpar++) {\r
391 Int_t index=track->GetIndex(ipar,jpar);\r
392 covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
393 covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
394 }\r
395 }\r
396 // chi2 distance TPC constrained and TPC+ITS\r
397 TMatrixD covarMInv = covarM.Invert();\r
398 TMatrixD mat2 = covarMInv*deltaT;\r
399 TMatrixD chi2 = delta*mat2; \r
400 //chi2.Print();\r
401\r
402 // chi2 distance TPC refitted constrained and TPC+ITS\r
403 TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
404 TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
405 TMatrixD chi2trackC = deltatrackC*mat2trackC; \r
406 //chi2trackC.Print();\r
407\r
408\r
409 //\r
410 // Propagate ITSout to TPC inner wall \r
411 // and calculate chi2 distance to track (InnerParams)\r
412 //\r
413 const Double_t kTPCRadius=85; \r
414 const Double_t kStep=3; \r
415\r
416 // clone track InnerParams has to be deleted\r
458d14c4 417 AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
a26e43aa 418 if (!trackInnerC2) continue;\r
419 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
420 {\r
421 if(trackInnerC2) delete trackInnerC2;\r
422 continue;\r
423 }\r
424\r
425 AliExternalTrackParam *outerITSc = new AliExternalTrackParam();\r
426 if(!outerITSc) continue;\r
427\r
428 TMatrixD chi2OuterITS(1,1);\r
429 chi2OuterITS(0,0) = 0;\r
430\r
431 if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
432 {\r
433 // propagate ITSout to TPC inner wall\r
434 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
435\r
436 if(friendTrack) \r
437 {\r
458d14c4 438 if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) ) \r
a26e43aa 439 {\r
440 if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
441 {\r
442 // transform ITSout to the track InnerParams reference frame \r
443 isOK = kFALSE;\r
444 isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
445 isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
446\r
447 if(!isOK) {\r
448 if(outerITSc) delete outerITSc;\r
449 if(trackInnerC2) delete trackInnerC2;\r
450 continue;\r
451 }\r
452 \r
453 //\r
454 // calculate chi2 between outerITS and innerParams\r
455 // cov without z-coordinate at the moment\r
456 //\r
457 TMatrixD deltaTouterITS(4,1);\r
458 TMatrixD deltaouterITS(1,4);\r
459 TMatrixD covarMouterITS(4,4);\r
460\r
461 Int_t kipar = 0;\r
462 Int_t kjpar = 0;\r
463 for (Int_t ipar=0; ipar<5; ipar++) {\r
464 if(ipar!=1) {\r
465 deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
466 deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
467 }\r
468\r
469 kjpar=0;\r
470 for (Int_t jpar=0; jpar<5; jpar++) {\r
471 Int_t index=outerITSc->GetIndex(ipar,jpar);\r
472 if(ipar !=1 || jpar!=1) {\r
473 covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
474 }\r
475 if(jpar!=1) kjpar++;\r
476 }\r
477 if(ipar!=1) kipar++;\r
478 }\r
479\r
480 // chi2 distance ITSout and InnerParams\r
481 TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
482 TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
483 chi2OuterITS = deltaouterITS*mat2outerITS; \r
484 //chi2OuterITS.Print();\r
485 } \r
486 }\r
487 }\r
488 }\r
489\r
490 //\r
491 // MC info\r
492 //\r
493 TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
494 TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
495 Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
496 Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
497 Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
498 Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
499 Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
500\r
501 AliTrackReference *refTPCIn = NULL;\r
502 AliTrackReference *refITS = NULL;\r
458d14c4 503 AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
a26e43aa 504 if(!trackInnerC3) continue;\r
505\r
506 if(IsUseMCInfo()) \r
507 {\r
508 if(!stack) return;\r
509\r
510 //\r
511 // global track\r
512 //\r
513 Int_t label = TMath::Abs(track->GetLabel()); \r
514 particle = stack->Particle(label);\r
515 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
516 {\r
517 particleMother = GetMother(particle,stack);\r
518 mech = particle->GetUniqueID();\r
519 isPrim = stack->IsPhysicalPrimary(label);\r
520 isFromStrangess = IsFromStrangeness(label,stack);\r
521 isFromConversion = IsFromConversion(label,stack);\r
522 isFromMaterial = IsFromMaterial(label,stack);\r
523 }\r
524\r
525 //\r
526 // TPC track\r
527 //\r
528 Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
529 particleTPC = stack->Particle(labelTPC);\r
530 if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
531 {\r
532 particleMotherTPC = GetMother(particleTPC,stack);\r
533 mechTPC = particleTPC->GetUniqueID();\r
534 isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
535 isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);\r
536 isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
537 isFromMaterialTPC = IsFromMaterial(labelTPC,stack);\r
538 }\r
539\r
540 //\r
541 // store first track reference\r
542 // for TPC track\r
543 //\r
544 TParticle *part=0;\r
545 TClonesArray *trefs=0;\r
546 Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
547\r
548 if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
549 {\r
550 Int_t nTrackRef = trefs->GetEntries();\r
551 //printf("nTrackRef %d \n",nTrackRef);\r
552\r
553 Int_t countITS = 0;\r
554 for (Int_t iref = 0; iref < nTrackRef; iref++) \r
555 {\r
556 AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
557 //printf("ref %p \n",ref);\r
558 //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
559 //printf("AliTrackReference::kTPC %d \n",AliTrackReference::kTPC);\r
560\r
561\r
562 // Ref. in the middle ITS \r
563 if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
564 {\r
565 if(!refITS && countITS==2) {\r
566 refITS = ref;\r
567 //printf("refITS %p \n",refITS);\r
568 }\r
569 countITS++;\r
570 }\r
571\r
572 // TPC\r
573 if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
574 {\r
575 if(!refTPCIn) {\r
576 refTPCIn = ref;\r
577 //printf("refTPCIn %p \n",refTPCIn);\r
578 //break;\r
579 }\r
580 }\r
581 }\r
582\r
583 // transform inner params to TrackRef\r
584 // reference frame\r
585 if(refTPCIn && trackInnerC3) \r
586 {\r
587 Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
588 isOK = kFALSE;\r
589 isOK = trackInnerC3->Rotate(kRefPhi);\r
590 isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
591\r
592 if(!isOK){\r
593 if(trackInnerC3) delete trackInnerC3;\r
594 if(refTPCIn) delete refTPCIn;\r
595 }\r
596 }\r
597 }\r
598\r
599 //\r
600 // ITS track\r
601 //\r
602 Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
603 particleITS = stack->Particle(labelITS);\r
604 if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
605 {\r
606 particleMotherITS = GetMother(particleITS,stack);\r
607 mechITS = particleITS->GetUniqueID();\r
608 isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
609 isFromStrangessITS = IsFromStrangeness(labelITS,stack);\r
610 isFromConversionITS = IsFromConversion(labelITS,stack);\r
611 isFromMaterialITS = IsFromMaterial(labelITS,stack);\r
612 }\r
613 }\r
614\r
615 //\r
616 Double_t vert[3] = {0}; \r
617 vert[0] = vtxESD->GetXv();\r
618 vert[1] = vtxESD->GetYv();\r
619 vert[2] = vtxESD->GetZv();\r
620 Int_t mult = vtxESD->GetNContributors();\r
621 Double_t bz = esdEvent->GetMagneticField();\r
622 Double_t runNumber = esdEvent->GetRunNumber();\r
08f8415e 623 Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
bdd49ee6 624 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
625\r
a26e43aa 626\r
627 //\r
628 if(!fTreeSRedirector) return;\r
629 (*fTreeSRedirector)<<"dNdPtTree"<<\r
bdd49ee6 630 "fileName.="<<&fileName<<\r
a26e43aa 631 "runNumber="<<runNumber<<\r
08f8415e 632 "evtTimeStamp="<<evtTimeStamp<<\r
bdd49ee6 633 "evtNumberInFile="<<evtNumberInFile<<\r
a26e43aa 634 "Bz="<<bz<<\r
08f8415e 635 "vertX="<<vert[0]<<\r
636 "vertY="<<vert[1]<<\r
637 "vertZ="<<vert[2]<<\r
a26e43aa 638 "mult="<<mult<<\r
639 "esdTrack.="<<track<<\r
640 "extTPCInnerC.="<<tpcInnerC<<\r
641 "extInnerParamC.="<<trackInnerC<<\r
642 "extInnerParam.="<<trackInnerC2<<\r
643 "extOuterITS.="<<outerITSc<<\r
644 "extInnerParamRef.="<<trackInnerC3<<\r
645 "refTPCIn.="<<refTPCIn<<\r
646 "refITS.="<<refITS<<\r
647 "chi2TPCInnerC="<<chi2(0,0)<<\r
648 "chi2InnerC="<<chi2trackC(0,0)<<\r
649 "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
650 "centralityF="<<centralityF<<\r
651 "particle.="<<particle<<\r
652 "particleMother.="<<particleMother<<\r
653 "mech="<<mech<<\r
654 "isPrim="<<isPrim<<\r
655 "isFromStrangess="<<isFromStrangess<<\r
656 "isFromConversion="<<isFromConversion<<\r
657 "isFromMaterial="<<isFromMaterial<<\r
658 "particleTPC.="<<particleTPC<<\r
659 "particleMotherTPC.="<<particleMotherTPC<<\r
660 "mechTPC="<<mechTPC<<\r
661 "isPrimTPC="<<isPrimTPC<<\r
662 "isFromStrangessTPC="<<isFromStrangessTPC<<\r
663 "isFromConversionTPC="<<isFromConversionTPC<<\r
664 "isFromMaterialTPC="<<isFromMaterialTPC<<\r
665 "particleITS.="<<particleITS<<\r
666 "particleMotherITS.="<<particleMotherITS<<\r
667 "mechITS="<<mechITS<<\r
668 "isPrimITS="<<isPrimITS<<\r
669 "isFromStrangessITS="<<isFromStrangessITS<<\r
670 "isFromConversionITS="<<isFromConversionITS<<\r
671 "isFromMaterialITS="<<isFromMaterialITS<<\r
672 "\n";\r
673\r
674 if(tpcInnerC) delete tpcInnerC;\r
675 if(trackInnerC) delete trackInnerC;\r
676 if(trackInnerC2) delete trackInnerC2;\r
677 if(outerITSc) delete outerITSc;\r
678\r
679 if(trackInnerC3) delete trackInnerC3;\r
680 }\r
681 }\r
682\r
683 PostData(0, fOutputSummary);\r
684 //PostData(1, fOutput);\r
685}\r
686\r
687\r
688//_____________________________________________________________________________\r
689Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
690{\r
691 // Constrain TPC inner params constrained\r
692 //\r
693 if(!tpcInnerC) return kFALSE; \r
694 if(!vtx) return kFALSE;\r
695\r
696 Double_t dz[2],cov[3];\r
697 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
698 //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
699 //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
700 if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
701\r
702\r
703 Double_t covar[6]; vtx->GetCovMatrix(covar);\r
704 Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
705 Double_t c[3]={covar[2],0.,covar[5]};\r
706 Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
707 if (chi2C>kVeryBig) return kFALSE; \r
708\r
709 if(!tpcInnerC->Update(p,c)) return kFALSE;\r
710\r
711 return kTRUE;\r
712}\r
713\r
714//_____________________________________________________________________________\r
715Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
716{\r
717 // Constrain TPC inner params constrained\r
718 //\r
719 if(!trackInnerC) return kFALSE; \r
720 if(!vtx) return kFALSE;\r
721\r
722 const Double_t kRadius = 2.8; \r
723 const Double_t kMaxStep = 1.0; \r
724\r
725 Double_t dz[2],cov[3];\r
726\r
727 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
728 //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
729 //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
730\r
731 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
732 if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
733\r
734 //\r
735 Double_t covar[6]; vtx->GetCovMatrix(covar);\r
736 Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
737 Double_t c[3]={covar[2],0.,covar[5]};\r
738 Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
739 if (chi2C>kVeryBig) return kFALSE; \r
740\r
741 if(!trackInnerC->Update(p,c)) return kFALSE;\r
742\r
743 return kTRUE;\r
744}\r
745\r
746\r
747//_____________________________________________________________________________\r
748TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
749{\r
750 if(!particle) return NULL;\r
751 if(!stack) return NULL;\r
752\r
753 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
754 TParticle* mother = NULL; \r
755 mother = stack->Particle(motherLabel); \r
756\r
757return mother;\r
758}\r
759\r
760//_____________________________________________________________________________\r
761Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
762{\r
763 Bool_t isFromConversion = kFALSE;\r
764\r
765 if(stack) {\r
766 TParticle* particle = stack->Particle(label);\r
767\r
768 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
769 {\r
770 Int_t mech = particle->GetUniqueID(); // production mechanism \r
771 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
772\r
773 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
774 TParticle* mother = stack->Particle(motherLabel); \r
775 if(mother) {\r
776 Int_t motherPdg = mother->GetPdgCode();\r
777\r
778 if(!isPrim && mech==5 && motherPdg==kGamma) { \r
779 isFromConversion=kTRUE; \r
780 }\r
781 }\r
782 } \r
783 } \r
784\r
785 return isFromConversion;\r
786}\r
787\r
788//_____________________________________________________________________________\r
789Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
790{\r
791 Bool_t isFromMaterial = kFALSE;\r
792\r
793 if(stack) {\r
794 TParticle* particle = stack->Particle(label);\r
795\r
796 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
797 {\r
798 Int_t mech = particle->GetUniqueID(); // production mechanism \r
799 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
800\r
801 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
802 TParticle* mother = stack->Particle(motherLabel); \r
803 if(mother) {\r
804 if(!isPrim && mech==13) { \r
805 isFromMaterial=kTRUE; \r
806 }\r
807 }\r
808 } \r
809 } \r
810\r
811 return isFromMaterial;\r
812}\r
813\r
814//_____________________________________________________________________________\r
815Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
816{\r
817 Bool_t isFromStrangeness = kFALSE;\r
818\r
819 if(stack) {\r
820 TParticle* particle = stack->Particle(label);\r
821\r
822 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
823 {\r
824 Int_t mech = particle->GetUniqueID(); // production mechanism \r
825 Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
826\r
827 Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
828 TParticle* mother = stack->Particle(motherLabel); \r
829 if(mother) {\r
830 Int_t motherPdg = mother->GetPdgCode();\r
831\r
832 // K+-, lambda, antilambda, K0s decays\r
833 if(!isPrim && mech==4 && \r
834 (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
835 {\r
836 isFromStrangeness = kTRUE;\r
837 } \r
838 }\r
839 } \r
840 } \r
841\r
842 return isFromStrangeness;\r
843}\r
844\r
845\r
846//_____________________________________________________________________________\r
847void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
848{\r
849 //\r
850 // Called one at the end \r
851 // locally on working node\r
852 //\r
853\r
854 // Post output data.\r
855 PostData(1, fOutput);\r
856 //PostData(0, fOutputSummary);\r
857}\r
858\r
859//_____________________________________________________________________________\r
860void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
861{\r
862 // Called one at the end \r
863 \r
864 // check output data\r
865 fOutputSummary = dynamic_cast<TTree*> (GetOutputData(0));\r
866 if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
867 if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0;\r
868\r
869 TChain* chain = new TChain("dNdPtTree");\r
870 if(!chain) return;\r
871 chain->Add("dNdPtOutliersAnalysisPbPb.root");\r
872 TTree *tree = chain->CopyTree("1");\r
873 if (chain) { delete chain; chain=0; }\r
874 if(!tree) return;\r
875 tree->Print();\r
876 fOutputSummary = tree;\r
877\r
878 if (!fOutputSummary) {\r
879 Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable GetOutputData(0)==0x0 ..." );\r
880 return;\r
881 }\r
882 \r
883\r
884\r
885 PostData(0, fOutputSummary);\r
886 //PostData(1, fOutput);\r
887}\r