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