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