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", Entry());
355 //_____________________________________________________________________________
356 void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliESDfriend* esdFriend)
359 // Find cosmic pairs (triggered or random)
362 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
363 AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
364 AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC();
365 const Double_t kMinPt=0.8;
366 const Double_t kMinPtMax=0.8;
367 const Double_t kMinNcl=50;
368 const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
369 Int_t ntracks=event->GetNumberOfTracks();
370 UInt_t specie = event->GetEventSpecie(); // skip laser events
371 if (specie==AliRecoParam::kCalib) return;
375 for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
376 AliESDtrack *track0 = event->GetTrack(itrack0);
377 if (!track0) continue;
378 if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
380 if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue;
381 if (track0->Pt() < kMinPt) continue;
382 if (track0->GetTPCncls() < kMinNcl) continue;
383 if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
384 if (track0->GetKinkIndex(0)>0) continue;
385 const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
388 //track0->GetImpactParametersTPC(dcaTPC,covTPC);
389 //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
390 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
391 // const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
392 AliESDfriendTrack* friendTrack0=NULL;
393 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack0 = esdFriend->GetTrack(itrack0);} //this guy can be NULL
395 for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
396 AliESDtrack *track1 = event->GetTrack(itrack1);
397 if (!track1) continue;
398 if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
399 if (track1->GetKinkIndex(0)>0) continue;
400 if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue;
401 if (track1->Pt() < kMinPt) continue;
402 if (track1->GetTPCncls()<kMinNcl) continue;
403 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
404 if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
405 //track1->GetImpactParametersTPC(dcaTPC,covTPC);
406 // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
407 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
409 const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
412 for (Int_t ipar=0; ipar<5; ipar++){
413 if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
414 if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
416 if (!isPair) continue;
417 if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
418 //delta with correct sign
420 TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
421 TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
422 TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
424 if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
425 if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
426 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
427 if (!isPair) continue;
428 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
429 Int_t eventNumber = event->GetEventNumberInFile();
432 Int_t ntracksSPD = vertexSPD->GetNContributors();
433 Int_t ntracksTPC = vertexTPC->GetNContributors();
434 Int_t runNumber = event->GetRunNumber();
435 Int_t timeStamp = event->GetTimeStamp();
436 ULong64_t triggerMask = event->GetTriggerMask();
437 Float_t magField = event->GetMagneticField();
438 TObjString triggerClass = event->GetFiredTriggerClasses().Data();
440 // Global event id calculation using orbitID, bunchCrossingID and periodID
441 ULong64_t orbitID = (ULong64_t)event->GetOrbitNumber();
442 ULong64_t bunchCrossID = (ULong64_t)event->GetBunchCrossNumber();
443 ULong64_t periodID = (ULong64_t)event->GetPeriodNumber();
444 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
447 AliESDfriendTrack* friendTrack1=NULL;
448 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack1 = esdFriend->GetTrack(itrack1);} //this guy can be NULL
450 AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
451 AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
452 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
453 if (gRandom->Rndm()>1./fFriendDownscaling){
454 friendTrackStore0 = 0;
455 friendTrackStore1 = 0;
458 if (fFriendDownscaling<=0){
459 if (((*fTreeSRedirector)<<"CosmicPairs").GetTree()){
460 TTree * tree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
462 Double_t sizeAll=tree->GetZipBytes();
463 TBranch * br= tree->GetBranch("friendTrack0.fPoints");
464 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
465 br= tree->GetBranch("friendTrack0.fCalibContainer");
466 if (br) sizeFriend+=br->GetZipBytes();
467 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
474 if(!fFillTree) return;
475 if(!fTreeSRedirector) return;
476 (*fTreeSRedirector)<<"CosmicPairs"<<
477 "gid="<<gid<< // global id of track
478 "fileName.="<<&fCurrentFileName<< // file name
479 "runNumber="<<runNumber<< // run number
480 "evtTimeStamp="<<timeStamp<< // time stamp of event
481 "evtNumberInFile="<<eventNumber<< // event number
482 "trigger="<<triggerMask<< // trigger mask
483 "triggerClass="<<&triggerClass<< // trigger class
484 "Bz="<<magField<< // magnetic field
486 "multSPD="<<ntracksSPD<< // event ultiplicity
487 "multTPC="<<ntracksTPC<< //
488 "vertSPD.="<<vertexSPD<< // primary vertex -SPD
489 "vertTPC.="<<vertexTPC<< // primary vertex -TPC
490 "t0.="<<track0<< // first half of comsic trak
491 "t1.="<<track1<< // second half of cosmic track
492 "friendTrack0.="<<friendTrackStore0<< // friend information first track + points
493 "friendTrack1.="<<friendTrackStore1<< // frined information first track + points
500 //_____________________________________________________________________________
501 void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
504 // Select real events with high-pT tracks
507 // get selection cuts
508 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
509 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
510 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
512 if(!evtCuts || !accCuts || !esdTrackCuts) {
513 Printf("ERROR cuts not available");
518 Bool_t isEventTriggered = kTRUE;
519 AliPhysicsSelection *physicsSelection = NULL;
520 AliTriggerAnalysis* triggerAnalysis = NULL;
523 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
525 if(evtCuts->IsTriggerRequired())
528 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
530 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
531 if(!physicsSelection) return;
532 //SetPhysicsTriggerSelection(physicsSelection);
534 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
535 // set trigger (V0AND)
536 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
537 if(!triggerAnalysis) return;
538 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
542 // centrality determination
543 Float_t centralityF = -1;
544 AliCentrality *esdCentrality = esdEvent->GetCentrality();
545 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
547 // use MC information
548 AliHeader* header = 0;
549 AliGenEventHeader* genHeader = 0;
553 Int_t multMCTrueTracks = 0;
556 // get MC event header
557 header = mcEvent->Header();
559 AliDebug(AliLog::kError, "Header not available");
563 stack = mcEvent->Stack();
565 AliDebug(AliLog::kError, "Stack not available");
570 genHeader = header->GenEventHeader();
572 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
575 genHeader->PrimaryVertex(vtxMC);
577 // multipliticy of all MC primary tracks
578 // in Zv, pt and eta ranges)
579 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
582 // get reconstructed vertex
583 //const AliESDVertex* vtxESD = 0;
584 AliESDVertex* vtxESD = 0;
585 if(GetAnalysisMode() == kTPCAnalysisMode) {
586 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
588 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
589 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
595 AliESDVertex* vtxTPC = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
596 AliESDVertex* vtxSPD = (AliESDVertex*)esdEvent->GetPrimaryVertexSPD();
602 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
603 //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZv());
604 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
605 Int_t ntracks = esdEvent->GetNumberOfTracks();
608 if(isEventOK && isEventTriggered)
612 // get IR information
614 AliESDHeader *esdHeader = 0;
615 esdHeader = esdEvent->GetHeader();
616 if(!esdHeader) return;
617 //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
618 //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
620 // Use when Peter commit the changes in the header
625 //Double_t vert[3] = {0};
626 //vert[0] = vtxESD->GetXv();
627 //vert[1] = vtxESD->GetYv();
628 //vert[2] = vtxESD->GetZv();
629 Int_t mult = vtxESD->GetNContributors();
630 Int_t multSPD = vtxSPD->GetNContributors();
631 Int_t multTPC = vtxTPC->GetNContributors();
633 Float_t bz = esdEvent->GetMagneticField();
634 Int_t runNumber = esdEvent->GetRunNumber();
635 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
636 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
638 // Global event id calculation using orbitID, bunchCrossingID and periodID
639 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
640 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
641 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
642 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
646 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
648 AliESDtrack *track = esdEvent->GetTrack(iTrack);
650 if(track->Charge()==0) continue;
651 if(!esdTrackCuts->AcceptTrack(track)) continue;
652 if(!accCuts->AcceptTrack(track)) continue;
654 // downscale low-pT tracks
655 Double_t scalempt= TMath::Min(track->Pt(),10.);
656 Double_t downscaleF = gRandom->Rndm();
657 downscaleF *= fLowPtTrackDownscaligF;
658 if(TMath::Exp(2*scalempt)<downscaleF) continue;
659 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
661 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
662 if (!tpcInner) continue;
663 // transform to the track reference frame
664 Bool_t isOK = kFALSE;
665 isOK = tpcInner->Rotate(track->GetAlpha());
666 isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
673 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
674 if(!fFillTree) return;
675 if(!fTreeSRedirector) return;
676 (*fTreeSRedirector)<<"highPt"<<
678 "fileName.="<<&fCurrentFileName<<
679 "runNumber="<<runNumber<<
680 "evtTimeStamp="<<evtTimeStamp<<
681 "evtNumberInFile="<<evtNumberInFile<<
682 "triggerClass="<<&triggerClass<< // trigger
683 "Bz="<<bz<< // magnetic field
685 "ntracksESD="<<ntracks<< // number of tracks in the ESD
686 "IRtot="<<ir1<< // interaction record history info
688 "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
689 "multSPD="<<multSPD<< // multiplicity of tracks pointing to the SPD primary vertex
690 "multTPC="<<multTPC<< // multiplicity of tracks pointing to the TPC primary vertex
691 "esdTrack.="<<track<<
692 "centralityF="<<centralityF<<
700 //_____________________________________________________________________________
701 void AliAnalysisTaskFilteredTree::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const esdFriend)
704 // Process laser events -> dump tracks and clusters to the special tree
706 if(!fFillTree) return;
707 if(!fTreeSRedirector) return;
708 const AliESDHeader* esdHeader = esdEvent->GetHeader();
709 if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) {
710 Int_t countLaserTracks = 0;
711 Int_t runNumber = esdEvent->GetRunNumber();
712 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
713 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
714 Float_t bz = esdEvent->GetMagneticField();
715 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
716 // Global event id calculation using orbitID, bunchCrossingID and periodID
717 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
718 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
719 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
720 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
721 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
722 AliESDtrack *track = esdEvent->GetTrack(iTrack);
724 if(track->GetTPCInnerParam()) countLaserTracks++;
725 AliESDfriendTrack* friendTrack=NULL;
726 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
727 (*fTreeSRedirector)<<"Laser"<<
728 "gid="<<gid<< // global identifier of event
729 "fileName.="<<&fCurrentFileName<< //
730 "runNumber="<<runNumber<<
731 "evtTimeStamp="<<evtTimeStamp<<
732 "evtNumberInFile="<<evtNumberInFile<<
733 "triggerClass="<<&triggerClass<< // trigger
734 "Bz="<<bz<< // magnetic field
735 "multTPCtracks="<<countLaserTracks<< // multiplicity of tracks
736 "track=."<<track<< // track parameters
737 "friendTrack.="<<friendTrack<< // friend track information
743 //_____________________________________________________________________________
744 void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
747 // Select real events with high-pT tracks
748 // Calculate and stor in the output tree:
749 // TPC constrained tracks
750 // InnerParams constrained tracks
752 // ITSout-InnerParams tracks
753 // chi2 distance between TPC constrained and TPC-ITS tracks
754 // chi2 distance between TPC refitted constrained and TPC-ITS tracks
755 // chi2 distance between ITSout and InnerParams tracks
757 // track references at ITSin, TPCin; InnerParam at first TPC track reference,
758 // particle ID, mother ID, production mechanism ...
760 // get selection cuts
761 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
762 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
763 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
765 if(!evtCuts || !accCuts || !esdTrackCuts) {
766 AliDebug(AliLog::kError, "cuts not available");
771 Bool_t isEventTriggered = kTRUE;
772 AliPhysicsSelection *physicsSelection = NULL;
773 AliTriggerAnalysis* triggerAnalysis = NULL;
776 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
780 if(evtCuts->IsTriggerRequired())
783 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
785 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
786 if(!physicsSelection) {AliInfo("no physics selection"); return;}
787 //SetPhysicsTriggerSelection(physicsSelection);
789 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
790 // set trigger (V0AND)
791 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
792 if(!triggerAnalysis) {AliInfo("no trigger analysis");return;}
793 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
797 // centrality determination
798 Float_t centralityF = -1;
799 AliCentrality *esdCentrality = esdEvent->GetCentrality();
800 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
802 // use MC information
803 AliHeader* header = 0;
804 AliGenEventHeader* genHeader = 0;
809 Int_t multMCTrueTracks = 0;
812 // get MC event header
813 header = mcEvent->Header();
815 AliDebug(AliLog::kError, "Header not available");
819 stack = mcEvent->Stack();
821 AliDebug(AliLog::kError, "Stack not available");
824 mcStackSize=stack->GetNtrack();
827 genHeader = header->GenEventHeader();
829 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
832 genHeader->PrimaryVertex(vtxMC);
834 // multipliticy of all MC primary tracks
835 // in Zv, pt and eta ranges)
836 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
840 // get reconstructed vertex
841 //const AliESDVertex* vtxESD = 0;
842 AliESDVertex* vtxESD = 0;
843 if(GetAnalysisMode() == kTPCAnalysisMode) {
844 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
846 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
847 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
850 AliInfo("no ESD vertex");
856 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
859 // Vertex info comparison and track multiplicity
861 AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
862 AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
863 Int_t contSPD = vertexSPD->GetNContributors();
864 Int_t contTPC = vertexTPC->GetNContributors();
865 TVectorD vertexPosTPC(3), vertexPosSPD(3);
866 vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
867 vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
870 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
871 AliESDtrack *track = esdEvent->GetTrack(iTrack);
873 if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
874 if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
877 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
878 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
882 if(isEventOK && isEventTriggered)
885 // get IR information
887 AliESDHeader *esdHeader = 0;
888 esdHeader = esdEvent->GetHeader();
889 if(!esdHeader) {AliInfo("no esdHeader");return;}
890 //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
891 //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
897 Double_t vert[3] = {0};
898 vert[0] = vtxESD->GetXv();
899 vert[1] = vtxESD->GetYv();
900 vert[2] = vtxESD->GetZv();
901 Int_t mult = vtxESD->GetNContributors();
902 Float_t bz = esdEvent->GetMagneticField();
903 Int_t runNumber = esdEvent->GetRunNumber();
904 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
905 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
907 Int_t numberOfTracks=esdEvent->GetNumberOfTracks();
909 for (Int_t iTrack = 0; iTrack < numberOfTracks; iTrack++)
911 AliESDtrack *track = esdEvent->GetTrack(iTrack);
912 AliESDfriendTrack* friendTrack=NULL;
913 Int_t numberOfFriendTracks=0;
914 if (esdFriend) numberOfFriendTracks=esdFriend->GetNumberOfTracks();
915 if (esdFriend && iTrack<numberOfFriendTracks) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
917 if(track->Charge()==0) continue;
918 if(!esdTrackCuts->AcceptTrack(track)) continue;
919 if(!accCuts->AcceptTrack(track)) continue;
921 // downscale low-pT tracks
922 Double_t scalempt= TMath::Min(track->Pt(),10.);
923 Double_t downscaleF = gRandom->Rndm();
924 downscaleF *= fLowPtTrackDownscaligF;
925 if(TMath::Exp(2*scalempt)<downscaleF) continue;
926 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
930 // TPC constrained tracks
931 // InnerParams constrained tracks
933 // ITSout-InnerParams tracks
934 // chi2 distance between TPC constrained and TPC-ITS tracks
935 // chi2 distance between TPC refitted constrained and TPC-ITS tracks
936 // chi2 distance between ITSout and InnerParams tracks
939 Double_t x[3]; track->GetXYZ(x);
940 Double_t b[3]; AliTracker::GetBxByBz(x,b);
943 // Transform TPC inner params to track reference frame
945 Bool_t isOKtpcInner = kFALSE;
946 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
948 // transform to the track reference frame
949 isOKtpcInner = tpcInner->Rotate(track->GetAlpha());
950 isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
954 // Constrain TPC inner to vertex
955 // clone TPCinner has to be deleted
957 Bool_t isOKtpcInnerC = kFALSE;
958 AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
960 isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);
961 isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());
962 isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
966 // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex
967 // Clone track InnerParams has to be deleted
969 Bool_t isOKtrackInnerC = kFALSE;
970 AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));
972 isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);
973 isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());
974 isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
978 // calculate chi2 between vi and vj vectors
979 // with covi and covj covariance matrices
980 // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
982 TMatrixD deltaT(5,1), deltaTtrackC(5,1);
983 TMatrixD delta(1,5), deltatrackC(1,5);
984 TMatrixD covarM(5,5), covarMtrackC(5,5);
986 TMatrixD chi2trackC(1,1);
988 if(isOKtpcInnerC && isOKtrackInnerC)
990 for (Int_t ipar=0; ipar<5; ipar++) {
991 deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
992 delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
994 deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
995 deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
997 for (Int_t jpar=0; jpar<5; jpar++) {
998 Int_t index=track->GetIndex(ipar,jpar);
999 covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
1000 covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];
1004 // chi2 distance TPC constrained and TPC+ITS
1005 TMatrixD covarMInv = covarM.Invert();
1006 TMatrixD mat2 = covarMInv*deltaT;
1010 // chi2 distance TPC refitted constrained and TPC+ITS
1011 TMatrixD covarMInvtrackC = covarMtrackC.Invert();
1012 TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;
1013 chi2trackC = deltatrackC*mat2trackC;
1014 //chi2trackC.Print();
1019 // Propagate ITSout to TPC inner wall
1020 // and calculate chi2 distance to track (InnerParams)
1022 const Double_t kTPCRadius=85;
1023 const Double_t kStep=3;
1025 // clone track InnerParams has to be deleted
1026 Bool_t isOKtrackInnerC2 = kFALSE;
1027 AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));
1029 isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);
1032 Bool_t isOKouterITSc = kFALSE;
1033 AliExternalTrackParam *outerITSc = NULL;
1034 TMatrixD chi2OuterITS(1,1);
1036 if(esdFriend && !esdFriend->TestSkipBit())
1038 // propagate ITSout to TPC inner wall
1043 if (friendTrack->GetITSOut()) outerITSc = new AliExternalTrackParam(*(friendTrack->GetITSOut()));
1046 isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);
1047 isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());
1048 isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());
1051 // calculate chi2 between outerITS and innerParams
1052 // cov without z-coordinate at the moment
1054 TMatrixD deltaTouterITS(4,1);
1055 TMatrixD deltaouterITS(1,4);
1056 TMatrixD covarMouterITS(4,4);
1058 if(isOKtrackInnerC2 && isOKouterITSc) {
1061 for (Int_t ipar=0; ipar<5; ipar++) {
1063 deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
1064 deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
1068 for (Int_t jpar=0; jpar<5; jpar++) {
1069 Int_t index=outerITSc->GetIndex(ipar,jpar);
1070 if(ipar !=1 || jpar!=1) {
1071 covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];
1073 if(jpar!=1) kjpar++;
1075 if(ipar!=1) kipar++;
1078 // chi2 distance ITSout and InnerParams
1079 TMatrixD covarMInvouterITS = covarMouterITS.Invert();
1080 TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;
1081 chi2OuterITS = deltaouterITS*mat2outerITS;
1082 //chi2OuterITS.Print();
1091 TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;
1092 TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;
1093 Int_t mech=-1, mechTPC=-1, mechITS=-1;
1094 Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;
1095 Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;
1096 Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;
1097 Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;
1099 AliTrackReference *refTPCIn = NULL;
1100 AliTrackReference *refTPCOut = NULL;
1101 AliTrackReference *refITS = NULL;
1102 AliTrackReference *refTRD = NULL;
1103 AliTrackReference *refTOF = NULL;
1104 AliTrackReference *refEMCAL = NULL;
1105 AliTrackReference *refPHOS = NULL;
1106 Int_t nrefTPC=0, nrefTRD=0, nrefTOF=0, nrefITS=0, nrefEMCAL=0, nrefPHOS=0;
1108 Bool_t isOKtrackInnerC3 = kFALSE;
1109 AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
1110 if(mcEvent && stack)
1112 do //artificial loop (once) to make the continue statements jump out of the MC part
1114 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1118 Int_t label = TMath::Abs(track->GetLabel());
1119 if (label >= mcStackSize) continue;
1120 particle = stack->Particle(label);
1121 if (!particle) continue;
1122 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)
1124 particleMother = GetMother(particle,stack);
1125 mech = particle->GetUniqueID();
1126 isPrim = stack->IsPhysicalPrimary(label);
1127 isFromStrangess = IsFromStrangeness(label,stack);
1128 isFromConversion = IsFromConversion(label,stack);
1129 isFromMaterial = IsFromMaterial(label,stack);
1135 Int_t labelTPC = TMath::Abs(track->GetTPCLabel());
1136 if (labelTPC >= mcStackSize) continue;
1137 particleTPC = stack->Particle(labelTPC);
1138 if (!particleTPC) continue;
1139 if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)
1141 particleMotherTPC = GetMother(particleTPC,stack);
1142 mechTPC = particleTPC->GetUniqueID();
1143 isPrimTPC = stack->IsPhysicalPrimary(labelTPC);
1144 isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);
1145 isFromConversionTPC = IsFromConversion(labelTPC,stack);
1146 isFromMaterialTPC = IsFromMaterial(labelTPC,stack);
1150 // store first track reference
1154 TClonesArray *trefs=0;
1155 Int_t status = mcEvent->GetParticleAndTR(TMath::Abs(labelTPC), part, trefs);
1157 if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.)
1159 Int_t nTrackRef = trefs->GetEntries();
1160 //printf("nTrackRef %d \n",nTrackRef);
1163 for (Int_t iref = 0; iref < nTrackRef; iref++)
1165 AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);
1167 // Ref. in the middle ITS
1168 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kITS)
1171 if(!refITS && countITS==2) {
1173 //printf("refITS %p \n",refITS);
1179 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC)
1185 //printf("refTPCIn %p \n",refTPCIn);
1190 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD)
1198 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF)
1206 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL)
1214 // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS)
1223 // transform inner params to TrackRef
1225 if(refTPCIn && trackInnerC3)
1227 Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());
1228 isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);
1229 isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);
1236 Int_t labelITS = TMath::Abs(track->GetITSLabel());
1237 if (labelITS >= mcStackSize) continue;
1238 particleITS = stack->Particle(labelITS);
1239 if (!particleITS) continue;
1240 if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)
1242 particleMotherITS = GetMother(particleITS,stack);
1243 mechITS = particleITS->GetUniqueID();
1244 isPrimITS = stack->IsPhysicalPrimary(labelITS);
1245 isFromStrangessITS = IsFromStrangeness(labelITS,stack);
1246 isFromConversionITS = IsFromConversion(labelITS,stack);
1247 isFromMaterialITS = IsFromMaterial(labelITS,stack);
1254 Bool_t dumpToTree=kFALSE;
1256 if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE;
1257 //if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1258 if(isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1259 if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE;
1260 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1263 // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks
1264 // Done only in case no MC info
1267 track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]);
1268 Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10;
1269 Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin));
1270 Bool_t keepPileUp=gRandom->Rndm()<0.05;
1271 if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){
1276 //init dummy objects
1277 static AliESDVertex dummyvtxESD;
1280 // dummyvtxESD=new AliESDVertex();
1281 // //dummyvtxESD->SetNContributors(1);
1282 // //UShort_t pindices[1]; pindices[0]=0;
1283 // //dummyvtxESD->SetIndices(1,pindices);
1284 // //dummyvtxESD->SetNContributors(0);
1286 static AliExternalTrackParam dummyexternaltrackparam;
1287 //if (!dummyexternaltrackparam) dummyexternaltrackparam=new AliExternalTrackParam();
1288 static AliTrackReference dummytrackreference;
1289 //if (!dummytrackreference) dummytrackreference=new AliTrackReference();
1290 static TParticle dummyparticle;
1291 //if (!dummyparticle) dummyparticle=new TParticle();
1293 //assign the dummy objects if needed
1294 if (!track) {track=fDummyTrack;}
1295 AliESDfriendTrack *friendTrackStore=friendTrack; // store friend track for later processing
1296 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
1297 friendTrackStore = (gRandom->Rndm()<1./fFriendDownscaling)? friendTrack:0;
1299 if (fFriendDownscaling<=0){
1300 if (((*fTreeSRedirector)<<"highPt").GetTree()){
1301 TTree * tree = ((*fTreeSRedirector)<<"highPt").GetTree();
1303 Double_t sizeAll=tree->GetZipBytes();
1304 TBranch * br= tree->GetBranch("friendTrack.fPoints");
1305 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
1306 br= tree->GetBranch("friendTrack.fCalibContainer");
1307 if (br) sizeFriend+=br->GetZipBytes();
1308 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) friendTrackStore=0;
1314 // if (!friendTrackStore && fFriendDownscaling<=1) {friendTrack=fDummyFriendTrack;}
1315 if (!vtxESD) {vtxESD=&dummyvtxESD;}
1318 if (!refTPCIn) {refTPCIn=&dummytrackreference;}
1319 if (!refITS) {refITS=&dummytrackreference;}
1320 if (!particle) {particle=&dummyparticle;}
1321 if (!particleMother) {particleMother=&dummyparticle;}
1322 if (!particleTPC) {particleTPC=&dummyparticle;}
1323 if (!particleMotherTPC) {particleMotherTPC=&dummyparticle;}
1324 if (!particleITS) {particleITS=&dummyparticle;}
1325 if (!particleMotherITS) {particleMotherITS=&dummyparticle;}
1328 Int_t ntracks = esdEvent->GetNumberOfTracks();
1330 // Global event id calculation using orbitID, bunchCrossingID and periodID
1331 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1332 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1333 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1334 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1337 FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0));
1339 if(fTreeSRedirector && dumpToTree && fFillTree) {
1340 (*fTreeSRedirector)<<"highPt"<<
1342 "fileName.="<<&fCurrentFileName<< // name of the chunk file (hopefully full)
1343 "runNumber="<<runNumber<< // runNumber
1344 "evtTimeStamp="<<evtTimeStamp<< // time stamp of event (in seconds)
1345 "evtNumberInFile="<<evtNumberInFile<< // event number
1346 "triggerClass="<<&triggerClass<< // trigger class as a string
1347 "Bz="<<bz<< // solenoid magnetic field in the z direction (in kGaus)
1348 "vtxESD.="<<vtxESD<< // vertexer ESD tracks (can be biased by TPC pileup tracks)
1349 "IRtot="<<ir1<< // interaction record (trigger) counters - coutner 1
1350 "IRint2="<<ir2<< // interaction record (trigger) coutners - counter 2
1351 "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
1352 "ntracks="<<ntracks<< // number of the esd tracks (to take into account the pileup in the TPC)
1353 // important variables for the pile-up studies
1354 "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
1355 "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
1356 "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
1357 "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
1358 "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
1359 "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
1361 "esdTrack.="<<track<< // esdTrack as used in the physical analysis
1362 // "friendTrack.="<<friendTrack<< // esdFriendTrack associated to the esdTrack
1363 "friendTrack.="<<friendTrackStore<< // esdFriendTrack associated to the esdTrack
1364 "extTPCInnerC.="<<tpcInnerC<< // ???
1365 "extInnerParamC.="<<trackInnerC<< // ???
1366 "extInnerParam.="<<trackInnerC2<< // ???
1367 "extOuterITS.="<<outerITSc<< // ???
1368 "extInnerParamRef.="<<trackInnerC3<< // ???
1369 "chi2TPCInnerC="<<chi2(0,0)<< // chi2 of tracks ???
1370 "chi2InnerC="<<chi2trackC(0,0)<< // chi2s of tracks TPCinner to the combined
1371 "chi2OuterITS="<<chi2OuterITS(0,0)<< // chi2s of tracks TPC at inner wall to the ITSout
1372 "centralityF="<<centralityF;
1374 static AliTrackReference refDummy;
1375 if (!refITS) refITS = &refDummy;
1376 if (!refTRD) refTRD = &refDummy;
1377 if (!refTOF) refTOF = &refDummy;
1378 if (!refEMCAL) refEMCAL = &refDummy;
1379 if (!refPHOS) refPHOS = &refDummy;
1380 (*fTreeSRedirector)<<"highPt"<<
1381 "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
1382 "nrefITS="<<nrefITS<< // number of track references in the ITS
1383 "nrefTPC="<<nrefTPC<< // number of track references in the TPC
1384 "nrefTRD="<<nrefTRD<< // number of track references in the TRD
1385 "nrefTOF="<<nrefTOF<< // number of track references in the TOF
1386 "nrefEMCAL="<<nrefEMCAL<< // number of track references in the TOF
1387 "nrefPHOS="<<nrefPHOS<< // number of track references in the TOF
1388 "refTPCIn.="<<refTPCIn<<
1389 "refTPCOut.="<<refTPCOut<<
1390 "refITS.="<<refITS<<
1391 "refTRD.="<<refTRD<<
1392 "refTOF.="<<refTOF<<
1393 "refEMCAL.="<<refEMCAL<<
1394 "refPHOS.="<<refPHOS<<
1395 "particle.="<<particle<<
1396 "particleMother.="<<particleMother<<
1399 "isFromStrangess="<<isFromStrangess<<
1400 "isFromConversion="<<isFromConversion<<
1401 "isFromMaterial="<<isFromMaterial<<
1402 "particleTPC.="<<particleTPC<<
1403 "particleMotherTPC.="<<particleMotherTPC<<
1404 "mechTPC="<<mechTPC<<
1405 "isPrimTPC="<<isPrimTPC<<
1406 "isFromStrangessTPC="<<isFromStrangessTPC<<
1407 "isFromConversionTPC="<<isFromConversionTPC<<
1408 "isFromMaterialTPC="<<isFromMaterialTPC<<
1409 "particleITS.="<<particleITS<<
1410 "particleMotherITS.="<<particleMotherITS<<
1411 "mechITS="<<mechITS<<
1412 "isPrimITS="<<isPrimITS<<
1413 "isFromStrangessITS="<<isFromStrangessITS<<
1414 "isFromConversionITS="<<isFromConversionITS<<
1415 "isFromMaterialITS="<<isFromMaterialITS;
1417 //finish writing the entry
1418 AliInfo("writing tree highPt");
1419 (*fTreeSRedirector)<<"highPt"<<"\n";
1421 AliSysInfo::AddStamp("filteringTask",iTrack,numberOfTracks,numberOfFriendTracks,(friendTrackStore)?0:1);
1424 delete trackInnerC2;
1426 delete trackInnerC3;
1432 //_____________________________________________________________________________
1433 void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1436 // Fill tree for efficiency studies MC only
1437 AliInfo("we start!");
1439 AliDebug(AliLog::kError, "mcEvent not available");
1443 // get selection cuts
1444 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1445 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1446 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1448 if(!evtCuts || !accCuts || !esdTrackCuts) {
1449 AliDebug(AliLog::kError, "cuts not available");
1453 // trigger selection
1454 Bool_t isEventTriggered = kTRUE;
1455 AliPhysicsSelection *physicsSelection = NULL;
1456 AliTriggerAnalysis* triggerAnalysis = NULL;
1459 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1463 if(evtCuts->IsTriggerRequired())
1466 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1467 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1468 if(!physicsSelection) return;
1470 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1471 // set trigger (V0AND)
1472 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1473 if(!triggerAnalysis) return;
1474 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1478 // centrality determination
1479 Float_t centralityF = -1;
1480 AliCentrality *esdCentrality = esdEvent->GetCentrality();
1481 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1483 // use MC information
1484 AliHeader* header = 0;
1485 AliGenEventHeader* genHeader = 0;
1486 AliStack* stack = 0;
1487 Int_t mcStackSize=0;
1490 Int_t multMCTrueTracks = 0;
1493 AliDebug(AliLog::kError, "mcEvent not available");
1496 // get MC event header
1497 header = mcEvent->Header();
1499 AliDebug(AliLog::kError, "Header not available");
1502 // MC particle stack
1503 stack = mcEvent->Stack();
1505 AliDebug(AliLog::kError, "Stack not available");
1508 mcStackSize=stack->GetNtrack();
1511 genHeader = header->GenEventHeader();
1513 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1516 genHeader->PrimaryVertex(vtxMC);
1518 // multipliticy of all MC primary tracks
1519 // in Zv, pt and eta ranges)
1520 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1523 // get reconstructed vertex
1524 //const AliESDVertex* vtxESD = 0;
1525 AliESDVertex* vtxESD = 0;
1526 if(GetAnalysisMode() == kTPCAnalysisMode) {
1527 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1529 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1530 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1538 // Vertex info comparison and track multiplicity
1540 AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
1541 AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
1542 Int_t contSPD = vertexSPD->GetNContributors();
1543 Int_t contTPC = vertexTPC->GetNContributors();
1544 TVectorD vertexPosTPC(3), vertexPosSPD(3);
1545 vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
1546 vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
1549 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
1550 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1551 if(!track) continue;
1552 if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
1553 if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
1556 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1557 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1558 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1560 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1563 if(isEventOK && isEventTriggered)
1565 //if(!stack) return;
1570 TParticle *particle=NULL;
1571 TParticle *particleMother=NULL;
1575 Double_t vert[3] = {0};
1576 vert[0] = vtxESD->GetXv();
1577 vert[1] = vtxESD->GetYv();
1578 vert[2] = vtxESD->GetZv();
1579 Int_t mult = vtxESD->GetNContributors();
1580 Double_t bz = esdEvent->GetMagneticField();
1581 Double_t runNumber = esdEvent->GetRunNumber();
1582 Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1583 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1584 // loop over MC stack
1585 for (Int_t iMc = 0; iMc < mcStackSize; ++iMc)
1587 particle = stack->Particle(iMc);
1591 // only charged particles
1592 if(!particle->GetPDG()) continue;
1593 Double_t charge = particle->GetPDG()->Charge()/3.;
1594 if (TMath::Abs(charge) < 0.001)
1597 // only primary particles
1598 Bool_t prim = stack->IsPhysicalPrimary(iMc);
1601 // downscale low-pT particles
1602 Double_t scalempt= TMath::Min(particle->Pt(),10.);
1603 Double_t downscaleF = gRandom->Rndm();
1604 downscaleF *= fLowPtTrackDownscaligF;
1605 if(TMath::Exp(2*scalempt)<downscaleF) continue;
1607 // is particle in acceptance
1608 if(!accCuts->AcceptTrack(particle)) continue;
1610 // check if particle reconstructed
1611 Bool_t isRec = kFALSE;
1612 Int_t trackIndex = -1;
1613 Int_t trackLoopIndex = -1;
1614 Int_t isESDtrackCut= 0;
1615 Int_t isAccCuts = 0;
1616 Int_t nRec = 0; // how many times reconstructed
1617 Int_t nFakes = 0; // how many times reconstructed as a fake track
1618 AliESDtrack *recTrack = NULL;
1620 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1622 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1623 if(!track) continue;
1624 if(track->Charge()==0) continue;
1626 Int_t label = TMath::Abs(track->GetLabel());
1627 if (label >= mcStackSize) continue;
1629 Bool_t isAcc=esdTrackCuts->AcceptTrack(track);
1630 if (isAcc) isESDtrackCut=1;
1631 if (accCuts->AcceptTrack(track)) isAccCuts=1;
1633 trackIndex = iTrack;
1636 if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track
1637 if (!isAcc) continue;
1638 trackLoopIndex = iTrack;
1640 recTrack = esdEvent->GetTrack(trackIndex);
1642 if(track->GetLabel()<0) nFakes++;
1648 // Store information in the output tree
1649 if (trackLoopIndex>-1) {
1650 recTrack = esdEvent->GetTrack(trackLoopIndex);
1651 } else if (trackIndex >-1) {
1652 recTrack = esdEvent->GetTrack(trackIndex);
1654 recTrack = fDummyTrack;
1657 particleMother = GetMother(particle,stack);
1658 mech = particle->GetUniqueID();
1660 //MC particle track length
1661 Double_t tpcTrackLength = 0.;
1662 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
1665 tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
1670 if(fTreeSRedirector && fFillTree) {
1671 (*fTreeSRedirector)<<"MCEffTree"<<
1672 "fileName.="<<&fCurrentFileName<<
1673 "triggerClass.="<<&triggerClass<<
1674 "runNumber="<<runNumber<<
1675 "evtTimeStamp="<<evtTimeStamp<<
1676 "evtNumberInFile="<<evtNumberInFile<< //
1677 "Bz="<<bz<< // magnetic field
1678 "vtxESD.="<<vtxESD<< // vertex info
1680 "mult="<<mult<< // primary vertex 9whatewe found) multiplicity
1681 "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
1682 // important variables for the pile-up studies
1683 "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
1684 "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
1685 "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
1686 "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
1687 "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
1688 "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
1691 "isAcc0="<<isESDtrackCut<< // track accepted by ESD track cuts
1692 "isAcc1="<<isAccCuts<< // track accepted by acceptance cuts flag
1693 "esdTrack.="<<recTrack<< // reconstructed track (only the longest from the loopers)
1694 "isRec="<<isRec<< // track was reconstructed
1695 "tpcTrackLength="<<tpcTrackLength<< // track length in the TPC r projection
1696 "particle.="<<particle<< // particle properties
1697 "particleMother.="<<particleMother<< // particle mother
1698 "mech="<<mech<< // production mechanizm
1699 "nRec="<<nRec<< // how many times reconstruted
1700 "nFakes="<<nFakes<< // how many times reconstructed as a fake track
1704 //if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
1710 //_____________________________________________________________________________
1711 Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {
1713 // check if particle is Z > 1
1715 if (track->GetTPCNcls() < 60) return kFALSE;
1716 Double_t mom = track->GetInnerParam()->GetP();
1717 if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
1718 Float_t dca[2], bCov[3];
1719 track->GetImpactParameters(dca,bCov);
1722 Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
1724 if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
1729 //_____________________________________________________________________________
1730 void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
1733 // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates
1736 AliDebug(AliLog::kError, "esdEvent not available");
1740 // get selection cuts
1741 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1742 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1743 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1745 if(!evtCuts || !accCuts || !esdTrackCuts) {
1746 AliDebug(AliLog::kError, "cuts not available");
1750 // trigger selection
1751 Bool_t isEventTriggered = kTRUE;
1752 AliPhysicsSelection *physicsSelection = NULL;
1753 AliTriggerAnalysis* triggerAnalysis = NULL;
1756 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1759 if(evtCuts->IsTriggerRequired())
1762 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1764 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1765 if(!physicsSelection) return;
1766 //SetPhysicsTriggerSelection(physicsSelection);
1768 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1769 // set trigger (V0AND)
1770 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1771 if(!triggerAnalysis) return;
1772 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1776 // centrality determination
1777 Float_t centralityF = -1;
1778 AliCentrality *esdCentrality = esdEvent->GetCentrality();
1779 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1782 // get reconstructed vertex
1783 //const AliESDVertex* vtxESD = 0;
1784 AliESDVertex* vtxESD = 0;
1785 if(GetAnalysisMode() == kTPCAnalysisMode) {
1786 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1788 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1789 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1797 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1798 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1799 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1802 if(isEventOK && isEventTriggered) {
1804 // Dump the pt downscaled V0 into the tree
1806 Int_t ntracks = esdEvent->GetNumberOfTracks();
1807 Int_t nV0s = esdEvent->GetNumberOfV0s();
1808 Int_t run = esdEvent->GetRunNumber();
1809 Int_t time= esdEvent->GetTimeStamp();
1810 Int_t evNr=esdEvent->GetEventNumberInFile();
1811 Double_t bz = esdEvent->GetMagneticField();
1813 // Global event id calculation using orbitID, bunchCrossingID and periodID
1814 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1815 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1816 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1817 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1820 for (Int_t iv0=0; iv0<nV0s; iv0++){
1822 AliESDv0 * v0 = esdEvent->GetV0(iv0);
1824 AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
1825 AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
1826 if (!track0) continue;
1827 if (!track1) continue;
1828 AliESDfriendTrack* friendTrack0=NULL;
1829 AliESDfriendTrack* friendTrack1=NULL;
1831 if (!esdFriend->TestSkipBit()){
1832 friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
1833 friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
1836 if (track0->GetSign()<0) {
1837 track1 = esdEvent->GetTrack(v0->GetIndex(0));
1838 track0 = esdEvent->GetTrack(v0->GetIndex(1));
1842 AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
1843 AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
1844 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
1845 if (gRandom->Rndm()>1./fFriendDownscaling){
1846 friendTrackStore0 = 0;
1847 friendTrackStore1 = 0;
1850 if (fFriendDownscaling<=0){
1851 if (((*fTreeSRedirector)<<"V0s").GetTree()){
1852 TTree * tree = ((*fTreeSRedirector)<<"V0s").GetTree();
1854 Double_t sizeAll=tree->GetZipBytes();
1855 TBranch * br= tree->GetBranch("friendTrack0.fPoints");
1856 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
1857 br= tree->GetBranch("friendTrack0.fCalibContainer");
1858 if (br) sizeFriend+=br->GetZipBytes();
1859 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
1860 friendTrackStore0=0;
1861 friendTrackStore1=0;
1868 Bool_t isDownscaled = IsV0Downscaled(v0);
1869 if (isDownscaled) continue;
1870 AliKFParticle kfparticle; //
1871 Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
1872 if (type==0) continue;
1873 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1875 if(!fFillTree) return;
1876 if(!fTreeSRedirector) return;
1877 (*fTreeSRedirector)<<"V0s"<<
1878 "gid="<<gid<< // global id of event
1879 "isDownscaled="<<isDownscaled<< //
1880 "triggerClass="<<&triggerClass<< // trigger
1882 "fileName.="<<&fCurrentFileName<< // full path - file name with ESD
1883 "runNumber="<<run<< //
1884 "evtTimeStamp="<<time<< // time stamp of event in secons
1885 "evtNumberInFile="<<evNr<< //
1886 "type="<<type<< // type of V0-
1887 "ntracks="<<ntracks<<
1889 "kf.="<<&kfparticle<<
1890 "track0.="<<track0<< // track
1891 "track1.="<<track1<<
1892 "friendTrack0.="<<friendTrackStore0<<
1893 "friendTrack1.="<<friendTrackStore1<<
1894 "centralityF="<<centralityF<<
1900 //_____________________________________________________________________________
1901 void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
1904 // Select real events with large TPC dEdx signal
1907 AliDebug(AliLog::kError, "esdEvent not available");
1911 // get selection cuts
1912 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1913 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1914 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1916 if(!evtCuts || !accCuts || !esdTrackCuts) {
1917 AliDebug(AliLog::kError, "cuts not available");
1923 Bool_t isEventTriggered = kTRUE;
1924 AliPhysicsSelection *physicsSelection = NULL;
1925 AliTriggerAnalysis* triggerAnalysis = NULL;
1928 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1930 if(evtCuts->IsTriggerRequired())
1933 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1935 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1936 if(!physicsSelection) return;
1938 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1939 // set trigger (V0AND)
1940 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1941 if(!triggerAnalysis) return;
1942 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1946 // get reconstructed vertex
1947 AliESDVertex* vtxESD = 0;
1948 if(GetAnalysisMode() == kTPCAnalysisMode) {
1949 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1951 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1952 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1959 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1960 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1961 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1965 if(isEventOK && isEventTriggered)
1967 Double_t vert[3] = {0};
1968 vert[0] = vtxESD->GetXv();
1969 vert[1] = vtxESD->GetYv();
1970 vert[2] = vtxESD->GetZv();
1971 Int_t mult = vtxESD->GetNContributors();
1972 Double_t bz = esdEvent->GetMagneticField();
1973 Double_t runNumber = esdEvent->GetRunNumber();
1974 Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1975 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
1977 // Global event id calculation using orbitID, bunchCrossingID and periodID
1978 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1979 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1980 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1981 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1984 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1986 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1987 if(!track) continue;
1988 AliESDfriendTrack* friendTrack=NULL;
1989 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
1990 if(track->Charge()==0) continue;
1991 if(!esdTrackCuts->AcceptTrack(track)) continue;
1992 if(!accCuts->AcceptTrack(track)) continue;
1994 if(!IsHighDeDxParticle(track)) continue;
1995 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1997 if(!fFillTree) return;
1998 if(!fTreeSRedirector) return;
1999 (*fTreeSRedirector)<<"dEdx"<< // high dEdx tree
2000 "gid="<<gid<< // global id
2001 "fileName.="<<&fCurrentFileName<< // file name
2002 "runNumber="<<runNumber<<
2003 "evtTimeStamp="<<evtTimeStamp<<
2004 "evtNumberInFile="<<evtNumberInFile<<
2005 "triggerClass="<<&triggerClass<< // trigger
2007 "vtxESD.="<<vtxESD<< //
2009 "esdTrack.="<<track<<
2010 "friendTrack.="<<friendTrack<<
2016 //_____________________________________________________________________________
2017 Int_t AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
2020 // Create KF particle in case the V0 fullfill selection criteria
2022 // Selection criteria
2026 // 4. rough mass cut
2027 // 5. Normalized pointing angle cut
2029 const Double_t cutMass=0.2;
2030 const Double_t kSigmaDCACut=3;
2032 // 0.) algo cut - accept only on the fly
2034 if (v0->GetOnFlyStatus() ==kFALSE) return 0;
2038 AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
2039 AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
2041 TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
2042 TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
2043 TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
2045 if (TMath::Abs(track0->GetTgl())>1) return 0;
2046 if (TMath::Abs(track1->GetTgl())>1) return 0;
2047 if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
2048 if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
2049 Float_t pos0[2]={0}, cov0[3]={0};
2050 Float_t pos1[2]={0}, cov1[3]={0};
2051 track0->GetImpactParameters(pos0,cov0);
2052 track0->GetImpactParameters(pos1,cov1);
2054 if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
2055 if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
2060 Double_t chi2KF = v0->GetKFInfo(2,2,2);
2061 if (chi2KF>25) return 0;
2063 // 4.) Rough mass cut - 0.200 GeV
2065 static Double_t masses[2]={-1};
2067 masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
2068 masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
2070 Double_t mass00= v0->GetEffMass(0,0);
2071 Double_t mass22= v0->GetEffMass(2,2);
2072 Double_t mass42= v0->GetEffMass(4,2);
2073 Double_t mass24= v0->GetEffMass(2,4);
2074 Bool_t massOK=kFALSE;
2079 if (TMath::Abs(mass00-0)<cutMass) {
2080 massOK=kTRUE; type+=1;
2081 if (TMath::Abs(mass00-0)<dmass) {
2083 dmass=TMath::Abs(mass00-0);
2087 if (TMath::Abs(mass24-masses[1])<cutMass) {
2088 massOK=kTRUE; type+=2;
2089 if (TMath::Abs(mass24-masses[1])<dmass){
2090 dmass = TMath::Abs(mass24-masses[1]);
2095 if (TMath::Abs(mass42-masses[1])<cutMass) {
2096 massOK=kTRUE; type+=4;
2097 if (TMath::Abs(mass42-masses[1])<dmass){
2098 dmass = TMath::Abs(mass42-masses[1]);
2103 if (TMath::Abs(mass22-masses[0])<cutMass) {
2104 massOK=kTRUE; type+=8;
2105 if (TMath::Abs(mass22-masses[0])<dmass){
2106 dmass = TMath::Abs(mass22-masses[0]);
2111 if (type==0) return 0;
2113 const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
2114 const AliExternalTrackParam *paramP = v0->GetParamP();
2115 const AliExternalTrackParam *paramN = v0->GetParamN();
2116 if (paramP->GetSign()<0){
2117 paramP=v0->GetParamP();
2118 paramN=v0->GetParamN();
2120 //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
2121 //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
2123 AliKFParticle kfp1( *paramP, spdg[p1] );
2124 AliKFParticle kfp2( *paramN, -1 *spdg[p2] );
2132 Double_t errPhi = V0KF.GetErrPhi();
2133 Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
2134 if (pointAngle/errPhi>10) return 0;
2139 //_____________________________________________________________________________
2140 Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)
2143 // Downscale randomly low pt V0
2146 Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
2147 Double_t scalempt= TMath::Min(maxPt,10.);
2148 Double_t downscaleF = gRandom->Rndm();
2149 downscaleF *= fLowPtV0DownscaligF;
2151 // Special treatment of the gamma conversion pt spectra is softer -
2152 Double_t mass00= v0->GetEffMass(0,0);
2153 const Double_t cutMass=0.2;
2154 if (TMath::Abs(mass00-0)<cutMass){
2155 downscaleF/=10.; // 10 times smaller downscaling for the gamma concersion candidate
2157 //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
2158 if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;
2162 TH1F his1("his1","his1",100,0,10);
2163 TH1F his2("his2","his2",100,0,10);
2164 {for (Int_t i=0; i<10000; i++){
2165 Double_t rnd=gRandom->Exp(1);
2166 Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
2168 if (!isDownscaled) his2->Fill(rnd);
2177 //_____________________________________________________________________________
2178 Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
2180 // Constrain TPC inner params constrained
2182 if(!tpcInnerC) return kFALSE;
2183 if(!vtx) return kFALSE;
2185 Double_t dz[2],cov[3];
2186 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2187 //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2188 //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2189 if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
2192 Double_t covar[6]; vtx->GetCovMatrix(covar);
2193 Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
2194 Double_t c[3]={covar[2],0.,covar[5]};
2195 Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
2196 if (chi2C>kVeryBig) return kFALSE;
2198 if(!tpcInnerC->Update(p,c)) return kFALSE;
2203 //_____________________________________________________________________________
2204 Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
2206 // Constrain TPC inner params constrained
2208 if(!trackInnerC) return kFALSE;
2209 if(!vtx) return kFALSE;
2211 const Double_t kRadius = 2.8;
2212 const Double_t kMaxStep = 1.0;
2214 Double_t dz[2],cov[3];
2216 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2217 //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2218 //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2220 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
2221 if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
2224 Double_t covar[6]; vtx->GetCovMatrix(covar);
2225 Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
2226 Double_t c[3]={covar[2],0.,covar[5]};
2227 Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
2228 if (chi2C>kVeryBig) return kFALSE;
2230 if(!trackInnerC->Update(p,c)) return kFALSE;
2236 //_____________________________________________________________________________
2237 TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack)
2239 if(!particle) return NULL;
2240 if(!stack) return NULL;
2242 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2243 TParticle* mother = NULL;
2244 Int_t mcStackSize=stack->GetNtrack();
2245 if (motherLabel>=mcStackSize) return NULL;
2246 mother = stack->Particle(motherLabel);
2251 //_____________________________________________________________________________
2252 Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(Int_t label, AliStack *const stack)
2254 Bool_t isFromConversion = kFALSE;
2257 Int_t mcStackSize=stack->GetNtrack();
2258 if (label>=mcStackSize) return kFALSE;
2259 TParticle* particle = stack->Particle(label);
2260 if (!particle) return kFALSE;
2262 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2264 Int_t mech = particle->GetUniqueID(); // production mechanism
2265 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2267 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2268 if (motherLabel>=mcStackSize) return kFALSE;
2269 TParticle* mother = stack->Particle(motherLabel);
2271 Int_t motherPdg = mother->GetPdgCode();
2273 if(!isPrim && mech==5 && motherPdg==kGamma) {
2274 isFromConversion=kTRUE;
2280 return isFromConversion;
2283 //_____________________________________________________________________________
2284 Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(Int_t label, AliStack *const stack)
2286 Bool_t isFromMaterial = kFALSE;
2289 Int_t mcStackSize=stack->GetNtrack();
2290 if (label>=mcStackSize) return kFALSE;
2291 TParticle* particle = stack->Particle(label);
2292 if (!particle) return kFALSE;
2294 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2296 Int_t mech = particle->GetUniqueID(); // production mechanism
2297 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2299 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2300 if (motherLabel>=mcStackSize) return kFALSE;
2301 TParticle* mother = stack->Particle(motherLabel);
2303 if(!isPrim && mech==13) {
2304 isFromMaterial=kTRUE;
2310 return isFromMaterial;
2313 //_____________________________________________________________________________
2314 Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(Int_t label, AliStack *const stack)
2316 Bool_t isFromStrangeness = kFALSE;
2319 Int_t mcStackSize=stack->GetNtrack();
2320 if (label>=mcStackSize) return kFALSE;
2321 TParticle* particle = stack->Particle(label);
2322 if (!particle) return kFALSE;
2324 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2326 Int_t mech = particle->GetUniqueID(); // production mechanism
2327 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2329 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2330 if (motherLabel>=mcStackSize) return kFALSE;
2331 TParticle* mother = stack->Particle(motherLabel);
2333 Int_t motherPdg = mother->GetPdgCode();
2335 // K+-, lambda, antilambda, K0s decays
2336 if(!isPrim && mech==4 &&
2337 (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
2339 isFromStrangeness = kTRUE;
2345 return isFromStrangeness;
2349 //_____________________________________________________________________________
2350 void AliAnalysisTaskFilteredTree::FinishTaskOutput()
2353 // Called one at the end
2354 // locally on working node
2356 Bool_t deleteTrees=kTRUE;
2357 if ((AliAnalysisManager::GetAnalysisManager()))
2359 if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
2360 AliAnalysisManager::kProofAnalysis)
2363 if (deleteTrees) delete fTreeSRedirector;
2364 fTreeSRedirector=NULL;
2367 //_____________________________________________________________________________
2368 void AliAnalysisTaskFilteredTree::Terminate(Option_t *)
2371 // Called one at the end
2375 //_____________________________________________________________________________
2376 Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
2379 // calculate mc event true track multiplicity
2381 if(!mcEvent) return 0;
2383 AliStack* stack = 0;
2386 // MC particle stack
2387 stack = mcEvent->Stack();
2388 if (!stack) return 0;
2391 //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
2394 Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
2395 if(!isEventOK) return 0;
2397 Int_t nPart = stack->GetNtrack();
2398 for (Int_t iMc = 0; iMc < nPart; ++iMc)
2400 TParticle* particle = stack->Particle(iMc);
2404 // only charged particles
2405 if(!particle->GetPDG()) continue;
2406 Double_t charge = particle->GetPDG()->Charge()/3.;
2407 if (TMath::Abs(charge) < 0.001)
2411 Bool_t prim = stack->IsPhysicalPrimary(iMc);
2414 // checked accepted without pt cut
2415 //if(accCuts->AcceptTrack(particle))
2416 if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
2425 //_____________________________________________________________________________
2426 void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC)
2429 // Fill pT relative resolution histograms for
2430 // TPC only, TPC only constrained to vertex and TPC+ITS tracking
2433 if(!ptpcInnerC) return;
2435 const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();
2436 if(!innerParam) return;
2439 ptrack->GetImpactParameters(dxy,dz);
2441 // TPC+ITS primary tracks
2442 if( abs(ptrack->Eta())<0.8 &&
2443 ptrack->GetTPCClusterInfo(3,1)>120 &&
2444 ptrack->IsOn(0x40) &&
2445 ptrack->GetTPCclusters(0)>0.0 &&
2446 ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2447 //abs(innerParam->GetX())>0.0 &&
2448 //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2449 //abs(innerParam->GetTgl())<0.85 &&
2450 ptrack->IsOn(0x0004) &&
2451 ptrack->GetNcls(0)>0 &&
2452 ptrack->GetITSchi2()>0 &&
2453 sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 &&
2454 sqrt(chi2TPCInnerC)<6 &&
2455 (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) &&
2457 abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )
2459 fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2460 fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2461 fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2464 // TPC primary tracks
2465 // and TPC constrained primary tracks
2467 AliExternalTrackParam *ptpcInner = (AliExternalTrackParam *) ptrack->GetTPCInnerParam();
2468 if(!ptpcInner) return;
2471 Float_t dxyTPC, dzTPC;
2472 ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);
2474 if( abs(ptrack->Eta())<0.8 &&
2475 ptrack->GetTPCClusterInfo(3,1)>120 &&
2476 ptrack->IsOn(0x40)&&
2477 ptrack->GetTPCclusters(0)>0.0 &&
2478 ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2479 //abs(innerParam->GetX())>0.0 &&
2480 //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2481 //abs(innerParam->GetTgl())<0.85 &&
2486 fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2487 fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2488 fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2490 // TPC constrained to vertex
2491 fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2492 fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2493 fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2498 void AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend){
2500 // Process ITS standalone tracks find match with closest TPC(or combined tracks) tracks
2501 // marian.ivanov@cern.ch
2502 // 0.) Init variables
2503 // 1.) GetTrack parameters at TPC inner wall
2504 // 2.) Match closest TPC track (STANDALONE/global) - chi2 match criteria
2506 // Logic to be used in reco:
2507 // 1.) Find matching ITSalone->TPCalone
2508 // 2.) if (!TPCalone.FindClose(TPCother)) TPCalone.Addopt(ITSalone)
2509 // 3.) ff ((ITSalone.FindClose(Global)==0) CreateGlobaltrack
2510 const Double_t radiusMatch=84.; // redius to propagate
2512 const Double_t dFastPhiCut=0.2; // 6 sigma (200 MeV) fast angular cut
2513 const Double_t dFastThetaCut=0.12; // 6 sigma (200 MeV) fast angular cut
2514 const Double_t dFastPosPhiCut=0.06; // 6 sigma (200 MeV) fast angular cut
2515 const Double_t dFastZCut=6; // 6 sigma (200 MeV) fast z difference cut
2516 const Double_t dFastPtCut=2.; // 6 sigma (200 MeV) fast 1/pt cut
2517 const Double_t chi2Cut=100; // chi2 matching cut
2519 if (!esdFriend) return; // not ITS standalone track
2520 if (esdFriend->TestSkipBit()) return; // friends tracks not stored
2521 Int_t ntracks=esdEvent->GetNumberOfTracks();
2522 Float_t bz = esdEvent->GetMagneticField();
2524 // 0.) Get parameters in reference radius TPC Inner wall
2527 TMatrixD vecPosR0(ntracks,6); // possition and momentum estimate at reference radius
2528 TMatrixD vecMomR0(ntracks,6); //
2529 TMatrixD vecPosR1(ntracks,6); // possition and momentum estimate at reference radius TPC track
2530 TMatrixD vecMomR1(ntracks,6); //
2531 Double_t xyz[3], pxyz[3]; //
2532 for (Int_t iTrack=0; iTrack<ntracks; iTrack++){
2533 AliESDtrack *track = esdEvent->GetTrack(iTrack);
2534 if(!track) continue;
2535 if (track->GetInnerParam()){
2536 const AliExternalTrackParam *trackTPC=track->GetInnerParam();
2537 trackTPC->GetXYZAt(radiusMatch,bz,xyz);
2538 trackTPC->GetPxPyPzAt(radiusMatch,bz,pxyz);
2539 for (Int_t i=0; i<3; i++){
2540 vecPosR1(iTrack,i)=xyz[i];
2541 vecMomR1(iTrack,i)=pxyz[i];
2543 vecPosR1(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); // phi pos angle
2544 vecMomR1(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); // phi mom angle
2545 vecMomR1(iTrack,4)= trackTPC->GetSigned1Pt();;
2546 vecMomR1(iTrack,5)= trackTPC->GetTgl();;
2548 AliESDfriendTrack* friendTrack=esdFriend->GetTrack(iTrack);
2549 if(!friendTrack) continue;
2550 if (friendTrack->GetITSOut()){
2551 const AliExternalTrackParam *trackITS=friendTrack->GetITSOut();
2552 trackITS->GetXYZAt(radiusMatch,bz,xyz);
2553 trackITS->GetPxPyPzAt(radiusMatch,bz,pxyz);
2554 for (Int_t i=0; i<3; i++){
2555 vecPosR0(iTrack,i)=xyz[i];
2556 vecMomR0(iTrack,i)=pxyz[i];
2558 vecPosR0(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]);
2559 vecMomR0(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]);
2560 vecMomR0(iTrack,4)= trackITS->GetSigned1Pt();;
2561 vecMomR0(iTrack,5)= trackITS->GetTgl();;
2565 // 1.) Find closest matching tracks, between the ITS standalone track
2566 // and the all other tracks
2567 // a.) caltegory - All
2568 // b.) category - without ITS
2571 Int_t ntracksPropagated=0;
2572 AliExternalTrackParam extTrackDummy;
2573 AliESDtrack esdTrackDummy;
2574 AliExternalTrackParam itsAtTPC;
2575 AliExternalTrackParam itsAtITSTPC;
2576 for (Int_t iTrack0=0; iTrack0<ntracks; iTrack0++){
2577 AliESDtrack *track0 = esdEvent->GetTrack(iTrack0);
2578 if(!track0) continue;
2579 if (track0->IsOn(AliVTrack::kTPCin)) continue;
2580 AliESDfriendTrack* friendTrack0=esdFriend->GetTrack(iTrack0);
2581 if (!friendTrack0) continue;
2582 //if (!track0->IsOn(AliVTrack::kITSpureSA)) continue;
2583 //if (!friendTrack0->GetITSOut()) continue; // is there flag for ITS standalone?
2584 ntracksPropagated++;
2586 // 2.) find clostest TPCtrack
2588 Double_t minChi2All=10000000;
2589 Double_t minChi2TPC=10000000;
2590 Double_t minChi2TPCITS=10000000;
2593 Int_t indexTPCITS=-1;
2594 Int_t ncandidates0=0; // n candidates - rough cut
2595 Int_t ncandidates1=0; // n candidates - rough + chi2 cut
2596 itsAtTPC=*(friendTrack0->GetITSOut());
2597 itsAtITSTPC=*(friendTrack0->GetITSOut());
2598 for (Int_t iTrack1=0; iTrack1<ntracks; iTrack1++){
2599 AliESDtrack *track1 = esdEvent->GetTrack(iTrack1);
2600 if(!track1) continue;
2601 if (!track1->IsOn(AliVTrack::kTPCin)) continue;
2604 if (TMath::Abs(vecPosR1(iTrack1,2)-vecPosR0(iTrack0,2))>dFastZCut) continue;
2605 if (TMath::Abs(vecPosR1(iTrack1,3)-vecPosR0(iTrack0,3))>dFastPosPhiCut) continue;
2606 if (TMath::Abs(vecMomR1(iTrack1,3)-vecMomR0(iTrack0,3))>dFastPhiCut) continue;
2607 if (TMath::Abs(vecMomR1(iTrack1,5)-vecMomR0(iTrack0,5))>dFastThetaCut) continue;
2608 if (TMath::Abs(vecMomR1(iTrack1,4)-vecMomR0(iTrack0,4))>dFastPtCut) continue;
2611 const AliExternalTrackParam * param1= track1->GetInnerParam();
2612 if (!friendTrack0->GetITSOut()) continue;
2613 AliExternalTrackParam outerITS = *(friendTrack0->GetITSOut());
2614 if (!outerITS.Rotate(param1->GetAlpha())) continue;
2615 if (!outerITS.PropagateTo(param1->GetX(),bz)) continue; // assume track close to the TPC inner wall
2616 Double_t chi2 = outerITS.GetPredictedChi2(param1);
2617 if (chi2>chi2Cut) continue;
2619 if (chi2<minChi2All){
2623 if (chi2<minChi2TPC && track1->IsOn(AliVTrack::kITSin)==0){
2628 if (chi2<minChi2TPCITS && track1->IsOn(AliVTrack::kITSin)){
2630 indexTPCITS=iTrack1;
2631 itsAtITSTPC=outerITS;
2635 AliESDtrack * trackAll= (indexAll>=0)? esdEvent->GetTrack(indexAll):&esdTrackDummy;
2636 AliESDtrack * trackTPC= (indexTPC>=0)? esdEvent->GetTrack(indexTPC):&esdTrackDummy;
2637 AliESDtrack * trackTPCITS= (indexTPCITS>=0)? esdEvent->GetTrack(indexTPCITS):&esdTrackDummy;
2638 (*fTreeSRedirector)<<"itsTPC"<<
2639 "indexAll="<<indexAll<< // index of closest track (chi2)
2640 "indexTPC="<<indexTPC<< // index of closest TPCalone tracks
2641 "indexTPCITS="<<indexTPCITS<< // index of closest cobined tracks
2642 "ncandidates0="<<ncandidates0<< // number of candidates
2643 "ncandidates1="<<ncandidates1<<
2645 "chi2All="<<minChi2All<< // chi2 of closest tracks
2646 "chi2TPC="<<minChi2TPC<<
2647 "chi2TPCITS="<<minChi2TPCITS<<
2649 "track0.="<<track0<< // ITS standalone tracks
2650 "trackAll.="<<trackAll<< // Closets other track
2651 "trackTPC.="<<trackTPC<< // Closest TPC only track
2652 "trackTPCITS.="<<trackTPCITS<< // closest combined track
2654 "itsAtTPC.="<<&itsAtTPC<< // ITS track parameters at the TPC alone track frame
2655 "itsAtITSTPC.="<<&itsAtITSTPC<< // ITS track parameters at the TPC combeined track frame
2660 void AliAnalysisTaskFilteredTree::ProcessTrackMatch(AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/){
2662 Use TPC standalone, ITS standalone and combined tracks to categorize/resp. recover track information.
2668 Topology categories:
2669 -Nice isolated tracks with full information
2670 -Overlapped tracks - Refit and sign them
2671 -Multiple found (check overlap factor) - Merge and sign
2672 -Charge particle (kink) decays - Change of direction - Sign them)
2674 -Array of indexes of closest tracks in each individual category
2675 -Chi2 distance to the closest tracks at reference radius of TPCin
2676 -Overlap factors - fraction of shared clusters or missing region
2677 -Chi2 distance at DCA
2679 -matrix closest tracks from each categories
2680 -characterization - chi2, index,chi2, overlap ratio
2683 0.) Kink decay or catastophic multiple scaterring
2684 (combining all track categories)
2685 - small chi2 at the DCA
2686 - significantly large deflection angle
2687 - Combinatorial algorithm - to decrease CPU time restriction of investigation to tracks with small DCA at refernce radius
2689 1.) if (TPConly && !(TPC+ITS) && ITSonly match ) {
2690 if (!kink) TPCOnly.addoptITS
2691 if (kink) TPConly sign
2694 2.) Overlap tracks - Refit with doUnfold