1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 Class to process/filter reconstruction information from ESD, ESD friends, MC and provide them for later reprocessing
18 Filtering schema - low pt part is downscaled - to have flat pt specatra of selected topologies (tracks and V0s)
19 Downscaling schema is controlled by downscaling factors
21 1.) Filtering on Lego train
22 2.) expert QA for tracking (resolution efficnecy)
23 3.) pt reoslution studies using V0s
24 4.) dEdx calibration using V0s
25 5.) pt resolution and dEdx studies using cosmic
27 6.) Info used for later raw data OFFLINE triggering (highPt, V0, laser, cosmic, high dEdx)
29 Exported trees (with full objects and dereived variables):
30 1.) "highPt" - filtered trees with esd tracks, derived variables(propagated tracks), optional MC info +optional space points
31 2.) "V0s" - - filtered trees with selected V0s (rough KF chi2 cut), KF particle and corresponding esd tracks + optionla space points
32 3.) "Laser" - dump laser tracks with space points if exests
33 4.) "CosmicTree" - cosmic track candidate (random or triggered) + esdTracks(up/down)+ optional points
34 5.) "dEdx" - tree with high dEdx tpc tracks
40 #include <TDatabasePDG.h>
43 #include "TTreeStream.h"
49 #include "TObjArray.h"
54 #include "AliHeader.h"
55 #include "AliGenEventHeader.h"
56 #include "AliInputEventHandler.h"
58 #include "AliTrackReference.h"
59 #include "AliTrackPointArray.h"
60 #include "AliSysInfo.h"
62 #include "AliPhysicsSelection.h"
63 #include "AliAnalysisTask.h"
64 #include "AliAnalysisManager.h"
65 #include "AliESDEvent.h"
66 #include "AliESDfriend.h"
67 #include "AliMCEvent.h"
68 #include "AliESDInputHandler.h"
69 #include "AliESDVertex.h"
70 #include "AliTracker.h"
71 #include "AliVTrack.h"
72 #include "AliGeomManager.h"
74 #include "AliCentrality.h"
75 #include "AliESDVZERO.h"
76 #include "AliMultiplicity.h"
78 #include "AliESDtrackCuts.h"
79 #include "AliMCEventHandler.h"
80 #include "AliFilteredTreeEventCuts.h"
81 #include "AliFilteredTreeAcceptanceCuts.h"
83 #include "AliAnalysisTaskFilteredTree.h"
84 #include "AliKFParticle.h"
90 ClassImp(AliAnalysisTaskFilteredTree)
92 //_____________________________________________________________________________
93 AliAnalysisTaskFilteredTree::AliAnalysisTaskFilteredTree(const char *name)
94 : AliAnalysisTaskSE(name)
101 , fUseESDfriends(kFALSE)
102 , fReducePileUp(kTRUE)
104 , fFilteredTreeEventCuts(0)
105 , fFilteredTreeAcceptanceCuts(0)
106 , fFilteredTreeRecAcceptanceCuts(0)
108 , fTrigger(AliTriggerAnalysis::kMB1)
109 , fAnalysisMode(kTPCAnalysisMode)
110 , fTreeSRedirector(0)
111 , fCentralityEstimator(0)
112 , fLowPtTrackDownscaligF(0)
113 , fLowPtV0DownscaligF(0)
114 , fFriendDownscaling(-3.)
115 , fProcessAll(kFALSE)
116 , fProcessCosmics(kFALSE)
117 , fProcessITSTPCmatchOut(kFALSE) // swittch to process ITS/TPC standalone tracks
123 , fCosmicPairsTree(0)
126 , fPtResPhiPtTPCITS(0)
129 , fPtResEtaPtTPCITS(0)
131 , fPtResCentPtTPCc(0)
132 , fPtResCentPtTPCITS(0)
133 , fCurrentFileName("")
138 // Define input and output slots here
139 DefineOutput(1, TTree::Class());
140 DefineOutput(2, TTree::Class());
141 DefineOutput(3, TTree::Class());
142 DefineOutput(4, TTree::Class());
143 DefineOutput(5, TTree::Class());
144 DefineOutput(6, TTree::Class());
145 DefineOutput(7, TList::Class());
148 //_____________________________________________________________________________
149 AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree()
154 delete fFilteredTreeEventCuts;
155 delete fFilteredTreeAcceptanceCuts;
156 delete fFilteredTreeRecAcceptanceCuts;
157 delete fEsdTrackCuts;
160 //____________________________________________________________________________
161 Bool_t AliAnalysisTaskFilteredTree::Notify()
165 static Int_t count = 0;
167 TTree *chain = (TChain*)GetInputData(0);
169 Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());
171 fCurrentFileName=chain->GetCurrentFile()->GetName();
175 //_____________________________________________________________________________
176 void AliAnalysisTaskFilteredTree::UserCreateOutputObjects()
182 //get the output file to make sure the trees will be associated to it
184 fTreeSRedirector = new TTreeSRedirector();
188 fV0Tree = ((*fTreeSRedirector)<<"V0s").GetTree();
189 fHighPtTree = ((*fTreeSRedirector)<<"highPt").GetTree();
190 fdEdxTree = ((*fTreeSRedirector)<<"dEdx").GetTree();
191 fLaserTree = ((*fTreeSRedirector)<<"Laser").GetTree();
192 fMCEffTree = ((*fTreeSRedirector)<<"MCEffTree").GetTree();
193 fCosmicPairsTree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
196 fDummyTrack=new AliESDtrack();
201 Double_t minPt = 0.1;
202 Double_t maxPt = 100.;
205 Double_t logminPt = TMath::Log10(minPt);
206 Double_t logmaxPt = TMath::Log10(maxPt);
207 Double_t binwidth = (logmaxPt-logminPt)/nbinsPt;
208 Double_t *binsPt = new Double_t[nbinsPt+1];
210 for (Int_t i=1;i<=nbinsPt;i++) {
211 binsPt[i] = minPt + TMath::Power(10,logminPt+i*binwidth);
214 // 1pT resol cov matrix bins
215 Double_t min1PtRes = 0.;
216 Double_t max1PtRes = 0.3;
217 Int_t nbins1PtRes = 300;
218 Double_t bins1PtRes[301];
219 for (Int_t i=0;i<=nbins1PtRes;i++) {
220 bins1PtRes[i] = min1PtRes + i*(max1PtRes-min1PtRes)/nbins1PtRes;
224 Double_t minPhi = 0.;
225 Double_t maxPhi = 6.5;
226 Int_t nbinsPhi = 100;
227 Double_t binsPhi[101];
228 for (Int_t i=0;i<=nbinsPhi;i++) {
229 binsPhi[i] = minPhi + i*(maxPhi-minPhi)/nbinsPhi;
233 Double_t minEta = -1.;
234 Double_t maxEta = 1.;
236 Double_t binsEta[21];
237 for (Int_t i=0;i<=nbinsEta;i++) {
238 binsEta[i] = minEta + i*(maxEta-minEta)/nbinsEta;
242 Double_t minCent = 0.;
243 Double_t maxCent = 100;
244 Int_t nbinsCent = 20;
245 Double_t binsCent[101];
246 for (Int_t i=0;i<=nbinsCent;i++) {
247 binsCent[i] = minCent + i*(maxCent-minCent)/nbinsCent;
250 fPtResPhiPtTPC = new TH3D("fPtResPhiPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
251 fPtResPhiPtTPCc = new TH3D("fPtResPhiPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
252 fPtResPhiPtTPCITS = new TH3D("fPtResPhiPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
254 fPtResEtaPtTPC = new TH3D("fPtResEtaPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
255 fPtResEtaPtTPCc = new TH3D("fPtResEtaPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
256 fPtResEtaPtTPCITS = new TH3D("fPtResEtaPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
258 fPtResCentPtTPC = new TH3D("fPtResCentPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
259 fPtResCentPtTPCc = new TH3D("fPtResCentPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
260 fPtResCentPtTPCITS = new TH3D("fPtResCentPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
267 fOutput->Add(fPtResPhiPtTPC);
268 fOutput->Add(fPtResPhiPtTPCc);
269 fOutput->Add(fPtResPhiPtTPCITS);
270 fOutput->Add(fPtResEtaPtTPC);
271 fOutput->Add(fPtResEtaPtTPCc);
272 fOutput->Add(fPtResEtaPtTPCITS);
273 fOutput->Add(fPtResCentPtTPC);
274 fOutput->Add(fPtResCentPtTPCc);
275 fOutput->Add(fPtResCentPtTPCITS);
277 // post data to outputs
280 PostData(2,fHighPtTree);
281 PostData(3,fdEdxTree);
282 PostData(4,fLaserTree);
283 PostData(5,fMCEffTree);
284 PostData(6,fCosmicPairsTree);
289 //_____________________________________________________________________________
290 void AliAnalysisTaskFilteredTree::UserExec(Option_t *)
293 // Called for each event
297 fESD = (AliESDEvent*) (InputEvent());
299 Printf("ERROR: ESD event not available");
302 //if MC info available - use it.
305 //fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
306 fESDfriend = ESDfriend();
308 Printf("ERROR: ESD friends not available");
311 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
316 //if set, use the environment variables to set the downscaling factors
317 //AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF
318 //AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF
319 //AliAnalysisTaskFilteredTree_fFriendDownscaling
321 env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF");
323 fLowPtTrackDownscaligF=env.Atof();
324 AliInfo(Form("fLowPtTrackDownscaligF=%f",fLowPtTrackDownscaligF));
326 env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF");
328 fLowPtV0DownscaligF=env.Atof();
329 AliInfo(Form("fLowPtV0DownscaligF=%f",fLowPtTrackDownscaligF));
331 env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fFriendDownscaling");
333 fFriendDownscaling=env.Atof();
334 AliInfo(Form(" fFriendDownscaling=%f",fFriendDownscaling));
340 ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC
343 Process(fESD,fMC,fESDfriend); // only global and TPC tracks
346 ProcessV0(fESD,fMC,fESDfriend);
347 ProcessLaser(fESD,fMC,fESDfriend);
348 ProcessdEdx(fESD,fMC,fESDfriend);
349 if (fProcessCosmics) { ProcessCosmics(fESD,fESDfriend); }
350 if(fMC) { ProcessMCEff(fESD,fMC,fESDfriend);}
351 if (fProcessITSTPCmatchOut) ProcessITSTPCmatchOut(fESD, fESDfriend);
352 printf("processed event %d\n", Int_t(Entry()));
355 //_____________________________________________________________________________
356 void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliESDfriend* esdFriend)
359 // Find cosmic pairs (triggered or random)
362 AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
363 AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC();
364 const Double_t kMinPt=0.8;
365 const Double_t kMinPtMax=0.8;
366 const Double_t kMinNcl=50;
367 const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
368 Int_t ntracks=event->GetNumberOfTracks();
369 UInt_t specie = event->GetEventSpecie(); // skip laser events
370 if (specie==AliRecoParam::kCalib) return;
371 Int_t ntracksFriend = esdFriend ? esdFriend->GetNumberOfTracks() : 0;
374 for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
375 AliESDtrack *track0 = event->GetTrack(itrack0);
376 if (!track0) continue;
377 if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
379 if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue;
380 if (track0->Pt() < kMinPt) continue;
381 if (track0->GetTPCncls() < kMinNcl) continue;
382 if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
383 if (track0->GetKinkIndex(0)>0) continue;
384 const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
387 //track0->GetImpactParametersTPC(dcaTPC,covTPC);
388 //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
389 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
390 // const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
391 AliESDfriendTrack* friendTrack0=NULL;
392 if (esdFriend &&!esdFriend->TestSkipBit()){
393 if (itrack0<ntracksFriend){
394 friendTrack0 = esdFriend->GetTrack(itrack0);
395 } //this guy can be NULL
398 for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
399 AliESDtrack *track1 = event->GetTrack(itrack1);
400 if (!track1) continue;
401 if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
402 if (track1->GetKinkIndex(0)>0) continue;
403 if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue;
404 if (track1->Pt() < kMinPt) continue;
405 if (track1->GetTPCncls()<kMinNcl) continue;
406 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
407 if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
408 //track1->GetImpactParametersTPC(dcaTPC,covTPC);
409 // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
410 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
412 const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
415 for (Int_t ipar=0; ipar<5; ipar++){
416 if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
417 if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
419 if (!isPair) continue;
420 if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
421 //delta with correct sign
423 TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
424 TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
425 TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
427 if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
428 if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
429 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
430 if (!isPair) continue;
431 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
432 Int_t eventNumber = event->GetEventNumberInFile();
435 Int_t ntracksSPD = vertexSPD->GetNContributors();
436 Int_t ntracksTPC = vertexTPC->GetNContributors();
437 Int_t runNumber = event->GetRunNumber();
438 Int_t timeStamp = event->GetTimeStamp();
439 ULong64_t triggerMask = event->GetTriggerMask();
440 Float_t magField = event->GetMagneticField();
441 TObjString triggerClass = event->GetFiredTriggerClasses().Data();
443 // Global event id calculation using orbitID, bunchCrossingID and periodID
444 ULong64_t orbitID = (ULong64_t)event->GetOrbitNumber();
445 ULong64_t bunchCrossID = (ULong64_t)event->GetBunchCrossNumber();
446 ULong64_t periodID = (ULong64_t)event->GetPeriodNumber();
447 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
450 AliESDfriendTrack* friendTrack1=NULL;
451 if (esdFriend &&!esdFriend->TestSkipBit()){
452 if (itrack1<ntracksFriend){
453 friendTrack1 = esdFriend->GetTrack(itrack1);
454 } //this guy can be NULL
458 AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
459 AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
460 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
461 if (gRandom->Rndm()>1./fFriendDownscaling){
462 friendTrackStore0 = 0;
463 friendTrackStore1 = 0;
466 if (fFriendDownscaling<=0){
467 if (((*fTreeSRedirector)<<"CosmicPairs").GetTree()){
468 TTree * tree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
470 Double_t sizeAll=tree->GetZipBytes();
471 TBranch * br= tree->GetBranch("friendTrack0.fPoints");
472 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
473 br= tree->GetBranch("friendTrack0.fCalibContainer");
474 if (br) sizeFriend+=br->GetZipBytes();
475 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
482 if(!fFillTree) return;
483 if(!fTreeSRedirector) return;
484 (*fTreeSRedirector)<<"CosmicPairs"<<
485 "gid="<<gid<< // global id of track
486 "fileName.="<<&fCurrentFileName<< // file name
487 "runNumber="<<runNumber<< // run number
488 "evtTimeStamp="<<timeStamp<< // time stamp of event
489 "evtNumberInFile="<<eventNumber<< // event number
490 "trigger="<<triggerMask<< // trigger mask
491 "triggerClass="<<&triggerClass<< // trigger class
492 "Bz="<<magField<< // magnetic field
494 "multSPD="<<ntracksSPD<< // event ultiplicity
495 "multTPC="<<ntracksTPC<< //
496 "vertSPD.="<<vertexSPD<< // primary vertex -SPD
497 "vertTPC.="<<vertexTPC<< // primary vertex -TPC
498 "t0.="<<track0<< // first half of comsic trak
499 "t1.="<<track1<< // second half of cosmic track
500 "friendTrack0.="<<friendTrackStore0<< // friend information first track + points
501 "friendTrack1.="<<friendTrackStore1<< // frined information first track + points
508 //_____________________________________________________________________________
509 void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
512 // Select real events with high-pT tracks
515 // get selection cuts
516 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
517 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
518 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
520 if(!evtCuts || !accCuts || !esdTrackCuts) {
521 Printf("ERROR cuts not available");
526 Bool_t isEventTriggered = kTRUE;
527 AliPhysicsSelection *physicsSelection = NULL;
528 AliTriggerAnalysis* triggerAnalysis = NULL;
531 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
533 if(evtCuts->IsTriggerRequired())
536 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
538 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
539 if(!physicsSelection) return;
540 //SetPhysicsTriggerSelection(physicsSelection);
542 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
543 // set trigger (V0AND)
544 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
545 if(!triggerAnalysis) return;
546 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
550 // centrality determination
551 Float_t centralityF = -1;
552 AliCentrality *esdCentrality = esdEvent->GetCentrality();
553 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
555 // use MC information
556 AliHeader* header = 0;
557 AliGenEventHeader* genHeader = 0;
561 Int_t multMCTrueTracks = 0;
564 // get MC event header
565 header = mcEvent->Header();
567 AliDebug(AliLog::kError, "Header not available");
571 stack = mcEvent->Stack();
573 AliDebug(AliLog::kError, "Stack not available");
578 genHeader = header->GenEventHeader();
580 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
583 genHeader->PrimaryVertex(vtxMC);
585 // multipliticy of all MC primary tracks
586 // in Zv, pt and eta ranges)
587 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
590 // get reconstructed vertex
591 //const AliESDVertex* vtxESD = 0;
592 AliESDVertex* vtxESD = 0;
593 if(GetAnalysisMode() == kTPCAnalysisMode) {
594 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
596 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
597 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
603 AliESDVertex* vtxTPC = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
604 AliESDVertex* vtxSPD = (AliESDVertex*)esdEvent->GetPrimaryVertexSPD();
610 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
611 //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZ());
612 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
613 Int_t ntracks = esdEvent->GetNumberOfTracks();
616 if(isEventOK && isEventTriggered)
620 // get IR information
622 AliESDHeader *esdHeader = 0;
623 esdHeader = esdEvent->GetHeader();
624 if(!esdHeader) return;
625 //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
626 //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
628 // Use when Peter commit the changes in the header
633 //Double_t vert[3] = {0};
634 //vert[0] = vtxESD->GetX();
635 //vert[1] = vtxESD->GetY();
636 //vert[2] = vtxESD->GetZ();
637 Int_t mult = vtxESD->GetNContributors();
638 Int_t multSPD = vtxSPD->GetNContributors();
639 Int_t multTPC = vtxTPC->GetNContributors();
641 Float_t bz = esdEvent->GetMagneticField();
642 Int_t runNumber = esdEvent->GetRunNumber();
643 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
644 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
646 // Global event id calculation using orbitID, bunchCrossingID and periodID
647 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
648 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
649 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
650 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
654 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
656 AliESDtrack *track = esdEvent->GetTrack(iTrack);
658 if(track->Charge()==0) continue;
659 if(!esdTrackCuts->AcceptTrack(track)) continue;
660 if(!accCuts->AcceptTrack(track)) continue;
662 // downscale low-pT tracks
663 Double_t scalempt= TMath::Min(track->Pt(),10.);
664 Double_t downscaleF = gRandom->Rndm();
665 downscaleF *= fLowPtTrackDownscaligF;
666 if(TMath::Exp(2*scalempt)<downscaleF) continue;
667 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
669 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
670 if (!tpcInner) continue;
671 // transform to the track reference frame
672 Bool_t isOK = kFALSE;
673 isOK = tpcInner->Rotate(track->GetAlpha());
674 isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
681 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
682 if(!fFillTree) return;
683 if(!fTreeSRedirector) return;
684 (*fTreeSRedirector)<<"highPt"<<
686 "fileName.="<<&fCurrentFileName<<
687 "runNumber="<<runNumber<<
688 "evtTimeStamp="<<evtTimeStamp<<
689 "evtNumberInFile="<<evtNumberInFile<<
690 "triggerClass="<<&triggerClass<< // trigger
691 "Bz="<<bz<< // magnetic field
693 "ntracksESD="<<ntracks<< // number of tracks in the ESD
694 "IRtot="<<ir1<< // interaction record history info
696 "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
697 "multSPD="<<multSPD<< // multiplicity of tracks pointing to the SPD primary vertex
698 "multTPC="<<multTPC<< // multiplicity of tracks pointing to the TPC primary vertex
699 "esdTrack.="<<track<<
700 "centralityF="<<centralityF<<
708 //_____________________________________________________________________________
709 void AliAnalysisTaskFilteredTree::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const esdFriend)
712 // Process laser events -> dump tracks and clusters to the special tree
714 const Double_t kMinPt = 5;
715 if(!fFillTree) return;
716 if(!fTreeSRedirector) return;
717 const AliESDHeader* esdHeader = esdEvent->GetHeader();
718 if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) {
719 Int_t countLaserTracks = 0;
720 Int_t runNumber = esdEvent->GetRunNumber();
721 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
722 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
723 Float_t bz = esdEvent->GetMagneticField();
724 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
725 // Global event id calculation using orbitID, bunchCrossingID and periodID
726 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
727 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
728 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
729 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
730 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
731 AliESDtrack *track = esdEvent->GetTrack(iTrack);
733 if(track->GetTPCInnerParam()) countLaserTracks++;
734 AliESDfriendTrack* friendTrack=NULL;
735 // suppress beam background and CE random reacks
736 if (track->GetInnerParam()->Pt()<kMinPt) continue;
737 Bool_t skipTrack=gRandom->Rndm()>1/(1+TMath::Abs(fFriendDownscaling));
738 if (skipTrack) continue;
739 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
740 (*fTreeSRedirector)<<"Laser"<<
741 "gid="<<gid<< // global identifier of event
742 "fileName.="<<&fCurrentFileName<< //
743 "runNumber="<<runNumber<<
744 "evtTimeStamp="<<evtTimeStamp<<
745 "evtNumberInFile="<<evtNumberInFile<<
746 "triggerClass="<<&triggerClass<< // trigger
747 "Bz="<<bz<< // magnetic field
748 "multTPCtracks="<<countLaserTracks<< // multiplicity of tracks
749 "track.="<<track<< // track parameters
750 "friendTrack.="<<friendTrack<< // friend track information
756 //_____________________________________________________________________________
757 void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
760 // Select real events with high-pT tracks
761 // Calculate and stor in the output tree:
762 // TPC constrained tracks
763 // InnerParams constrained tracks
765 // ITSout-InnerParams tracks
766 // chi2 distance between TPC constrained and TPC-ITS tracks
767 // chi2 distance between TPC refitted constrained and TPC-ITS tracks
768 // chi2 distance between ITSout and InnerParams tracks
770 // track references at ITSin, TPCin; InnerParam at first TPC track reference,
771 // particle ID, mother ID, production mechanism ...
773 // get selection cuts
774 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
775 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
776 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
778 if(!evtCuts || !accCuts || !esdTrackCuts) {
779 AliDebug(AliLog::kError, "cuts not available");
784 Bool_t isEventTriggered = kTRUE;
785 AliPhysicsSelection *physicsSelection = NULL;
786 AliTriggerAnalysis* triggerAnalysis = NULL;
789 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
793 if(evtCuts->IsTriggerRequired())
796 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
798 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
799 if(!physicsSelection) {AliInfo("no physics selection"); return;}
800 //SetPhysicsTriggerSelection(physicsSelection);
802 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
803 // set trigger (V0AND)
804 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
805 if(!triggerAnalysis) {AliInfo("no trigger analysis");return;}
806 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
810 // centrality determination
811 Float_t centralityF = -1;
812 AliCentrality *esdCentrality = esdEvent->GetCentrality();
813 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
815 // use MC information
816 AliHeader* header = 0;
817 AliGenEventHeader* genHeader = 0;
822 Int_t multMCTrueTracks = 0;
825 // get MC event header
826 header = mcEvent->Header();
828 AliDebug(AliLog::kError, "Header not available");
832 stack = mcEvent->Stack();
834 AliDebug(AliLog::kError, "Stack not available");
837 mcStackSize=stack->GetNtrack();
840 genHeader = header->GenEventHeader();
842 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
845 genHeader->PrimaryVertex(vtxMC);
847 // multipliticy of all MC primary tracks
848 // in Zv, pt and eta ranges)
849 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
853 // get reconstructed vertex
854 //const AliESDVertex* vtxESD = 0;
855 AliESDVertex* vtxESD = 0;
856 if(GetAnalysisMode() == kTPCAnalysisMode) {
857 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
859 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
860 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
863 AliInfo("no ESD vertex");
869 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
872 // Vertex info comparison and track multiplicity
874 AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
875 AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
876 Int_t contSPD = vertexSPD->GetNContributors();
877 Int_t contTPC = vertexTPC->GetNContributors();
878 TVectorD vertexPosTPC(3), vertexPosSPD(3);
879 vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
880 vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
883 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
884 AliESDtrack *track = esdEvent->GetTrack(iTrack);
886 if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
887 if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
890 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
891 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
895 if(isEventOK && isEventTriggered)
898 // get IR information
900 AliESDHeader *esdHeader = 0;
901 esdHeader = esdEvent->GetHeader();
902 if(!esdHeader) {AliInfo("no esdHeader");return;}
903 //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
904 //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
910 Double_t vert[3] = {0};
911 vert[0] = vtxESD->GetX();
912 vert[1] = vtxESD->GetY();
913 vert[2] = vtxESD->GetZ();
914 Int_t mult = vtxESD->GetNContributors();
915 Float_t bz = esdEvent->GetMagneticField();
916 Int_t runNumber = esdEvent->GetRunNumber();
917 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
918 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
920 Int_t numberOfTracks=esdEvent->GetNumberOfTracks();
922 for (Int_t iTrack = 0; iTrack < numberOfTracks; iTrack++)
924 AliESDtrack *track = esdEvent->GetTrack(iTrack);
925 AliESDfriendTrack* friendTrack=NULL;
926 Int_t numberOfFriendTracks=0;
927 if (esdFriend) numberOfFriendTracks=esdFriend->GetNumberOfTracks();
928 if (esdFriend && iTrack<numberOfFriendTracks) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
930 if(track->Charge()==0) continue;
931 if(!esdTrackCuts->AcceptTrack(track)) continue;
932 if(!accCuts->AcceptTrack(track)) continue;
934 // downscale low-pT tracks
935 Double_t scalempt= TMath::Min(track->Pt(),10.);
936 Double_t downscaleF = gRandom->Rndm();
937 downscaleF *= fLowPtTrackDownscaligF;
938 if(TMath::Exp(2*scalempt)<downscaleF) continue;
939 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
943 // TPC constrained tracks
944 // InnerParams constrained tracks
946 // ITSout-InnerParams tracks
947 // chi2 distance between TPC constrained and TPC-ITS tracks
948 // chi2 distance between TPC refitted constrained and TPC-ITS tracks
949 // chi2 distance between ITSout and InnerParams tracks
952 Double_t x[3]; track->GetXYZ(x);
953 Double_t b[3]; AliTracker::GetBxByBz(x,b);
956 // Transform TPC inner params to track reference frame
958 Bool_t isOKtpcInner = kFALSE;
959 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
961 // transform to the track reference frame
962 isOKtpcInner = tpcInner->Rotate(track->GetAlpha());
963 isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
967 // Constrain TPC inner to vertex
968 // clone TPCinner has to be deleted
970 Bool_t isOKtpcInnerC = kFALSE;
971 AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
973 isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);
974 isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());
975 isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
979 // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex
980 // Clone track InnerParams has to be deleted
982 Bool_t isOKtrackInnerC = kFALSE;
983 AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));
985 isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);
986 isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());
987 isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
991 // calculate chi2 between vi and vj vectors
992 // with covi and covj covariance matrices
993 // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
995 TMatrixD deltaT(5,1), deltaTtrackC(5,1);
996 TMatrixD delta(1,5), deltatrackC(1,5);
997 TMatrixD covarM(5,5), covarMtrackC(5,5);
999 TMatrixD chi2trackC(1,1);
1001 if(isOKtpcInnerC && isOKtrackInnerC)
1003 for (Int_t ipar=0; ipar<5; ipar++) {
1004 deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
1005 delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
1007 deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
1008 deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
1010 for (Int_t jpar=0; jpar<5; jpar++) {
1011 Int_t index=track->GetIndex(ipar,jpar);
1012 covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
1013 covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];
1017 // chi2 distance TPC constrained and TPC+ITS
1018 TMatrixD covarMInv = covarM.Invert();
1019 TMatrixD mat2 = covarMInv*deltaT;
1023 // chi2 distance TPC refitted constrained and TPC+ITS
1024 TMatrixD covarMInvtrackC = covarMtrackC.Invert();
1025 TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;
1026 chi2trackC = deltatrackC*mat2trackC;
1027 //chi2trackC.Print();
1032 // Propagate ITSout to TPC inner wall
1033 // and calculate chi2 distance to track (InnerParams)
1035 const Double_t kTPCRadius=85;
1036 const Double_t kStep=3;
1038 // clone track InnerParams has to be deleted
1039 Bool_t isOKtrackInnerC2 = kFALSE;
1040 AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));
1042 isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);
1045 Bool_t isOKouterITSc = kFALSE;
1046 AliExternalTrackParam *outerITSc = NULL;
1047 TMatrixD chi2OuterITS(1,1);
1049 if(esdFriend && !esdFriend->TestSkipBit())
1051 // propagate ITSout to TPC inner wall
1056 if (friendTrack->GetITSOut()) outerITSc = new AliExternalTrackParam(*(friendTrack->GetITSOut()));
1059 isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);
1060 isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());
1061 isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());
1064 // calculate chi2 between outerITS and innerParams
1065 // cov without z-coordinate at the moment
1067 TMatrixD deltaTouterITS(4,1);
1068 TMatrixD deltaouterITS(1,4);
1069 TMatrixD covarMouterITS(4,4);
1071 if(isOKtrackInnerC2 && isOKouterITSc) {
1074 for (Int_t ipar=0; ipar<5; ipar++) {
1076 deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
1077 deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
1081 for (Int_t jpar=0; jpar<5; jpar++) {
1082 Int_t index=outerITSc->GetIndex(ipar,jpar);
1083 if(ipar !=1 || jpar!=1) {
1084 covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];
1086 if(jpar!=1) kjpar++;
1088 if(ipar!=1) kipar++;
1091 // chi2 distance ITSout and InnerParams
1092 TMatrixD covarMInvouterITS = covarMouterITS.Invert();
1093 TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;
1094 chi2OuterITS = deltaouterITS*mat2outerITS;
1095 //chi2OuterITS.Print();
1104 TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;
1105 TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;
1106 Int_t mech=-1, mechTPC=-1, mechITS=-1;
1107 Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;
1108 Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;
1109 Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;
1110 Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;
1112 AliTrackReference *refTPCIn = NULL;
1113 AliTrackReference *refTPCOut = NULL;
1114 AliTrackReference *refITS = NULL;
1115 AliTrackReference *refTRD = NULL;
1116 AliTrackReference *refTOF = NULL;
1117 AliTrackReference *refEMCAL = NULL;
1118 AliTrackReference *refPHOS = NULL;
1119 Int_t nrefTPC=0, nrefTRD=0, nrefTOF=0, nrefITS=0, nrefEMCAL=0, nrefPHOS=0;
1121 Bool_t isOKtrackInnerC3 = kFALSE;
1122 AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
1123 if(mcEvent && stack)
1125 do //artificial loop (once) to make the continue statements jump out of the MC part
1127 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1131 Int_t label = TMath::Abs(track->GetLabel());
1132 if (label >= mcStackSize) continue;
1133 particle = stack->Particle(label);
1134 if (!particle) continue;
1135 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)
1137 particleMother = GetMother(particle,stack);
1138 mech = particle->GetUniqueID();
1139 isPrim = stack->IsPhysicalPrimary(label);
1140 isFromStrangess = IsFromStrangeness(label,stack);
1141 isFromConversion = IsFromConversion(label,stack);
1142 isFromMaterial = IsFromMaterial(label,stack);
1148 Int_t labelTPC = TMath::Abs(track->GetTPCLabel());
1149 if (labelTPC >= mcStackSize) continue;
1150 particleTPC = stack->Particle(labelTPC);
1151 if (!particleTPC) continue;
1152 if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)
1154 particleMotherTPC = GetMother(particleTPC,stack);
1155 mechTPC = particleTPC->GetUniqueID();
1156 isPrimTPC = stack->IsPhysicalPrimary(labelTPC);
1157 isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);
1158 isFromConversionTPC = IsFromConversion(labelTPC,stack);
1159 isFromMaterialTPC = IsFromMaterial(labelTPC,stack);
1163 // store first track reference
1167 TClonesArray *trefs=0;
1168 Int_t status = mcEvent->GetParticleAndTR(TMath::Abs(labelTPC), part, trefs);
1170 if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.)
1172 Int_t nTrackRef = trefs->GetEntries();
1173 //printf("nTrackRef %d \n",nTrackRef);
1176 for (Int_t iref = 0; iref < nTrackRef; iref++)
1178 AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);
1180 // Ref. in the middle ITS
1181 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kITS)
1184 if(!refITS && countITS==2) {
1186 //printf("refITS %p \n",refITS);
1192 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC)
1198 //printf("refTPCIn %p \n",refTPCIn);
1203 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD)
1211 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF)
1219 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL)
1227 // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS)
1236 // transform inner params to TrackRef
1238 if(refTPCIn && trackInnerC3)
1240 Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());
1241 isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);
1242 isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);
1249 Int_t labelITS = TMath::Abs(track->GetITSLabel());
1250 if (labelITS >= mcStackSize) continue;
1251 particleITS = stack->Particle(labelITS);
1252 if (!particleITS) continue;
1253 if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)
1255 particleMotherITS = GetMother(particleITS,stack);
1256 mechITS = particleITS->GetUniqueID();
1257 isPrimITS = stack->IsPhysicalPrimary(labelITS);
1258 isFromStrangessITS = IsFromStrangeness(labelITS,stack);
1259 isFromConversionITS = IsFromConversion(labelITS,stack);
1260 isFromMaterialITS = IsFromMaterial(labelITS,stack);
1267 Bool_t dumpToTree=kFALSE;
1269 if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE;
1270 //if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1271 if(isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1272 if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE;
1273 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1276 // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks
1277 // Done only in case no MC info
1280 track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]);
1281 Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10;
1282 Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin));
1283 Bool_t keepPileUp=gRandom->Rndm()<0.05;
1284 if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){
1289 //init dummy objects
1290 static AliESDVertex dummyvtxESD;
1293 // dummyvtxESD=new AliESDVertex();
1294 // //dummyvtxESD->SetNContributors(1);
1295 // //UShort_t pindices[1]; pindices[0]=0;
1296 // //dummyvtxESD->SetIndices(1,pindices);
1297 // //dummyvtxESD->SetNContributors(0);
1299 static AliExternalTrackParam dummyexternaltrackparam;
1300 //if (!dummyexternaltrackparam) dummyexternaltrackparam=new AliExternalTrackParam();
1301 static AliTrackReference dummytrackreference;
1302 //if (!dummytrackreference) dummytrackreference=new AliTrackReference();
1303 static TParticle dummyparticle;
1304 //if (!dummyparticle) dummyparticle=new TParticle();
1306 //assign the dummy objects if needed
1307 if (!track) {track=fDummyTrack;}
1308 AliESDfriendTrack *friendTrackStore=friendTrack; // store friend track for later processing
1309 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
1310 friendTrackStore = (gRandom->Rndm()<1./fFriendDownscaling)? friendTrack:0;
1312 if (fFriendDownscaling<=0){
1313 if (((*fTreeSRedirector)<<"highPt").GetTree()){
1314 TTree * tree = ((*fTreeSRedirector)<<"highPt").GetTree();
1316 Double_t sizeAll=tree->GetZipBytes();
1317 TBranch * br= tree->GetBranch("friendTrack.fPoints");
1318 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
1319 br= tree->GetBranch("friendTrack.fCalibContainer");
1320 if (br) sizeFriend+=br->GetZipBytes();
1321 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) friendTrackStore=0;
1327 // if (!friendTrackStore && fFriendDownscaling<=1) {friendTrack=fDummyFriendTrack;}
1328 if (!vtxESD) {vtxESD=&dummyvtxESD;}
1331 if (!refTPCIn) {refTPCIn=&dummytrackreference;}
1332 if (!refITS) {refITS=&dummytrackreference;}
1333 if (!particle) {particle=&dummyparticle;}
1334 if (!particleMother) {particleMother=&dummyparticle;}
1335 if (!particleTPC) {particleTPC=&dummyparticle;}
1336 if (!particleMotherTPC) {particleMotherTPC=&dummyparticle;}
1337 if (!particleITS) {particleITS=&dummyparticle;}
1338 if (!particleMotherITS) {particleMotherITS=&dummyparticle;}
1341 Int_t ntracks = esdEvent->GetNumberOfTracks();
1343 // Global event id calculation using orbitID, bunchCrossingID and periodID
1344 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1345 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1346 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1347 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1350 FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0));
1352 if(fTreeSRedirector && dumpToTree && fFillTree) {
1353 (*fTreeSRedirector)<<"highPt"<<
1355 "fileName.="<<&fCurrentFileName<< // name of the chunk file (hopefully full)
1356 "runNumber="<<runNumber<< // runNumber
1357 "evtTimeStamp="<<evtTimeStamp<< // time stamp of event (in seconds)
1358 "evtNumberInFile="<<evtNumberInFile<< // event number
1359 "triggerClass="<<&triggerClass<< // trigger class as a string
1360 "Bz="<<bz<< // solenoid magnetic field in the z direction (in kGaus)
1361 "vtxESD.="<<vtxESD<< // vertexer ESD tracks (can be biased by TPC pileup tracks)
1362 "IRtot="<<ir1<< // interaction record (trigger) counters - coutner 1
1363 "IRint2="<<ir2<< // interaction record (trigger) coutners - counter 2
1364 "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
1365 "ntracks="<<ntracks<< // number of the esd tracks (to take into account the pileup in the TPC)
1366 // important variables for the pile-up studies
1367 "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
1368 "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
1369 "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
1370 "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
1371 "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
1372 "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
1374 "esdTrack.="<<track<< // esdTrack as used in the physical analysis
1375 // "friendTrack.="<<friendTrack<< // esdFriendTrack associated to the esdTrack
1376 "friendTrack.="<<friendTrackStore<< // esdFriendTrack associated to the esdTrack
1377 "extTPCInnerC.="<<tpcInnerC<< // ???
1378 "extInnerParamC.="<<trackInnerC<< // ???
1379 "extInnerParam.="<<trackInnerC2<< // ???
1380 "extOuterITS.="<<outerITSc<< // ???
1381 "extInnerParamRef.="<<trackInnerC3<< // ???
1382 "chi2TPCInnerC="<<chi2(0,0)<< // chi2 of tracks ???
1383 "chi2InnerC="<<chi2trackC(0,0)<< // chi2s of tracks TPCinner to the combined
1384 "chi2OuterITS="<<chi2OuterITS(0,0)<< // chi2s of tracks TPC at inner wall to the ITSout
1385 "centralityF="<<centralityF;
1387 static AliTrackReference refDummy;
1388 if (!refITS) refITS = &refDummy;
1389 if (!refTRD) refTRD = &refDummy;
1390 if (!refTOF) refTOF = &refDummy;
1391 if (!refEMCAL) refEMCAL = &refDummy;
1392 if (!refPHOS) refPHOS = &refDummy;
1393 (*fTreeSRedirector)<<"highPt"<<
1394 "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
1395 "nrefITS="<<nrefITS<< // number of track references in the ITS
1396 "nrefTPC="<<nrefTPC<< // number of track references in the TPC
1397 "nrefTRD="<<nrefTRD<< // number of track references in the TRD
1398 "nrefTOF="<<nrefTOF<< // number of track references in the TOF
1399 "nrefEMCAL="<<nrefEMCAL<< // number of track references in the TOF
1400 "nrefPHOS="<<nrefPHOS<< // number of track references in the TOF
1401 "refTPCIn.="<<refTPCIn<<
1402 "refTPCOut.="<<refTPCOut<<
1403 "refITS.="<<refITS<<
1404 "refTRD.="<<refTRD<<
1405 "refTOF.="<<refTOF<<
1406 "refEMCAL.="<<refEMCAL<<
1407 "refPHOS.="<<refPHOS<<
1408 "particle.="<<particle<<
1409 "particleMother.="<<particleMother<<
1412 "isFromStrangess="<<isFromStrangess<<
1413 "isFromConversion="<<isFromConversion<<
1414 "isFromMaterial="<<isFromMaterial<<
1415 "particleTPC.="<<particleTPC<<
1416 "particleMotherTPC.="<<particleMotherTPC<<
1417 "mechTPC="<<mechTPC<<
1418 "isPrimTPC="<<isPrimTPC<<
1419 "isFromStrangessTPC="<<isFromStrangessTPC<<
1420 "isFromConversionTPC="<<isFromConversionTPC<<
1421 "isFromMaterialTPC="<<isFromMaterialTPC<<
1422 "particleITS.="<<particleITS<<
1423 "particleMotherITS.="<<particleMotherITS<<
1424 "mechITS="<<mechITS<<
1425 "isPrimITS="<<isPrimITS<<
1426 "isFromStrangessITS="<<isFromStrangessITS<<
1427 "isFromConversionITS="<<isFromConversionITS<<
1428 "isFromMaterialITS="<<isFromMaterialITS;
1430 //finish writing the entry
1431 AliInfo("writing tree highPt");
1432 (*fTreeSRedirector)<<"highPt"<<"\n";
1434 AliSysInfo::AddStamp("filteringTask",iTrack,numberOfTracks,numberOfFriendTracks,(friendTrackStore)?0:1);
1437 delete trackInnerC2;
1439 delete trackInnerC3;
1445 //_____________________________________________________________________________
1446 void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1449 // Fill tree for efficiency studies MC only
1450 AliInfo("we start!");
1452 AliDebug(AliLog::kError, "mcEvent not available");
1456 // get selection cuts
1457 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1458 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1459 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1461 if(!evtCuts || !accCuts || !esdTrackCuts) {
1462 AliDebug(AliLog::kError, "cuts not available");
1466 // trigger selection
1467 Bool_t isEventTriggered = kTRUE;
1468 AliPhysicsSelection *physicsSelection = NULL;
1469 AliTriggerAnalysis* triggerAnalysis = NULL;
1472 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1476 if(evtCuts->IsTriggerRequired())
1479 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1480 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1481 if(!physicsSelection) return;
1483 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1484 // set trigger (V0AND)
1485 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1486 if(!triggerAnalysis) return;
1487 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1491 // centrality determination
1492 Float_t centralityF = -1;
1493 AliCentrality *esdCentrality = esdEvent->GetCentrality();
1494 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1496 // use MC information
1497 AliHeader* header = 0;
1498 AliGenEventHeader* genHeader = 0;
1499 AliStack* stack = 0;
1500 Int_t mcStackSize=0;
1503 Int_t multMCTrueTracks = 0;
1506 AliDebug(AliLog::kError, "mcEvent not available");
1509 // get MC event header
1510 header = mcEvent->Header();
1512 AliDebug(AliLog::kError, "Header not available");
1515 // MC particle stack
1516 stack = mcEvent->Stack();
1518 AliDebug(AliLog::kError, "Stack not available");
1521 mcStackSize=stack->GetNtrack();
1524 genHeader = header->GenEventHeader();
1526 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1529 genHeader->PrimaryVertex(vtxMC);
1531 // multipliticy of all MC primary tracks
1532 // in Zv, pt and eta ranges)
1533 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1536 // get reconstructed vertex
1537 //const AliESDVertex* vtxESD = 0;
1538 AliESDVertex* vtxESD = 0;
1539 if(GetAnalysisMode() == kTPCAnalysisMode) {
1540 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1542 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1543 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1551 // Vertex info comparison and track multiplicity
1553 AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
1554 AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
1555 Int_t contSPD = vertexSPD->GetNContributors();
1556 Int_t contTPC = vertexTPC->GetNContributors();
1557 TVectorD vertexPosTPC(3), vertexPosSPD(3);
1558 vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
1559 vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
1562 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
1563 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1564 if(!track) continue;
1565 if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
1566 if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
1569 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1570 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1571 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1573 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1576 if(isEventOK && isEventTriggered)
1578 //if(!stack) return;
1583 TParticle *particle=NULL;
1584 TParticle *particleMother=NULL;
1588 Double_t vert[3] = {0};
1589 vert[0] = vtxESD->GetX();
1590 vert[1] = vtxESD->GetY();
1591 vert[2] = vtxESD->GetZ();
1592 Int_t mult = vtxESD->GetNContributors();
1593 Double_t bz = esdEvent->GetMagneticField();
1594 Double_t runNumber = esdEvent->GetRunNumber();
1595 Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1596 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1597 // loop over MC stack
1598 for (Int_t iMc = 0; iMc < mcStackSize; ++iMc)
1600 particle = stack->Particle(iMc);
1604 // only charged particles
1605 if(!particle->GetPDG()) continue;
1606 Double_t charge = particle->GetPDG()->Charge()/3.;
1607 if (TMath::Abs(charge) < 0.001)
1610 // only primary particles
1611 Bool_t prim = stack->IsPhysicalPrimary(iMc);
1614 // downscale low-pT particles
1615 Double_t scalempt= TMath::Min(particle->Pt(),10.);
1616 Double_t downscaleF = gRandom->Rndm();
1617 downscaleF *= fLowPtTrackDownscaligF;
1618 if(TMath::Exp(2*scalempt)<downscaleF) continue;
1620 // is particle in acceptance
1621 if(!accCuts->AcceptTrack(particle)) continue;
1623 // check if particle reconstructed
1624 Bool_t isRec = kFALSE;
1625 Int_t trackIndex = -1;
1626 Int_t trackLoopIndex = -1;
1627 Int_t isESDtrackCut= 0;
1628 Int_t isAccCuts = 0;
1629 Int_t nRec = 0; // how many times reconstructed
1630 Int_t nFakes = 0; // how many times reconstructed as a fake track
1631 AliESDtrack *recTrack = NULL;
1633 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1635 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1636 if(!track) continue;
1637 if(track->Charge()==0) continue;
1639 Int_t label = TMath::Abs(track->GetLabel());
1640 if (label >= mcStackSize) continue;
1642 Bool_t isAcc=esdTrackCuts->AcceptTrack(track);
1643 if (isAcc) isESDtrackCut=1;
1644 if (accCuts->AcceptTrack(track)) isAccCuts=1;
1646 trackIndex = iTrack;
1649 if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track
1650 if (!isAcc) continue;
1651 trackLoopIndex = iTrack;
1653 recTrack = esdEvent->GetTrack(trackIndex);
1655 if(track->GetLabel()<0) nFakes++;
1661 // Store information in the output tree
1662 if (trackLoopIndex>-1) {
1663 recTrack = esdEvent->GetTrack(trackLoopIndex);
1664 } else if (trackIndex >-1) {
1665 recTrack = esdEvent->GetTrack(trackIndex);
1667 recTrack = fDummyTrack;
1670 particleMother = GetMother(particle,stack);
1671 mech = particle->GetUniqueID();
1673 //MC particle track length
1674 Double_t tpcTrackLength = 0.;
1675 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
1678 tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
1683 if(fTreeSRedirector && fFillTree) {
1684 (*fTreeSRedirector)<<"MCEffTree"<<
1685 "fileName.="<<&fCurrentFileName<<
1686 "triggerClass.="<<&triggerClass<<
1687 "runNumber="<<runNumber<<
1688 "evtTimeStamp="<<evtTimeStamp<<
1689 "evtNumberInFile="<<evtNumberInFile<< //
1690 "Bz="<<bz<< // magnetic field
1691 "vtxESD.="<<vtxESD<< // vertex info
1693 "mult="<<mult<< // primary vertex 9whatewe found) multiplicity
1694 "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
1695 // important variables for the pile-up studies
1696 "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
1697 "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
1698 "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
1699 "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
1700 "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
1701 "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
1704 "isAcc0="<<isESDtrackCut<< // track accepted by ESD track cuts
1705 "isAcc1="<<isAccCuts<< // track accepted by acceptance cuts flag
1706 "esdTrack.="<<recTrack<< // reconstructed track (only the longest from the loopers)
1707 "isRec="<<isRec<< // track was reconstructed
1708 "tpcTrackLength="<<tpcTrackLength<< // track length in the TPC r projection
1709 "particle.="<<particle<< // particle properties
1710 "particleMother.="<<particleMother<< // particle mother
1711 "mech="<<mech<< // production mechanizm
1712 "nRec="<<nRec<< // how many times reconstruted
1713 "nFakes="<<nFakes<< // how many times reconstructed as a fake track
1717 //if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
1723 //_____________________________________________________________________________
1724 Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {
1726 // check if particle is Z > 1
1728 if (track->GetTPCNcls() < 60) return kFALSE;
1729 Double_t mom = track->GetInnerParam()->GetP();
1730 if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
1731 Float_t dca[2], bCov[3];
1732 track->GetImpactParameters(dca,bCov);
1735 Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
1737 if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
1742 //_____________________________________________________________________________
1743 void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
1746 // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates
1749 AliDebug(AliLog::kError, "esdEvent not available");
1753 // get selection cuts
1754 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1755 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1756 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1758 if(!evtCuts || !accCuts || !esdTrackCuts) {
1759 AliDebug(AliLog::kError, "cuts not available");
1763 // trigger selection
1764 Bool_t isEventTriggered = kTRUE;
1765 AliPhysicsSelection *physicsSelection = NULL;
1766 AliTriggerAnalysis* triggerAnalysis = NULL;
1769 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1772 if(evtCuts->IsTriggerRequired())
1775 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1777 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1778 if(!physicsSelection) return;
1779 //SetPhysicsTriggerSelection(physicsSelection);
1781 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1782 // set trigger (V0AND)
1783 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1784 if(!triggerAnalysis) return;
1785 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1789 // centrality determination
1790 Float_t centralityF = -1;
1791 AliCentrality *esdCentrality = esdEvent->GetCentrality();
1792 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1795 // get reconstructed vertex
1796 //const AliESDVertex* vtxESD = 0;
1797 AliESDVertex* vtxESD = 0;
1798 if(GetAnalysisMode() == kTPCAnalysisMode) {
1799 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1801 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1802 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1810 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1811 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1812 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1815 if(isEventOK && isEventTriggered) {
1817 // Dump the pt downscaled V0 into the tree
1819 Int_t ntracks = esdEvent->GetNumberOfTracks();
1820 Int_t nV0s = esdEvent->GetNumberOfV0s();
1821 Int_t run = esdEvent->GetRunNumber();
1822 Int_t time= esdEvent->GetTimeStamp();
1823 Int_t evNr=esdEvent->GetEventNumberInFile();
1824 Double_t bz = esdEvent->GetMagneticField();
1826 // Global event id calculation using orbitID, bunchCrossingID and periodID
1827 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1828 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1829 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1830 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1833 for (Int_t iv0=0; iv0<nV0s; iv0++){
1835 AliESDv0 * v0 = esdEvent->GetV0(iv0);
1837 AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
1838 AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
1839 if (!track0) continue;
1840 if (!track1) continue;
1841 AliESDfriendTrack* friendTrack0=NULL;
1842 AliESDfriendTrack* friendTrack1=NULL;
1844 if (!esdFriend->TestSkipBit()){
1845 Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
1846 if (v0->GetIndex(0)<ntracksFriend){
1847 friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
1849 if (v0->GetIndex(1)<ntracksFriend){
1850 friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
1854 if (track0->GetSign()<0) {
1855 track1 = esdEvent->GetTrack(v0->GetIndex(0));
1856 track0 = esdEvent->GetTrack(v0->GetIndex(1));
1860 AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
1861 AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
1862 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
1863 if (gRandom->Rndm()>1./fFriendDownscaling){
1864 friendTrackStore0 = 0;
1865 friendTrackStore1 = 0;
1868 if (fFriendDownscaling<=0){
1869 if (((*fTreeSRedirector)<<"V0s").GetTree()){
1870 TTree * tree = ((*fTreeSRedirector)<<"V0s").GetTree();
1872 Double_t sizeAll=tree->GetZipBytes();
1873 TBranch * br= tree->GetBranch("friendTrack0.fPoints");
1874 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
1875 br= tree->GetBranch("friendTrack0.fCalibContainer");
1876 if (br) sizeFriend+=br->GetZipBytes();
1877 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
1878 friendTrackStore0=0;
1879 friendTrackStore1=0;
1886 Bool_t isDownscaled = IsV0Downscaled(v0);
1887 if (isDownscaled) continue;
1888 AliKFParticle kfparticle; //
1889 Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
1890 if (type==0) continue;
1891 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1893 if(!fFillTree) return;
1894 if(!fTreeSRedirector) return;
1895 (*fTreeSRedirector)<<"V0s"<<
1896 "gid="<<gid<< // global id of event
1897 "isDownscaled="<<isDownscaled<< //
1898 "triggerClass="<<&triggerClass<< // trigger
1900 "fileName.="<<&fCurrentFileName<< // full path - file name with ESD
1901 "runNumber="<<run<< //
1902 "evtTimeStamp="<<time<< // time stamp of event in secons
1903 "evtNumberInFile="<<evNr<< //
1904 "type="<<type<< // type of V0-
1905 "ntracks="<<ntracks<<
1907 "kf.="<<&kfparticle<<
1908 "track0.="<<track0<< // track
1909 "track1.="<<track1<<
1910 "friendTrack0.="<<friendTrackStore0<<
1911 "friendTrack1.="<<friendTrackStore1<<
1912 "centralityF="<<centralityF<<
1918 //_____________________________________________________________________________
1919 void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
1922 // Select real events with large TPC dEdx signal
1925 AliDebug(AliLog::kError, "esdEvent not available");
1929 // get selection cuts
1930 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1931 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1932 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1934 if(!evtCuts || !accCuts || !esdTrackCuts) {
1935 AliDebug(AliLog::kError, "cuts not available");
1941 Bool_t isEventTriggered = kTRUE;
1942 AliPhysicsSelection *physicsSelection = NULL;
1943 AliTriggerAnalysis* triggerAnalysis = NULL;
1946 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1948 if(evtCuts->IsTriggerRequired())
1951 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1953 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1954 if(!physicsSelection) return;
1956 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1957 // set trigger (V0AND)
1958 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1959 if(!triggerAnalysis) return;
1960 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1964 // get reconstructed vertex
1965 AliESDVertex* vtxESD = 0;
1966 if(GetAnalysisMode() == kTPCAnalysisMode) {
1967 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1969 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1970 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1977 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1978 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1979 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1983 if(isEventOK && isEventTriggered)
1985 Double_t vert[3] = {0};
1986 vert[0] = vtxESD->GetX();
1987 vert[1] = vtxESD->GetY();
1988 vert[2] = vtxESD->GetZ();
1989 Int_t mult = vtxESD->GetNContributors();
1990 Double_t bz = esdEvent->GetMagneticField();
1991 Double_t runNumber = esdEvent->GetRunNumber();
1992 Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1993 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1995 // Global event id calculation using orbitID, bunchCrossingID and periodID
1996 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1997 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1998 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1999 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
2002 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
2004 AliESDtrack *track = esdEvent->GetTrack(iTrack);
2005 if(!track) continue;
2006 AliESDfriendTrack* friendTrack=NULL;
2007 if (esdFriend && !esdFriend->TestSkipBit()) {
2008 Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
2009 if (iTrack<ntracksFriend){
2010 friendTrack = esdFriend->GetTrack(iTrack);
2011 } //this guy can be NULL
2014 if(track->Charge()==0) continue;
2015 if(!esdTrackCuts->AcceptTrack(track)) continue;
2016 if(!accCuts->AcceptTrack(track)) continue;
2018 if(!IsHighDeDxParticle(track)) continue;
2019 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2021 if(!fFillTree) return;
2022 if(!fTreeSRedirector) return;
2023 (*fTreeSRedirector)<<"dEdx"<< // high dEdx tree
2024 "gid="<<gid<< // global id
2025 "fileName.="<<&fCurrentFileName<< // file name
2026 "runNumber="<<runNumber<<
2027 "evtTimeStamp="<<evtTimeStamp<<
2028 "evtNumberInFile="<<evtNumberInFile<<
2029 "triggerClass="<<&triggerClass<< // trigger
2031 "vtxESD.="<<vtxESD<< //
2033 "esdTrack.="<<track<<
2034 "friendTrack.="<<friendTrack<<
2040 //_____________________________________________________________________________
2041 Int_t AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
2044 // Create KF particle in case the V0 fullfill selection criteria
2046 // Selection criteria
2050 // 4. rough mass cut
2051 // 5. Normalized pointing angle cut
2053 const Double_t cutMass=0.2;
2054 const Double_t kSigmaDCACut=3;
2056 // 0.) algo cut - accept only on the fly
2058 if (v0->GetOnFlyStatus() ==kFALSE) return 0;
2062 AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
2063 AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
2065 TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
2066 TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
2067 TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
2069 if (TMath::Abs(track0->GetTgl())>1) return 0;
2070 if (TMath::Abs(track1->GetTgl())>1) return 0;
2071 if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
2072 if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
2073 Float_t pos0[2]={0}, cov0[3]={0};
2074 Float_t pos1[2]={0}, cov1[3]={0};
2075 track0->GetImpactParameters(pos0,cov0);
2076 track0->GetImpactParameters(pos1,cov1);
2078 if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
2079 if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
2084 Double_t chi2KF = v0->GetKFInfo(2,2,2);
2085 if (chi2KF>25) return 0;
2087 // 4.) Rough mass cut - 0.200 GeV
2089 static Double_t masses[2]={-1};
2091 masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
2092 masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
2094 Double_t mass00= v0->GetEffMass(0,0);
2095 Double_t mass22= v0->GetEffMass(2,2);
2096 Double_t mass42= v0->GetEffMass(4,2);
2097 Double_t mass24= v0->GetEffMass(2,4);
2098 Bool_t massOK=kFALSE;
2103 if (TMath::Abs(mass00-0)<cutMass) {
2104 massOK=kTRUE; type+=1;
2105 if (TMath::Abs(mass00-0)<dmass) {
2107 dmass=TMath::Abs(mass00-0);
2111 if (TMath::Abs(mass24-masses[1])<cutMass) {
2112 massOK=kTRUE; type+=2;
2113 if (TMath::Abs(mass24-masses[1])<dmass){
2114 dmass = TMath::Abs(mass24-masses[1]);
2119 if (TMath::Abs(mass42-masses[1])<cutMass) {
2120 massOK=kTRUE; type+=4;
2121 if (TMath::Abs(mass42-masses[1])<dmass){
2122 dmass = TMath::Abs(mass42-masses[1]);
2127 if (TMath::Abs(mass22-masses[0])<cutMass) {
2128 massOK=kTRUE; type+=8;
2129 if (TMath::Abs(mass22-masses[0])<dmass){
2130 dmass = TMath::Abs(mass22-masses[0]);
2135 if (type==0) return 0;
2137 const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
2138 const AliExternalTrackParam *paramP = v0->GetParamP();
2139 const AliExternalTrackParam *paramN = v0->GetParamN();
2140 if (paramP->GetSign()<0){
2141 paramP=v0->GetParamP();
2142 paramN=v0->GetParamN();
2144 //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
2145 //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
2147 AliKFParticle kfp1( *paramP, spdg[p1] );
2148 AliKFParticle kfp2( *paramN, -1 *spdg[p2] );
2156 Double_t errPhi = V0KF.GetErrPhi();
2157 Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
2158 if (pointAngle/errPhi>10) return 0;
2163 //_____________________________________________________________________________
2164 Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)
2167 // Downscale randomly low pt V0
2170 Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
2171 Double_t scalempt= TMath::Min(maxPt,10.);
2172 Double_t downscaleF = gRandom->Rndm();
2173 downscaleF *= fLowPtV0DownscaligF;
2175 // Special treatment of the gamma conversion pt spectra is softer -
2176 Double_t mass00= v0->GetEffMass(0,0);
2177 const Double_t cutMass=0.2;
2178 if (TMath::Abs(mass00-0)<cutMass){
2179 downscaleF/=10.; // 10 times smaller downscaling for the gamma concersion candidate
2181 //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
2182 if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;
2186 TH1F his1("his1","his1",100,0,10);
2187 TH1F his2("his2","his2",100,0,10);
2188 {for (Int_t i=0; i<10000; i++){
2189 Double_t rnd=gRandom->Exp(1);
2190 Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
2192 if (!isDownscaled) his2->Fill(rnd);
2201 //_____________________________________________________________________________
2202 Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
2204 // Constrain TPC inner params constrained
2206 if(!tpcInnerC) return kFALSE;
2207 if(!vtx) return kFALSE;
2209 Double_t dz[2],cov[3];
2210 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2211 //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2212 //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2213 if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
2216 Double_t covar[6]; vtx->GetCovMatrix(covar);
2217 Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
2218 Double_t c[3]={covar[2],0.,covar[5]};
2219 Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
2220 if (chi2C>kVeryBig) return kFALSE;
2222 if(!tpcInnerC->Update(p,c)) return kFALSE;
2227 //_____________________________________________________________________________
2228 Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
2230 // Constrain TPC inner params constrained
2232 if(!trackInnerC) return kFALSE;
2233 if(!vtx) return kFALSE;
2235 const Double_t kRadius = 2.8;
2236 const Double_t kMaxStep = 1.0;
2238 Double_t dz[2],cov[3];
2240 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2241 //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2242 //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2244 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
2245 if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
2248 Double_t covar[6]; vtx->GetCovMatrix(covar);
2249 Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
2250 Double_t c[3]={covar[2],0.,covar[5]};
2251 Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
2252 if (chi2C>kVeryBig) return kFALSE;
2254 if(!trackInnerC->Update(p,c)) return kFALSE;
2260 //_____________________________________________________________________________
2261 TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack)
2263 if(!particle) return NULL;
2264 if(!stack) return NULL;
2266 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2267 TParticle* mother = NULL;
2268 Int_t mcStackSize=stack->GetNtrack();
2269 if (motherLabel>=mcStackSize) return NULL;
2270 mother = stack->Particle(motherLabel);
2275 //_____________________________________________________________________________
2276 Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(Int_t label, AliStack *const stack)
2278 Bool_t isFromConversion = kFALSE;
2281 Int_t mcStackSize=stack->GetNtrack();
2282 if (label>=mcStackSize) return kFALSE;
2283 TParticle* particle = stack->Particle(label);
2284 if (!particle) return kFALSE;
2286 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2288 Int_t mech = particle->GetUniqueID(); // production mechanism
2289 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2291 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2292 if (motherLabel>=mcStackSize) return kFALSE;
2293 TParticle* mother = stack->Particle(motherLabel);
2295 Int_t motherPdg = mother->GetPdgCode();
2297 if(!isPrim && mech==5 && motherPdg==kGamma) {
2298 isFromConversion=kTRUE;
2304 return isFromConversion;
2307 //_____________________________________________________________________________
2308 Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(Int_t label, AliStack *const stack)
2310 Bool_t isFromMaterial = kFALSE;
2313 Int_t mcStackSize=stack->GetNtrack();
2314 if (label>=mcStackSize) return kFALSE;
2315 TParticle* particle = stack->Particle(label);
2316 if (!particle) return kFALSE;
2318 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2320 Int_t mech = particle->GetUniqueID(); // production mechanism
2321 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2323 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2324 if (motherLabel>=mcStackSize) return kFALSE;
2325 TParticle* mother = stack->Particle(motherLabel);
2327 if(!isPrim && mech==13) {
2328 isFromMaterial=kTRUE;
2334 return isFromMaterial;
2337 //_____________________________________________________________________________
2338 Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(Int_t label, AliStack *const stack)
2340 Bool_t isFromStrangeness = kFALSE;
2343 Int_t mcStackSize=stack->GetNtrack();
2344 if (label>=mcStackSize) return kFALSE;
2345 TParticle* particle = stack->Particle(label);
2346 if (!particle) return kFALSE;
2348 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2350 Int_t mech = particle->GetUniqueID(); // production mechanism
2351 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2353 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2354 if (motherLabel>=mcStackSize) return kFALSE;
2355 TParticle* mother = stack->Particle(motherLabel);
2357 Int_t motherPdg = mother->GetPdgCode();
2359 // K+-, lambda, antilambda, K0s decays
2360 if(!isPrim && mech==4 &&
2361 (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
2363 isFromStrangeness = kTRUE;
2369 return isFromStrangeness;
2373 //_____________________________________________________________________________
2374 void AliAnalysisTaskFilteredTree::FinishTaskOutput()
2377 // Called one at the end
2378 // locally on working node
2380 Bool_t deleteTrees=kTRUE;
2381 if ((AliAnalysisManager::GetAnalysisManager()))
2383 if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
2384 AliAnalysisManager::kProofAnalysis)
2387 if (deleteTrees) delete fTreeSRedirector;
2388 fTreeSRedirector=NULL;
2391 //_____________________________________________________________________________
2392 void AliAnalysisTaskFilteredTree::Terminate(Option_t *)
2395 // Called one at the end
2399 //_____________________________________________________________________________
2400 Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
2403 // calculate mc event true track multiplicity
2405 if(!mcEvent) return 0;
2407 AliStack* stack = 0;
2410 // MC particle stack
2411 stack = mcEvent->Stack();
2412 if (!stack) return 0;
2415 //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
2418 Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
2419 if(!isEventOK) return 0;
2421 Int_t nPart = stack->GetNtrack();
2422 for (Int_t iMc = 0; iMc < nPart; ++iMc)
2424 TParticle* particle = stack->Particle(iMc);
2428 // only charged particles
2429 if(!particle->GetPDG()) continue;
2430 Double_t charge = particle->GetPDG()->Charge()/3.;
2431 if (TMath::Abs(charge) < 0.001)
2435 Bool_t prim = stack->IsPhysicalPrimary(iMc);
2438 // checked accepted without pt cut
2439 //if(accCuts->AcceptTrack(particle))
2440 if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
2449 //_____________________________________________________________________________
2450 void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC)
2453 // Fill pT relative resolution histograms for
2454 // TPC only, TPC only constrained to vertex and TPC+ITS tracking
2457 if(!ptpcInnerC) return;
2459 const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();
2460 if(!innerParam) return;
2463 ptrack->GetImpactParameters(dxy,dz);
2465 // TPC+ITS primary tracks
2466 if( abs(ptrack->Eta())<0.8 &&
2467 ptrack->GetTPCClusterInfo(3,1)>120 &&
2468 ptrack->IsOn(0x40) &&
2469 ptrack->GetTPCclusters(0)>0.0 &&
2470 ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2471 //abs(innerParam->GetX())>0.0 &&
2472 //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2473 //abs(innerParam->GetTgl())<0.85 &&
2474 ptrack->IsOn(0x0004) &&
2475 ptrack->GetNcls(0)>0 &&
2476 ptrack->GetITSchi2()>0 &&
2477 sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 &&
2478 sqrt(chi2TPCInnerC)<6 &&
2479 (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) &&
2481 abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )
2483 fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2484 fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2485 fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2488 // TPC primary tracks
2489 // and TPC constrained primary tracks
2491 AliExternalTrackParam *ptpcInner = (AliExternalTrackParam *) ptrack->GetTPCInnerParam();
2492 if(!ptpcInner) return;
2495 Float_t dxyTPC, dzTPC;
2496 ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);
2498 if( abs(ptrack->Eta())<0.8 &&
2499 ptrack->GetTPCClusterInfo(3,1)>120 &&
2500 ptrack->IsOn(0x40)&&
2501 ptrack->GetTPCclusters(0)>0.0 &&
2502 ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2503 //abs(innerParam->GetX())>0.0 &&
2504 //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2505 //abs(innerParam->GetTgl())<0.85 &&
2510 fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2511 fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2512 fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2514 // TPC constrained to vertex
2515 fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2516 fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2517 fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2522 void AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend){
2524 // Process ITS standalone tracks find match with closest TPC(or combined tracks) tracks
2525 // marian.ivanov@cern.ch
2526 // 0.) Init variables
2527 // 1.) GetTrack parameters at TPC inner wall
2528 // 2.) Match closest TPC track (STANDALONE/global) - chi2 match criteria
2530 // Logic to be used in reco:
2531 // 1.) Find matching ITSalone->TPCalone
2532 // 2.) if (!TPCalone.FindClose(TPCother)) TPCalone.Addopt(ITSalone)
2533 // 3.) ff ((ITSalone.FindClose(Global)==0) CreateGlobaltrack
2534 const Double_t radiusMatch=84.; // redius to propagate
2536 const Double_t dFastPhiCut=0.2; // 6 sigma (200 MeV) fast angular cut
2537 const Double_t dFastThetaCut=0.12; // 6 sigma (200 MeV) fast angular cut
2538 const Double_t dFastPosPhiCut=0.06; // 6 sigma (200 MeV) fast angular cut
2539 const Double_t dFastZCut=6; // 6 sigma (200 MeV) fast z difference cut
2540 const Double_t dFastPtCut=2.; // 6 sigma (200 MeV) fast 1/pt cut
2541 const Double_t chi2Cut=100; // chi2 matching cut
2543 if (!esdFriend) return; // not ITS standalone track
2544 if (esdFriend->TestSkipBit()) return; // friends tracks not stored
2545 Int_t ntracks=esdEvent->GetNumberOfTracks();
2546 Float_t bz = esdEvent->GetMagneticField();
2548 // 0.) Get parameters in reference radius TPC Inner wall
2551 TMatrixD vecPosR0(ntracks,6); // possition and momentum estimate at reference radius
2552 TMatrixD vecMomR0(ntracks,6); //
2553 TMatrixD vecPosR1(ntracks,6); // possition and momentum estimate at reference radius TPC track
2554 TMatrixD vecMomR1(ntracks,6); //
2555 Double_t xyz[3], pxyz[3]; //
2556 for (Int_t iTrack=0; iTrack<ntracks; iTrack++){
2557 AliESDtrack *track = esdEvent->GetTrack(iTrack);
2558 if(!track) continue;
2559 if (track->GetInnerParam()){
2560 const AliExternalTrackParam *trackTPC=track->GetInnerParam();
2561 trackTPC->GetXYZAt(radiusMatch,bz,xyz);
2562 trackTPC->GetPxPyPzAt(radiusMatch,bz,pxyz);
2563 for (Int_t i=0; i<3; i++){
2564 vecPosR1(iTrack,i)=xyz[i];
2565 vecMomR1(iTrack,i)=pxyz[i];
2567 vecPosR1(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); // phi pos angle
2568 vecMomR1(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); // phi mom angle
2569 vecMomR1(iTrack,4)= trackTPC->GetSigned1Pt();;
2570 vecMomR1(iTrack,5)= trackTPC->GetTgl();;
2572 AliESDfriendTrack* friendTrack=esdFriend->GetTrack(iTrack);
2573 if(!friendTrack) continue;
2574 if (friendTrack->GetITSOut()){
2575 const AliExternalTrackParam *trackITS=friendTrack->GetITSOut();
2576 trackITS->GetXYZAt(radiusMatch,bz,xyz);
2577 trackITS->GetPxPyPzAt(radiusMatch,bz,pxyz);
2578 for (Int_t i=0; i<3; i++){
2579 vecPosR0(iTrack,i)=xyz[i];
2580 vecMomR0(iTrack,i)=pxyz[i];
2582 vecPosR0(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]);
2583 vecMomR0(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]);
2584 vecMomR0(iTrack,4)= trackITS->GetSigned1Pt();;
2585 vecMomR0(iTrack,5)= trackITS->GetTgl();;
2589 // 1.) Find closest matching tracks, between the ITS standalone track
2590 // and the all other tracks
2591 // a.) caltegory - All
2592 // b.) category - without ITS
2595 Int_t ntracksPropagated=0;
2596 AliExternalTrackParam extTrackDummy;
2597 AliESDtrack esdTrackDummy;
2598 AliExternalTrackParam itsAtTPC;
2599 AliExternalTrackParam itsAtITSTPC;
2600 for (Int_t iTrack0=0; iTrack0<ntracks; iTrack0++){
2601 AliESDtrack *track0 = esdEvent->GetTrack(iTrack0);
2602 if(!track0) continue;
2603 if (track0->IsOn(AliVTrack::kTPCin)) continue;
2604 AliESDfriendTrack* friendTrack0=esdFriend->GetTrack(iTrack0);
2605 if (!friendTrack0) continue;
2606 //if (!track0->IsOn(AliVTrack::kITSpureSA)) continue;
2607 //if (!friendTrack0->GetITSOut()) continue; // is there flag for ITS standalone?
2608 ntracksPropagated++;
2610 // 2.) find clostest TPCtrack
2612 Double_t minChi2All=10000000;
2613 Double_t minChi2TPC=10000000;
2614 Double_t minChi2TPCITS=10000000;
2617 Int_t indexTPCITS=-1;
2618 Int_t ncandidates0=0; // n candidates - rough cut
2619 Int_t ncandidates1=0; // n candidates - rough + chi2 cut
2620 itsAtTPC=*(friendTrack0->GetITSOut());
2621 itsAtITSTPC=*(friendTrack0->GetITSOut());
2622 for (Int_t iTrack1=0; iTrack1<ntracks; iTrack1++){
2623 AliESDtrack *track1 = esdEvent->GetTrack(iTrack1);
2624 if(!track1) continue;
2625 if (!track1->IsOn(AliVTrack::kTPCin)) continue;
2628 if (TMath::Abs(vecPosR1(iTrack1,2)-vecPosR0(iTrack0,2))>dFastZCut) continue;
2629 if (TMath::Abs(vecPosR1(iTrack1,3)-vecPosR0(iTrack0,3))>dFastPosPhiCut) continue;
2630 if (TMath::Abs(vecMomR1(iTrack1,3)-vecMomR0(iTrack0,3))>dFastPhiCut) continue;
2631 if (TMath::Abs(vecMomR1(iTrack1,5)-vecMomR0(iTrack0,5))>dFastThetaCut) continue;
2632 if (TMath::Abs(vecMomR1(iTrack1,4)-vecMomR0(iTrack0,4))>dFastPtCut) continue;
2635 const AliExternalTrackParam * param1= track1->GetInnerParam();
2636 if (!friendTrack0->GetITSOut()) continue;
2637 AliExternalTrackParam outerITS = *(friendTrack0->GetITSOut());
2638 if (!outerITS.Rotate(param1->GetAlpha())) continue;
2639 if (!outerITS.PropagateTo(param1->GetX(),bz)) continue; // assume track close to the TPC inner wall
2640 Double_t chi2 = outerITS.GetPredictedChi2(param1);
2641 if (chi2>chi2Cut) continue;
2643 if (chi2<minChi2All){
2647 if (chi2<minChi2TPC && track1->IsOn(AliVTrack::kITSin)==0){
2652 if (chi2<minChi2TPCITS && track1->IsOn(AliVTrack::kITSin)){
2654 indexTPCITS=iTrack1;
2655 itsAtITSTPC=outerITS;
2659 AliESDtrack * trackAll= (indexAll>=0)? esdEvent->GetTrack(indexAll):&esdTrackDummy;
2660 AliESDtrack * trackTPC= (indexTPC>=0)? esdEvent->GetTrack(indexTPC):&esdTrackDummy;
2661 AliESDtrack * trackTPCITS= (indexTPCITS>=0)? esdEvent->GetTrack(indexTPCITS):&esdTrackDummy;
2662 (*fTreeSRedirector)<<"itsTPC"<<
2663 "indexAll="<<indexAll<< // index of closest track (chi2)
2664 "indexTPC="<<indexTPC<< // index of closest TPCalone tracks
2665 "indexTPCITS="<<indexTPCITS<< // index of closest cobined tracks
2666 "ncandidates0="<<ncandidates0<< // number of candidates
2667 "ncandidates1="<<ncandidates1<<
2669 "chi2All="<<minChi2All<< // chi2 of closest tracks
2670 "chi2TPC="<<minChi2TPC<<
2671 "chi2TPCITS="<<minChi2TPCITS<<
2673 "track0.="<<track0<< // ITS standalone tracks
2674 "trackAll.="<<trackAll<< // Closets other track
2675 "trackTPC.="<<trackTPC<< // Closest TPC only track
2676 "trackTPCITS.="<<trackTPCITS<< // closest combined track
2678 "itsAtTPC.="<<&itsAtTPC<< // ITS track parameters at the TPC alone track frame
2679 "itsAtITSTPC.="<<&itsAtITSTPC<< // ITS track parameters at the TPC combeined track frame
2684 void AliAnalysisTaskFilteredTree::ProcessTrackMatch(AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/){
2686 Use TPC standalone, ITS standalone and combined tracks to categorize/resp. recover track information.
2692 Topology categories:
2693 -Nice isolated tracks with full information
2694 -Overlapped tracks - Refit and sign them
2695 -Multiple found (check overlap factor) - Merge and sign
2696 -Charge particle (kink) decays - Change of direction - Sign them)
2698 -Array of indexes of closest tracks in each individual category
2699 -Chi2 distance to the closest tracks at reference radius of TPCin
2700 -Overlap factors - fraction of shared clusters or missing region
2701 -Chi2 distance at DCA
2703 -matrix closest tracks from each categories
2704 -characterization - chi2, index,chi2, overlap ratio
2707 0.) Kink decay or catastophic multiple scaterring
2708 (combining all track categories)
2709 - small chi2 at the DCA
2710 - significantly large deflection angle
2711 - Combinatorial algorithm - to decrease CPU time restriction of investigation to tracks with small DCA at refernce radius
2713 1.) if (TPConly && !(TPC+ITS) && ITSonly match ) {
2714 if (!kink) TPCOnly.addoptITS
2715 if (kink) TPConly sign
2718 2.) Overlap tracks - Refit with doUnfold