]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/AliAnalysisTaskFilteredTree.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
PMD info on da/amore updated
[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskFilteredTree.cxx
CommitLineData
c683985a 1/**************************************************************************
234cae2f 2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
c683985a 15
1a30582a 16/*
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
20 Usage:
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
26 +
27 6.) Info used for later raw data OFFLINE triggering (highPt, V0, laser, cosmic, high dEdx)
28
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
35*/
36
c683985a 37#include "iostream"
38#include "TSystem.h"
c683985a 39#include <TPDGCode.h>
40#include <TDatabasePDG.h>
41
42#include "TChain.h"
43#include "TTreeStream.h"
44#include "TTree.h"
45#include "TH1F.h"
46#include "TH3.h"
47#include "TCanvas.h"
48#include "TList.h"
49#include "TObjArray.h"
50#include "TFile.h"
51#include "TMatrixD.h"
52#include "TRandom3.h"
53
54#include "AliHeader.h"
55#include "AliGenEventHeader.h"
56#include "AliInputEventHandler.h"
57#include "AliStack.h"
58#include "AliTrackReference.h"
2e16763c 59#include "AliTrackPointArray.h"
60#include "AliSysInfo.h"
c683985a 61
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"
73
74#include "AliCentrality.h"
75#include "AliESDVZERO.h"
76#include "AliMultiplicity.h"
77
78#include "AliESDtrackCuts.h"
79#include "AliMCEventHandler.h"
80#include "AliFilteredTreeEventCuts.h"
81#include "AliFilteredTreeAcceptanceCuts.h"
82
83#include "AliAnalysisTaskFilteredTree.h"
84#include "AliKFParticle.h"
85#include "AliESDv0.h"
86#include "TVectorD.h"
87
88using namespace std;
89
90ClassImp(AliAnalysisTaskFilteredTree)
91
234cae2f 92 //_____________________________________________________________________________
93 AliAnalysisTaskFilteredTree::AliAnalysisTaskFilteredTree(const char *name)
c683985a 94 : AliAnalysisTaskSE(name)
95 , fESD(0)
96 , fMC(0)
97 , fESDfriend(0)
98 , fOutput(0)
99 , fPitList(0)
100 , fUseMCInfo(kFALSE)
101 , fUseESDfriends(kFALSE)
102 , fReducePileUp(kTRUE)
103 , fFillTree(kTRUE)
104 , fFilteredTreeEventCuts(0)
105 , fFilteredTreeAcceptanceCuts(0)
106 , fFilteredTreeRecAcceptanceCuts(0)
107 , fEsdTrackCuts(0)
108 , fTrigger(AliTriggerAnalysis::kMB1)
109 , fAnalysisMode(kTPCAnalysisMode)
110 , fTreeSRedirector(0)
111 , fCentralityEstimator(0)
112 , fLowPtTrackDownscaligF(0)
113 , fLowPtV0DownscaligF(0)
b4990577 114 , fFriendDownscaling(-3.)
c683985a 115 , fProcessAll(kFALSE)
116 , fProcessCosmics(kFALSE)
4a00e62d 117 , fProcessITSTPCmatchOut(kFALSE) // swittch to process ITS/TPC standalone tracks
c683985a 118 , fHighPtTree(0)
119 , fV0Tree(0)
120 , fdEdxTree(0)
121 , fLaserTree(0)
122 , fMCEffTree(0)
123 , fCosmicPairsTree(0)
124 , fPtResPhiPtTPC(0)
125 , fPtResPhiPtTPCc(0)
126 , fPtResPhiPtTPCITS(0)
127 , fPtResEtaPtTPC(0)
128 , fPtResEtaPtTPCc(0)
129 , fPtResEtaPtTPCITS(0)
130 , fPtResCentPtTPC(0)
131 , fPtResCentPtTPCc(0)
2e16763c 132 , fPtResCentPtTPCITS(0)
1a30582a 133 , fCurrentFileName("")
2e16763c 134 , fDummyTrack(0)
c683985a 135{
136 // Constructor
137
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());
146}
147
148//_____________________________________________________________________________
149AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree()
150{
1a30582a 151 //
152 // Destructor
153 //
c683985a 154 delete fFilteredTreeEventCuts;
155 delete fFilteredTreeAcceptanceCuts;
156 delete fFilteredTreeRecAcceptanceCuts;
157 delete fEsdTrackCuts;
158}
159
160//____________________________________________________________________________
161Bool_t AliAnalysisTaskFilteredTree::Notify()
162{
1a30582a 163 //
164 //
c683985a 165 static Int_t count = 0;
166 count++;
167 TTree *chain = (TChain*)GetInputData(0);
1a30582a 168 if(chain){
c683985a 169 Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());
1a30582a 170 }
171 fCurrentFileName=chain->GetCurrentFile()->GetName();
172 return kTRUE;
c683985a 173}
174
175//_____________________________________________________________________________
176void AliAnalysisTaskFilteredTree::UserCreateOutputObjects()
177{
178 // Create histograms
179 // Called once
180
181 //
182 //get the output file to make sure the trees will be associated to it
183 OpenFile(1);
184 fTreeSRedirector = new TTreeSRedirector();
234cae2f 185
c683985a 186 //
187 // Create trees
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();
194
82522b46 195 if (!fDummyTrack) {
2e16763c 196 fDummyTrack=new AliESDtrack();
197 }
c683985a 198
199 // histogram booking
200
201 Double_t minPt = 0.1;
202 Double_t maxPt = 100.;
203 Int_t nbinsPt = 30;
204
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];
209 binsPt[0] = minPt;
210 for (Int_t i=1;i<=nbinsPt;i++) {
211 binsPt[i] = minPt + TMath::Power(10,logminPt+i*binwidth);
212 }
213
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;
221 }
222
223 // phi bins
224 Double_t minPhi = 0.;
225 Double_t maxPhi = 6.5;
226 Int_t nbinsPhi = 100;
227 Double_t binsPhi[101];
234cae2f 228 for (Int_t i=0;i<=nbinsPhi;i++) {
c683985a 229 binsPhi[i] = minPhi + i*(maxPhi-minPhi)/nbinsPhi;
230 }
231
232 // eta bins
233 Double_t minEta = -1.;
234 Double_t maxEta = 1.;
235 Int_t nbinsEta = 20;
236 Double_t binsEta[21];
237 for (Int_t i=0;i<=nbinsEta;i++) {
238 binsEta[i] = minEta + i*(maxEta-minEta)/nbinsEta;
239 }
240
241 // mult bins
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;
248 }
234cae2f 249
c683985a 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);
234cae2f 253
254 fPtResEtaPtTPC = new TH3D("fPtResEtaPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
c683985a 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);
234cae2f 257
258 fPtResCentPtTPC = new TH3D("fPtResCentPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
c683985a 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);
261
234cae2f 262
c683985a 263 fOutput = new TList;
264 if(!fOutput) return;
265 fOutput->SetOwner();
266
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);
276
277 // post data to outputs
278
279 PostData(1,fV0Tree);
280 PostData(2,fHighPtTree);
281 PostData(3,fdEdxTree);
282 PostData(4,fLaserTree);
283 PostData(5,fMCEffTree);
284 PostData(6,fCosmicPairsTree);
285
286 PostData(7,fOutput);
287}
288
289//_____________________________________________________________________________
290void AliAnalysisTaskFilteredTree::UserExec(Option_t *)
291{
292 //
293 // Called for each event
294 //
295
296 // ESD event
297 fESD = (AliESDEvent*) (InputEvent());
298 if (!fESD) {
299 Printf("ERROR: ESD event not available");
300 return;
301 }
c683985a 302 //if MC info available - use it.
303 fMC = MCEvent();
c683985a 304 if(fUseESDfriends) {
234cae2f 305 //fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
306 fESDfriend = ESDfriend();
234cae2f 307 if(!fESDfriend) {
308 Printf("ERROR: ESD friends not available");
c683985a 309 }
310 }
1a30582a 311 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
312 if (!inputHandler){
313 return;
314 }
c683985a 315
316 //if set, use the environment variables to set the downscaling factors
317 //AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF
318 //AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF
287a52d0 319 //AliAnalysisTaskFilteredTree_fFriendDownscaling
c683985a 320 TString env;
321 env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF");
1a30582a 322 if (!env.IsNull()){
c683985a 323 fLowPtTrackDownscaligF=env.Atof();
324 AliInfo(Form("fLowPtTrackDownscaligF=%f",fLowPtTrackDownscaligF));
325 }
326 env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF");
1a30582a 327 if (!env.IsNull()){
c683985a 328 fLowPtV0DownscaligF=env.Atof();
329 AliInfo(Form("fLowPtV0DownscaligF=%f",fLowPtTrackDownscaligF));
330 }
287a52d0 331 env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fFriendDownscaling");
1a30582a 332 if (!env.IsNull()){
287a52d0 333 fFriendDownscaling=env.Atof();
334 AliInfo(Form(" fFriendDownscaling=%f",fFriendDownscaling));
335 }
1a30582a 336 //
337 //
c683985a 338 //
339 if(fProcessAll) {
340 ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC
341 }
342 else {
343 Process(fESD,fMC,fESDfriend); // only global and TPC tracks
344 }
c683985a 345 //
346 ProcessV0(fESD,fMC,fESDfriend);
347 ProcessLaser(fESD,fMC,fESDfriend);
348 ProcessdEdx(fESD,fMC,fESDfriend);
234cae2f 349 if (fProcessCosmics) { ProcessCosmics(fESD,fESDfriend); }
c683985a 350 if(fMC) { ProcessMCEff(fESD,fMC,fESDfriend);}
4a00e62d 351 if (fProcessITSTPCmatchOut) ProcessITSTPCmatchOut(fESD, fESDfriend);
1236ed2e 352 printf("processed event %d\n", Int_t(Entry()));
c683985a 353}
354
355//_____________________________________________________________________________
234cae2f 356void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliESDfriend* esdFriend)
c683985a 357{
358 //
1a30582a 359 // Find cosmic pairs (triggered or random)
234cae2f 360 //
361 //
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();
234cae2f 369 UInt_t specie = event->GetEventSpecie(); // skip laser events
370 if (specie==AliRecoParam::kCalib) return;
c0e33396 371 Int_t ntracksFriend = esdFriend ? esdFriend->GetNumberOfTracks() : 0;
234cae2f 372
373
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;
378
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
385 //rm primaries
c683985a 386 //
234cae2f 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();
234cae2f 391 AliESDfriendTrack* friendTrack0=NULL;
dfd5cb6f 392 if (esdFriend &&!esdFriend->TestSkipBit()){
393 if (itrack0<ntracksFriend){
394 friendTrack0 = esdFriend->GetTrack(itrack0);
395 } //this guy can be NULL
396 }
234cae2f 397
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;
c683985a 410 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
234cae2f 411 //
412 const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
413 //
414 Bool_t isPair=kTRUE;
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;
418 }
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
422 /*
82522b46 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"
426 */
234cae2f 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();
234cae2f 433 //
434 //
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();
442
b55a470b 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);
448
c683985a 449
234cae2f 450 AliESDfriendTrack* friendTrack1=NULL;
dfd5cb6f 451 if (esdFriend &&!esdFriend->TestSkipBit()){
452 if (itrack1<ntracksFriend){
453 friendTrack1 = esdFriend->GetTrack(itrack1);
454 } //this guy can be NULL
455 }
456
82522b46 457 //
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;
464 }
465 }
466 if (fFriendDownscaling<=0){
467 if (((*fTreeSRedirector)<<"CosmicPairs").GetTree()){
468 TTree * tree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
469 if (tree){
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) {
476 friendTrackStore0=0;
477 friendTrackStore1=0;
478 }
479 }
480 }
481 }
c683985a 482 if(!fFillTree) return;
483 if(!fTreeSRedirector) return;
234cae2f 484 (*fTreeSRedirector)<<"CosmicPairs"<<
1a30582a 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
234cae2f 493 //
1a30582a 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
234cae2f 502 "\n";
234cae2f 503 }
234cae2f 504 }
c683985a 505}
506
507
508//_____________________________________________________________________________
509void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
510{
511 //
512 // Select real events with high-pT tracks
513 //
c683985a 514
515 // get selection cuts
516 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
517 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
518 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
519
520 if(!evtCuts || !accCuts || !esdTrackCuts) {
521 Printf("ERROR cuts not available");
522 return;
523 }
524
525 // trigger selection
526 Bool_t isEventTriggered = kTRUE;
527 AliPhysicsSelection *physicsSelection = NULL;
528 AliTriggerAnalysis* triggerAnalysis = NULL;
529
530 //
531 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
c683985a 532 // trigger
533 if(evtCuts->IsTriggerRequired())
534 {
535 // always MB
536 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
537
538 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
539 if(!physicsSelection) return;
540 //SetPhysicsTriggerSelection(physicsSelection);
541
542 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
543 // set trigger (V0AND)
544 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
545 if(!triggerAnalysis) return;
546 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
547 }
548 }
549
550 // centrality determination
551 Float_t centralityF = -1;
552 AliCentrality *esdCentrality = esdEvent->GetCentrality();
553 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
554
555 // use MC information
556 AliHeader* header = 0;
557 AliGenEventHeader* genHeader = 0;
558 AliStack* stack = 0;
559 TArrayF vtxMC(3);
560
561 Int_t multMCTrueTracks = 0;
562 if(mcEvent)
563 {
564 // get MC event header
565 header = mcEvent->Header();
566 if (!header) {
567 AliDebug(AliLog::kError, "Header not available");
568 return;
569 }
570 // MC particle stack
571 stack = mcEvent->Stack();
572 if (!stack) {
573 AliDebug(AliLog::kError, "Stack not available");
574 return;
575 }
576
577 // get MC vertex
578 genHeader = header->GenEventHeader();
579 if (!genHeader) {
580 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
581 return;
582 }
583 genHeader->PrimaryVertex(vtxMC);
584
585 // multipliticy of all MC primary tracks
586 // in Zv, pt and eta ranges)
587 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
588 } // end bUseMC
589
590 // get reconstructed vertex
591 //const AliESDVertex* vtxESD = 0;
592 AliESDVertex* vtxESD = 0;
593 if(GetAnalysisMode() == kTPCAnalysisMode) {
234cae2f 594 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
c683985a 595 }
596 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
234cae2f 597 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
c683985a 598 }
599 else {
234cae2f 600 return;
c683985a 601 }
602
603 AliESDVertex* vtxTPC = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
604 AliESDVertex* vtxSPD = (AliESDVertex*)esdEvent->GetPrimaryVertexSPD();
605
606 if(!vtxESD) return;
607 if(!vtxTPC) return;
608 if(!vtxSPD) return;
609
610 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
e690d4d0 611 //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZ());
c683985a 612 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
613 Int_t ntracks = esdEvent->GetNumberOfTracks();
614
615 // check event cuts
616 if(isEventOK && isEventTriggered)
617 {
618
619 //
620 // get IR information
621 //
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
627
628 // Use when Peter commit the changes in the header
629 Int_t ir1 = 0;
630 Int_t ir2 = 0;
631
632 //
07aef7d1 633 //Double_t vert[3] = {0};
e690d4d0 634 //vert[0] = vtxESD->GetX();
635 //vert[1] = vtxESD->GetY();
636 //vert[2] = vtxESD->GetZ();
c683985a 637 Int_t mult = vtxESD->GetNContributors();
638 Int_t multSPD = vtxSPD->GetNContributors();
639 Int_t multTPC = vtxTPC->GetNContributors();
640
641 Float_t bz = esdEvent->GetMagneticField();
642 Int_t runNumber = esdEvent->GetRunNumber();
643 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
644 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
b55a470b 645
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);
651
c683985a 652
653 // high pT tracks
654 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
655 {
656 AliESDtrack *track = esdEvent->GetTrack(iTrack);
657 if(!track) continue;
658 if(track->Charge()==0) continue;
659 if(!esdTrackCuts->AcceptTrack(track)) continue;
660 if(!accCuts->AcceptTrack(track)) continue;
234cae2f 661
c683985a 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);
668
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());
675 if(!isOK) continue;
676
677 // Dump to the tree
678 // vertex
679 // TPC-ITS tracks
680 //
681 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
c683985a 682 if(!fFillTree) return;
683 if(!fTreeSRedirector) return;
684 (*fTreeSRedirector)<<"highPt"<<
b55a470b 685 "gid="<<gid<<
1a30582a 686 "fileName.="<<&fCurrentFileName<<
c683985a 687 "runNumber="<<runNumber<<
688 "evtTimeStamp="<<evtTimeStamp<<
689 "evtNumberInFile="<<evtNumberInFile<<
234cae2f 690 "triggerClass="<<&triggerClass<< // trigger
c683985a 691 "Bz="<<bz<< // magnetic field
692 "vtxESD.="<<vtxESD<<
234cae2f 693 "ntracksESD="<<ntracks<< // number of tracks in the ESD
234cae2f 694 "IRtot="<<ir1<< // interaction record history info
695 "IRint2="<<ir2<<
c683985a 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<<
701 "\n";
702 }
703 }
234cae2f 704
c683985a 705}
706
707
708//_____________________________________________________________________________
234cae2f 709void AliAnalysisTaskFilteredTree::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const esdFriend)
c683985a 710{
711 //
82522b46 712 // Process laser events -> dump tracks and clusters to the special tree
c683985a 713 //
1236ed2e 714 const Double_t kMinPt = 5;
234cae2f 715 if(!fFillTree) return;
716 if(!fTreeSRedirector) return;
c683985a 717 const AliESDHeader* esdHeader = esdEvent->GetHeader();
1a30582a 718 if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) {
c683985a 719 Int_t countLaserTracks = 0;
234cae2f 720 Int_t runNumber = esdEvent->GetRunNumber();
721 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
722 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
723 Float_t bz = esdEvent->GetMagneticField();
82522b46 724 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
b55a470b 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);
1a30582a 730 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
c683985a 731 AliESDtrack *track = esdEvent->GetTrack(iTrack);
732 if(!track) continue;
1a30582a 733 if(track->GetTPCInnerParam()) countLaserTracks++;
234cae2f 734 AliESDfriendTrack* friendTrack=NULL;
1236ed2e 735 // suppress beam background and CE random reacks
15120578 736 if (track->GetInnerParam()->Pt()<kMinPt) continue;
1236ed2e 737 Bool_t skipTrack=gRandom->Rndm()>1/(1+TMath::Abs(fFriendDownscaling));
738 if (skipTrack) continue;
1a30582a 739 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
c683985a 740 (*fTreeSRedirector)<<"Laser"<<
82522b46 741 "gid="<<gid<< // global identifier of event
1a30582a 742 "fileName.="<<&fCurrentFileName<< //
c683985a 743 "runNumber="<<runNumber<<
744 "evtTimeStamp="<<evtTimeStamp<<
745 "evtNumberInFile="<<evtNumberInFile<<
82522b46 746 "triggerClass="<<&triggerClass<< // trigger
747 "Bz="<<bz<< // magnetic field
748 "multTPCtracks="<<countLaserTracks<< // multiplicity of tracks
1236ed2e 749 "track.="<<track<< // track parameters
82522b46 750 "friendTrack.="<<friendTrack<< // friend track information
c683985a 751 "\n";
752 }
753 }
754}
755
756//_____________________________________________________________________________
757void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
758{
759 //
760 // Select real events with high-pT tracks
761 // Calculate and stor in the output tree:
762 // TPC constrained tracks
763 // InnerParams constrained tracks
764 // TPC-ITS 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
769 // MC information:
770 // track references at ITSin, TPCin; InnerParam at first TPC track reference,
771 // particle ID, mother ID, production mechanism ...
772 //
c683985a 773 // get selection cuts
774 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
775 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
776 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
777
778 if(!evtCuts || !accCuts || !esdTrackCuts) {
779 AliDebug(AliLog::kError, "cuts not available");
780 return;
781 }
782
783 // trigger selection
784 Bool_t isEventTriggered = kTRUE;
785 AliPhysicsSelection *physicsSelection = NULL;
786 AliTriggerAnalysis* triggerAnalysis = NULL;
787
788 //
789 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
234cae2f 790
c683985a 791
792 // trigger
793 if(evtCuts->IsTriggerRequired())
794 {
795 // always MB
796 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
797
798 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
799 if(!physicsSelection) {AliInfo("no physics selection"); return;}
800 //SetPhysicsTriggerSelection(physicsSelection);
801
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());
807 }
808 }
809
810 // centrality determination
811 Float_t centralityF = -1;
812 AliCentrality *esdCentrality = esdEvent->GetCentrality();
813 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
814
815 // use MC information
816 AliHeader* header = 0;
817 AliGenEventHeader* genHeader = 0;
818 AliStack* stack = 0;
819 TArrayF vtxMC(3);
820 Int_t mcStackSize=0;
821
822 Int_t multMCTrueTracks = 0;
823 if(mcEvent)
824 {
825 // get MC event header
826 header = mcEvent->Header();
827 if (!header) {
828 AliDebug(AliLog::kError, "Header not available");
829 return;
830 }
831 // MC particle stack
832 stack = mcEvent->Stack();
833 if (!stack) {
834 AliDebug(AliLog::kError, "Stack not available");
835 return;
836 }
837 mcStackSize=stack->GetNtrack();
838
839 // get MC vertex
840 genHeader = header->GenEventHeader();
841 if (!genHeader) {
842 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
843 return;
844 }
845 genHeader->PrimaryVertex(vtxMC);
846
847 // multipliticy of all MC primary tracks
848 // in Zv, pt and eta ranges)
849 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
850
851 } // end bUseMC
852
853 // get reconstructed vertex
854 //const AliESDVertex* vtxESD = 0;
855 AliESDVertex* vtxESD = 0;
856 if(GetAnalysisMode() == kTPCAnalysisMode) {
234cae2f 857 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
c683985a 858 }
859 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
234cae2f 860 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
c683985a 861 }
862 else {
863 AliInfo("no ESD vertex");
234cae2f 864 return;
c683985a 865 }
866
867 if(!vtxESD) return;
868
869 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
234cae2f 870
c683985a 871 //
872 // Vertex info comparison and track multiplicity
873 //
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());
881 Int_t ntracksTPC=0;
882 Int_t ntracksITS=0;
883 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
884 AliESDtrack *track = esdEvent->GetTrack(iTrack);
885 if(!track) continue;
886 if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
887 if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
888 }
234cae2f 889
c683985a 890 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
891 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
892
893
894 // check event cuts
895 if(isEventOK && isEventTriggered)
896 {
897 //
898 // get IR information
899 //
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
905 //
906 Int_t ir1 = 0;
907 Int_t ir2 = 0;
908
909 //
910 Double_t vert[3] = {0};
e690d4d0 911 vert[0] = vtxESD->GetX();
912 vert[1] = vtxESD->GetY();
913 vert[2] = vtxESD->GetZ();
c683985a 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();
919
2e16763c 920 Int_t numberOfTracks=esdEvent->GetNumberOfTracks();
c683985a 921 // high pT tracks
2e16763c 922 for (Int_t iTrack = 0; iTrack < numberOfTracks; iTrack++)
c683985a 923 {
924 AliESDtrack *track = esdEvent->GetTrack(iTrack);
234cae2f 925 AliESDfriendTrack* friendTrack=NULL;
2e16763c 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
c683985a 929 if(!track) continue;
930 if(track->Charge()==0) continue;
931 if(!esdTrackCuts->AcceptTrack(track)) continue;
932 if(!accCuts->AcceptTrack(track)) continue;
234cae2f 933
c683985a 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);
940
941 // Dump to the tree
942 // vertex
943 // TPC constrained tracks
944 // InnerParams constrained tracks
945 // TPC-ITS 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
950 // MC information
234cae2f 951
c683985a 952 Double_t x[3]; track->GetXYZ(x);
953 Double_t b[3]; AliTracker::GetBxByBz(x,b);
954
955 //
956 // Transform TPC inner params to track reference frame
957 //
958 Bool_t isOKtpcInner = kFALSE;
959 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
960 if (tpcInner) {
961 // transform to the track reference frame
962 isOKtpcInner = tpcInner->Rotate(track->GetAlpha());
963 isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
964 }
965
966 //
967 // Constrain TPC inner to vertex
968 // clone TPCinner has to be deleted
969 //
970 Bool_t isOKtpcInnerC = kFALSE;
971 AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
972 if (tpcInnerC) {
973 isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);
974 isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());
975 isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
976 }
977
978 //
979 // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex
980 // Clone track InnerParams has to be deleted
981 //
982 Bool_t isOKtrackInnerC = kFALSE;
983 AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));
984 if (trackInnerC) {
985 isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);
986 isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());
987 isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
988 }
234cae2f 989
c683985a 990 //
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)
994 //
995 TMatrixD deltaT(5,1), deltaTtrackC(5,1);
996 TMatrixD delta(1,5), deltatrackC(1,5);
997 TMatrixD covarM(5,5), covarMtrackC(5,5);
998 TMatrixD chi2(1,1);
999 TMatrixD chi2trackC(1,1);
1000
1001 if(isOKtpcInnerC && isOKtrackInnerC)
1002 {
1003 for (Int_t ipar=0; ipar<5; ipar++) {
1004 deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
234cae2f 1005 delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
c683985a 1006
1007 deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
234cae2f 1008 deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
c683985a 1009
1010 for (Int_t jpar=0; jpar<5; jpar++) {
234cae2f 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];
c683985a 1014 }
1015 }
1016
1017 // chi2 distance TPC constrained and TPC+ITS
1018 TMatrixD covarMInv = covarM.Invert();
1019 TMatrixD mat2 = covarMInv*deltaT;
1020 chi2 = delta*mat2;
1021 //chi2.Print();
1022
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();
1028 }
1029
1030
1031 //
1032 // Propagate ITSout to TPC inner wall
1033 // and calculate chi2 distance to track (InnerParams)
1034 //
1035 const Double_t kTPCRadius=85;
1036 const Double_t kStep=3;
1037
1038 // clone track InnerParams has to be deleted
1039 Bool_t isOKtrackInnerC2 = kFALSE;
1040 AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));
1041 if (trackInnerC2) {
1042 isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);
1043 }
1044
1045 Bool_t isOKouterITSc = kFALSE;
1046 AliExternalTrackParam *outerITSc = NULL;
1047 TMatrixD chi2OuterITS(1,1);
1048
234cae2f 1049 if(esdFriend && !esdFriend->TestSkipBit())
c683985a 1050 {
1051 // propagate ITSout to TPC inner wall
c683985a 1052 if(friendTrack)
234cae2f 1053 {
1054
1055 outerITSc = NULL;
1056 if (friendTrack->GetITSOut()) outerITSc = new AliExternalTrackParam(*(friendTrack->GetITSOut()));
c683985a 1057 if(outerITSc)
234cae2f 1058 {
c683985a 1059 isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);
1060 isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());
1061 isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());
1062
234cae2f 1063 //
c683985a 1064 // calculate chi2 between outerITS and innerParams
234cae2f 1065 // cov without z-coordinate at the moment
1066 //
c683985a 1067 TMatrixD deltaTouterITS(4,1);
1068 TMatrixD deltaouterITS(1,4);
1069 TMatrixD covarMouterITS(4,4);
1070
1071 if(isOKtrackInnerC2 && isOKouterITSc) {
234cae2f 1072 Int_t kipar = 0;
1073 Int_t kjpar = 0;
c683985a 1074 for (Int_t ipar=0; ipar<5; ipar++) {
234cae2f 1075 if(ipar!=1) {
c683985a 1076 deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
234cae2f 1077 deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
1078 }
c683985a 1079
1080 kjpar=0;
1081 for (Int_t jpar=0; jpar<5; jpar++) {
234cae2f 1082 Int_t index=outerITSc->GetIndex(ipar,jpar);
1083 if(ipar !=1 || jpar!=1) {
1084 covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];
1085 }
c683985a 1086 if(jpar!=1) kjpar++;
234cae2f 1087 }
1088 if(ipar!=1) kipar++;
1089 }
c683985a 1090
1091 // chi2 distance ITSout and InnerParams
1092 TMatrixD covarMInvouterITS = covarMouterITS.Invert();
1093 TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;
1094 chi2OuterITS = deltaouterITS*mat2outerITS;
1095 //chi2OuterITS.Print();
234cae2f 1096 }
c683985a 1097 }
1098 }
1099 }
1100
1101 //
1102 // MC info
1103 //
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;
1111
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;
1120
1121 Bool_t isOKtrackInnerC3 = kFALSE;
1122 AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
2e16763c 1123 if(mcEvent && stack)
c683985a 1124 {
2e16763c 1125 do //artificial loop (once) to make the continue statements jump out of the MC part
234cae2f 1126 {
2e16763c 1127 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1128 //
1129 // global track
1130 //
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.)
1136 {
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);
1143 }
c683985a 1144
2e16763c 1145 //
1146 // TPC track
1147 //
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.)
1153 {
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);
1160 }
c683985a 1161
2e16763c 1162 //
1163 // store first track reference
1164 // for TPC track
1165 //
1166 TParticle *part=0;
1167 TClonesArray *trefs=0;
1168 Int_t status = mcEvent->GetParticleAndTR(TMath::Abs(labelTPC), part, trefs);
c683985a 1169
2e16763c 1170 if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.)
c683985a 1171 {
2e16763c 1172 Int_t nTrackRef = trefs->GetEntries();
1173 //printf("nTrackRef %d \n",nTrackRef);
c683985a 1174
2e16763c 1175 Int_t countITS = 0;
1176 for (Int_t iref = 0; iref < nTrackRef; iref++)
c683985a 1177 {
2e16763c 1178 AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);
1179
1180 // Ref. in the middle ITS
1181 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kITS)
1182 {
1183 nrefITS++;
1184 if(!refITS && countITS==2) {
1185 refITS = ref;
1186 //printf("refITS %p \n",refITS);
1187 }
1188 countITS++;
234cae2f 1189 }
c683985a 1190
2e16763c 1191 // TPC
1192 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC)
1193 {
1194 nrefTPC++;
1195 refTPCOut=ref;
1196 if(!refTPCIn) {
1197 refTPCIn = ref;
1198 //printf("refTPCIn %p \n",refTPCIn);
1199 //break;
1200 }
234cae2f 1201 }
2e16763c 1202 // TRD
1203 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD)
1204 {
1205 nrefTRD++;
1206 if(!refTRD) {
1207 refTRD = ref;
1208 }
234cae2f 1209 }
2e16763c 1210 // TOF
1211 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF)
1212 {
1213 nrefTOF++;
1214 if(!refTOF) {
1215 refTOF = ref;
1216 }
234cae2f 1217 }
2e16763c 1218 // EMCAL
1219 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL)
1220 {
1221 nrefEMCAL++;
1222 if(!refEMCAL) {
1223 refEMCAL = ref;
1224 }
1225 }
1226 // PHOS
1227 // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS)
1228 // {
1229 // nrefPHOS++;
1230 // if(!refPHOS) {
1231 // refPHOS = ref;
1232 // }
1233 // }
c683985a 1234 }
2e16763c 1235
1236 // transform inner params to TrackRef
1237 // reference frame
1238 if(refTPCIn && trackInnerC3)
c683985a 1239 {
2e16763c 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);
c683985a 1243 }
234cae2f 1244 }
c683985a 1245
2e16763c 1246 //
1247 // ITS track
1248 //
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.)
234cae2f 1254 {
2e16763c 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);
234cae2f 1261 }
c683985a 1262 }
2e16763c 1263 while (0);
c683985a 1264 }
1265
234cae2f 1266 //
2e16763c 1267 Bool_t dumpToTree=kFALSE;
1268
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();
1274 if (fReducePileUp){
1275 //
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
1278 //
1279 Float_t dcaTPC[2];
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)){
1285 dumpToTree=kFALSE;
1286 }
234cae2f 1287 }
c683985a 1288
2e16763c 1289 //init dummy objects
b3982e10 1290 static AliESDVertex dummyvtxESD;
1291 //if (!dummyvtxESD)
1292 //{
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);
1298 //}
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();
2e16763c 1305
1306 //assign the dummy objects if needed
b3982e10 1307 if (!track) {track=fDummyTrack;}
b4990577 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;
1311 }
1312 if (fFriendDownscaling<=0){
1313 if (((*fTreeSRedirector)<<"highPt").GetTree()){
1314 TTree * tree = ((*fTreeSRedirector)<<"highPt").GetTree();
1315 if (tree){
1316 Double_t sizeAll=tree->GetZipBytes();
1317 TBranch * br= tree->GetBranch("friendTrack.fPoints");
1318 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
82522b46 1319 br= tree->GetBranch("friendTrack.fCalibContainer");
b4990577 1320 if (br) sizeFriend+=br->GetZipBytes();
1321 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) friendTrackStore=0;
1322 }
1323 }
1324 }
1325
1326
1327 // if (!friendTrackStore && fFriendDownscaling<=1) {friendTrack=fDummyFriendTrack;}
b3982e10 1328 if (!vtxESD) {vtxESD=&dummyvtxESD;}
c683985a 1329 if (mcEvent)
1330 {
b3982e10 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;}
c683985a 1339 }
2e16763c 1340
1341 Int_t ntracks = esdEvent->GetNumberOfTracks();
1342
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);
1348
1349 // fill histograms
1350 FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0));
1351
1a30582a 1352 if(fTreeSRedirector && dumpToTree && fFillTree) {
2e16763c 1353 (*fTreeSRedirector)<<"highPt"<<
1354 "gid="<<gid<<
1a30582a 1355 "fileName.="<<&fCurrentFileName<< // name of the chunk file (hopefully full)
2e16763c 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)
2e16763c 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
1373 //
1374 "esdTrack.="<<track<< // esdTrack as used in the physical analysis
b4990577 1375 // "friendTrack.="<<friendTrack<< // esdFriendTrack associated to the esdTrack
1376 "friendTrack.="<<friendTrackStore<< // esdFriendTrack associated to the esdTrack
2e16763c 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;
1a30582a 1386 if (mcEvent){
b3982e10 1387 static AliTrackReference refDummy;
2e16763c 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<<
1410 "mech="<<mech<<
1411 "isPrim="<<isPrim<<
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;
1429 }
1430 //finish writing the entry
1431 AliInfo("writing tree highPt");
1432 (*fTreeSRedirector)<<"highPt"<<"\n";
1433 }
82522b46 1434 AliSysInfo::AddStamp("filteringTask",iTrack,numberOfTracks,numberOfFriendTracks,(friendTrackStore)?0:1);
2e16763c 1435 delete tpcInnerC;
1436 delete trackInnerC;
1437 delete trackInnerC2;
1438 delete outerITSc;
1439 delete trackInnerC3;
1440 }
c683985a 1441 }
c683985a 1442}
1443
1444
1445//_____________________________________________________________________________
1446void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1447{
1448 //
1449 // Fill tree for efficiency studies MC only
1450 AliInfo("we start!");
c683985a 1451 if(!mcEvent) {
1452 AliDebug(AliLog::kError, "mcEvent not available");
1453 return;
1454 }
1455
1456 // get selection cuts
1457 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1458 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1459 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1460
1461 if(!evtCuts || !accCuts || !esdTrackCuts) {
1462 AliDebug(AliLog::kError, "cuts not available");
1463 return;
1464 }
1465
1466 // trigger selection
1467 Bool_t isEventTriggered = kTRUE;
1468 AliPhysicsSelection *physicsSelection = NULL;
1469 AliTriggerAnalysis* triggerAnalysis = NULL;
1470
1471 //
1472 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
c683985a 1473
c683985a 1474
1475 // trigger
1476 if(evtCuts->IsTriggerRequired())
1477 {
1478 // always MB
1479 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
c683985a 1480 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1481 if(!physicsSelection) return;
1482
1483 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1484 // set trigger (V0AND)
1485 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1486 if(!triggerAnalysis) return;
1487 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1488 }
1489 }
1490
1491 // centrality determination
1492 Float_t centralityF = -1;
1493 AliCentrality *esdCentrality = esdEvent->GetCentrality();
1494 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1495
1496 // use MC information
1497 AliHeader* header = 0;
1498 AliGenEventHeader* genHeader = 0;
1499 AliStack* stack = 0;
1500 Int_t mcStackSize=0;
1501 TArrayF vtxMC(3);
1502
1503 Int_t multMCTrueTracks = 0;
1504 //
1505 if(!mcEvent) {
1506 AliDebug(AliLog::kError, "mcEvent not available");
1507 return;
1508 }
1509 // get MC event header
1510 header = mcEvent->Header();
1511 if (!header) {
1512 AliDebug(AliLog::kError, "Header not available");
1513 return;
1514 }
1515 // MC particle stack
1516 stack = mcEvent->Stack();
1517 if (!stack) {
1518 AliDebug(AliLog::kError, "Stack not available");
1519 return;
1520 }
1521 mcStackSize=stack->GetNtrack();
1522
1523 // get MC vertex
1524 genHeader = header->GenEventHeader();
1525 if (!genHeader) {
1526 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1527 return;
1528 }
1529 genHeader->PrimaryVertex(vtxMC);
1530
1531 // multipliticy of all MC primary tracks
1532 // in Zv, pt and eta ranges)
1533 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1534
1535
1536 // get reconstructed vertex
1537 //const AliESDVertex* vtxESD = 0;
1538 AliESDVertex* vtxESD = 0;
1539 if(GetAnalysisMode() == kTPCAnalysisMode) {
1540 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1541 }
1542 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1543 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1544 }
1545 else {
1546 return;
1547 }
1548
1549 if(!vtxESD) return;
1550 //
1551 // Vertex info comparison and track multiplicity
1552 //
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());
1560 Int_t ntracksTPC=0;
1561 Int_t ntracksITS=0;
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++;
1567 }
1568
1569 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1570 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1571 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1572
1573 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1574
1575 // check event cuts
1576 if(isEventOK && isEventTriggered)
1577 {
1578 //if(!stack) return;
1579
1580 //
1581 // MC info
1582 //
1583 TParticle *particle=NULL;
1584 TParticle *particleMother=NULL;
1585 Int_t mech=-1;
1586
1587 // reco event info
1588 Double_t vert[3] = {0};
e690d4d0 1589 vert[0] = vtxESD->GetX();
1590 vert[1] = vtxESD->GetY();
1591 vert[2] = vtxESD->GetZ();
c683985a 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();
234cae2f 1597 // loop over MC stack
c683985a 1598 for (Int_t iMc = 0; iMc < mcStackSize; ++iMc)
1599 {
1600 particle = stack->Particle(iMc);
1601 if (!particle)
1602 continue;
1603
1604 // only charged particles
1605 if(!particle->GetPDG()) continue;
1606 Double_t charge = particle->GetPDG()->Charge()/3.;
1607 if (TMath::Abs(charge) < 0.001)
1608 continue;
1609
1610 // only primary particles
1611 Bool_t prim = stack->IsPhysicalPrimary(iMc);
1612 if(!prim) continue;
1613
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;
1619
1620 // is particle in acceptance
1621 if(!accCuts->AcceptTrack(particle)) continue;
1622
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;
1632
1633 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1634 {
1635 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1636 if(!track) continue;
1637 if(track->Charge()==0) continue;
234cae2f 1638 //
1639 Int_t label = TMath::Abs(track->GetLabel());
1640 if (label >= mcStackSize) continue;
1641 if(label == iMc) {
1642 Bool_t isAcc=esdTrackCuts->AcceptTrack(track);
1643 if (isAcc) isESDtrackCut=1;
1644 if (accCuts->AcceptTrack(track)) isAccCuts=1;
1645 isRec = kTRUE;
c683985a 1646 trackIndex = iTrack;
1647
234cae2f 1648 if (recTrack){
1649 if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track
1650 if (!isAcc) continue;
1651 trackLoopIndex = iTrack;
1652 }
1653 recTrack = esdEvent->GetTrack(trackIndex);
c683985a 1654 nRec++;
1655 if(track->GetLabel()<0) nFakes++;
1656
234cae2f 1657 continue;
1658 }
c683985a 1659 }
234cae2f 1660
c683985a 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);
1666 } else {
2e16763c 1667 recTrack = fDummyTrack;
c683985a 1668 }
1669
1670 particleMother = GetMother(particle,stack);
1671 mech = particle->GetUniqueID();
1672
1673 //MC particle track length
1674 Double_t tpcTrackLength = 0.;
1675 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
1676 if(mcParticle) {
1677 Int_t counter;
1678 tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
1679 }
1680
1681
1682 //
1683 if(fTreeSRedirector && fFillTree) {
1684 (*fTreeSRedirector)<<"MCEffTree"<<
1a30582a 1685 "fileName.="<<&fCurrentFileName<<
c683985a 1686 "triggerClass.="<<&triggerClass<<
1687 "runNumber="<<runNumber<<
1688 "evtTimeStamp="<<evtTimeStamp<<
1689 "evtNumberInFile="<<evtNumberInFile<< //
1690 "Bz="<<bz<< // magnetic field
1691 "vtxESD.="<<vtxESD<< // vertex info
234cae2f 1692 //
c683985a 1693 "mult="<<mult<< // primary vertex 9whatewe found) multiplicity
234cae2f 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
1702 //
1703 //
1704 "isAcc0="<<isESDtrackCut<< // track accepted by ESD track cuts
1705 "isAcc1="<<isAccCuts<< // track accepted by acceptance cuts flag
c683985a 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
1714 "\n";
1715 }
1716
2e16763c 1717 //if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
c683985a 1718 }
1719 }
1720
1721}
1722
1723//_____________________________________________________________________________
1724Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {
1725 //
1726 // check if particle is Z > 1
1727 //
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);
1733 //
1734
1735 Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
1736
1737 if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
1738
1739 return kFALSE;
1740}
1741
1742//_____________________________________________________________________________
234cae2f 1743void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
c683985a 1744{
1745 //
1746 // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates
1747 //
1748 if(!esdEvent) {
1749 AliDebug(AliLog::kError, "esdEvent not available");
1750 return;
1751 }
1752
1753 // get selection cuts
1754 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1755 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1756 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1757
1758 if(!evtCuts || !accCuts || !esdTrackCuts) {
1759 AliDebug(AliLog::kError, "cuts not available");
1760 return;
1761 }
1762
1763 // trigger selection
1764 Bool_t isEventTriggered = kTRUE;
1765 AliPhysicsSelection *physicsSelection = NULL;
1766 AliTriggerAnalysis* triggerAnalysis = NULL;
1767
1768 //
1769 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
c683985a 1770
1771 // trigger
1772 if(evtCuts->IsTriggerRequired())
1773 {
1774 // always MB
1775 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1776
1777 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1778 if(!physicsSelection) return;
1779 //SetPhysicsTriggerSelection(physicsSelection);
1780
1781 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1782 // set trigger (V0AND)
1783 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1784 if(!triggerAnalysis) return;
1785 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1786 }
1787 }
1788
1789 // centrality determination
1790 Float_t centralityF = -1;
1791 AliCentrality *esdCentrality = esdEvent->GetCentrality();
1792 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1793
1794
1795 // get reconstructed vertex
1796 //const AliESDVertex* vtxESD = 0;
1797 AliESDVertex* vtxESD = 0;
1798 if(GetAnalysisMode() == kTPCAnalysisMode) {
234cae2f 1799 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
c683985a 1800 }
1801 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
234cae2f 1802 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
c683985a 1803 }
1804 else {
234cae2f 1805 return;
c683985a 1806 }
1807
1808 if(!vtxESD) return;
1809
1810 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1811 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1812 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1813
1814 // check event cuts
1815 if(isEventOK && isEventTriggered) {
c683985a 1816 //
234cae2f 1817 // Dump the pt downscaled V0 into the tree
1818 //
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();
c683985a 1825
b55a470b 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);
1831
234cae2f 1832
1833 for (Int_t iv0=0; iv0<nV0s; iv0++){
1834
1835 AliESDv0 * v0 = esdEvent->GetV0(iv0);
1836 if (!v0) continue;
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;
82522b46 1843 if (esdFriend) {
1844 if (!esdFriend->TestSkipBit()){
dfd5cb6f 1845 Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
1846 if (v0->GetIndex(0)<ntracksFriend){
1847 friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
1848 }
1849 if (v0->GetIndex(1)<ntracksFriend){
1850 friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
1851 }
234cae2f 1852 }
1853 }
1a30582a 1854 if (track0->GetSign()<0) {
1855 track1 = esdEvent->GetTrack(v0->GetIndex(0));
1856 track0 = esdEvent->GetTrack(v0->GetIndex(1));
1857 }
1858
82522b46 1859 //
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;
1866 }
1867 }
1868 if (fFriendDownscaling<=0){
1869 if (((*fTreeSRedirector)<<"V0s").GetTree()){
1870 TTree * tree = ((*fTreeSRedirector)<<"V0s").GetTree();
1871 if (tree){
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;
1880 }
1881 }
1882 }
1883 }
1884
234cae2f 1885 //
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();
1892
1893 if(!fFillTree) return;
1894 if(!fTreeSRedirector) return;
1895 (*fTreeSRedirector)<<"V0s"<<
1a30582a 1896 "gid="<<gid<< // global id of event
1897 "isDownscaled="<<isDownscaled<< //
234cae2f 1898 "triggerClass="<<&triggerClass<< // trigger
1a30582a 1899 "Bz="<<bz<< //
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-
234cae2f 1905 "ntracks="<<ntracks<<
1906 "v0.="<<v0<<
1907 "kf.="<<&kfparticle<<
1a30582a 1908 "track0.="<<track0<< // track
234cae2f 1909 "track1.="<<track1<<
82522b46 1910 "friendTrack0.="<<friendTrackStore0<<
1911 "friendTrack1.="<<friendTrackStore1<<
234cae2f 1912 "centralityF="<<centralityF<<
1913 "\n";
234cae2f 1914 }
c683985a 1915 }
1916}
1917
1918//_____________________________________________________________________________
234cae2f 1919void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
c683985a 1920{
1921 //
1922 // Select real events with large TPC dEdx signal
1923 //
1924 if(!esdEvent) {
1925 AliDebug(AliLog::kError, "esdEvent not available");
1926 return;
1927 }
1928
1929 // get selection cuts
1930 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1931 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1932 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1933
1934 if(!evtCuts || !accCuts || !esdTrackCuts) {
1935 AliDebug(AliLog::kError, "cuts not available");
1936 return;
1937 }
1938
c683985a 1939
1940 // trigger
1941 Bool_t isEventTriggered = kTRUE;
1942 AliPhysicsSelection *physicsSelection = NULL;
1943 AliTriggerAnalysis* triggerAnalysis = NULL;
1944
1945 //
1946 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
c683985a 1947
1948 if(evtCuts->IsTriggerRequired())
1949 {
1950 // always MB
1951 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1952
1953 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1954 if(!physicsSelection) return;
1955
1956 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1957 // set trigger (V0AND)
1958 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1959 if(!triggerAnalysis) return;
1960 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1961 }
1962 }
1963
1964 // get reconstructed vertex
1965 AliESDVertex* vtxESD = 0;
1966 if(GetAnalysisMode() == kTPCAnalysisMode) {
234cae2f 1967 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
c683985a 1968 }
1969 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
234cae2f 1970 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
c683985a 1971 }
1972 else {
234cae2f 1973 return;
c683985a 1974 }
1975 if(!vtxESD) return;
1976
1977 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1978 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1979 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1980
1981
1982 // check event cuts
1983 if(isEventOK && isEventTriggered)
1984 {
1985 Double_t vert[3] = {0};
e690d4d0 1986 vert[0] = vtxESD->GetX();
1987 vert[1] = vtxESD->GetY();
1988 vert[2] = vtxESD->GetZ();
c683985a 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();
b55a470b 1994
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);
2000
c683985a 2001 // large dEdx
2002 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
2003 {
2004 AliESDtrack *track = esdEvent->GetTrack(iTrack);
2005 if(!track) continue;
234cae2f 2006 AliESDfriendTrack* friendTrack=NULL;
dfd5cb6f 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
2012 }
2013
c683985a 2014 if(track->Charge()==0) continue;
2015 if(!esdTrackCuts->AcceptTrack(track)) continue;
2016 if(!accCuts->AcceptTrack(track)) continue;
2017
2018 if(!IsHighDeDxParticle(track)) continue;
2019 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2020
2021 if(!fFillTree) return;
2022 if(!fTreeSRedirector) return;
1a30582a 2023 (*fTreeSRedirector)<<"dEdx"<< // high dEdx tree
2024 "gid="<<gid<< // global id
2025 "fileName.="<<&fCurrentFileName<< // file name
234cae2f 2026 "runNumber="<<runNumber<<
2027 "evtTimeStamp="<<evtTimeStamp<<
2028 "evtNumberInFile="<<evtNumberInFile<<
2029 "triggerClass="<<&triggerClass<< // trigger
2030 "Bz="<<bz<<
1a30582a 2031 "vtxESD.="<<vtxESD<< //
234cae2f 2032 "mult="<<mult<<
2033 "esdTrack.="<<track<<
2034 "friendTrack.="<<friendTrack<<
2035 "\n";
c683985a 2036 }
2037 }
2038}
2039
2040//_____________________________________________________________________________
2041Int_t AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
2042{
2043 //
2044 // Create KF particle in case the V0 fullfill selection criteria
2045 //
2046 // Selection criteria
2047 // 0. algorithm cut
2048 // 1. track cut
2049 // 3. chi2 cut
2050 // 4. rough mass cut
2051 // 5. Normalized pointing angle cut
2052 //
2053 const Double_t cutMass=0.2;
2054 const Double_t kSigmaDCACut=3;
2055 //
2056 // 0.) algo cut - accept only on the fly
2057 //
2058 if (v0->GetOnFlyStatus() ==kFALSE) return 0;
2059 //
2060 // 1.) track cut
2061 //
2062 AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
2063 AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
2064 /*
234cae2f 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";
2068 */
c683985a 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;
c683985a 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);
2077 //
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;
2080 //
2081 //
2082 // 3.) Chi2 cut
2083 //
2084 Double_t chi2KF = v0->GetKFInfo(2,2,2);
2085 if (chi2KF>25) return 0;
2086 //
2087 // 4.) Rough mass cut - 0.200 GeV
2088 //
2089 static Double_t masses[2]={-1};
2090 if (masses[0]<0){
2091 masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
2092 masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
2093 }
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;
2099 Int_t type=0;
2100 Int_t ptype=0;
2101 Double_t dmass=1;
2102 Int_t p1=0, p2=0;
2103 if (TMath::Abs(mass00-0)<cutMass) {
2104 massOK=kTRUE; type+=1;
2105 if (TMath::Abs(mass00-0)<dmass) {
2106 ptype=1;
2107 dmass=TMath::Abs(mass00-0);
2108 p1=0; p2=0;
2109 }
2110 }
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]);
2115 ptype=2;
2116 p1=2; p2=4;
2117 }
2118 }
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]);
2123 ptype=4;
2124 p1=4; p2=2;
2125 }
2126 }
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]);
2131 ptype=8;
2132 p1=2; p2=2;
2133 }
2134 }
2135 if (type==0) return 0;
2136 //
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();
2143 }
2144 //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
2145 //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
2146 //
2147 AliKFParticle kfp1( *paramP, spdg[p1] );
2148 AliKFParticle kfp2( *paramN, -1 *spdg[p2] );
2149 AliKFParticle V0KF;
2150 (V0KF)+=kfp1;
2151 (V0KF)+=kfp2;
2152 kfparticle=V0KF;
2153 //
2154 // Pointing angle
2155 //
2156 Double_t errPhi = V0KF.GetErrPhi();
2157 Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
2158 if (pointAngle/errPhi>10) return 0;
2159 //
2160 return ptype;
2161}
2162
2163//_____________________________________________________________________________
2164Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)
2165{
2166 //
2167 // Downscale randomly low pt V0
2168 //
2169 //return kFALSE;
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;
2174 //
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
2180 }
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;
2183 return kFALSE;
2184
2185 /*
234cae2f 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();
2191 his1->Fill(rnd);
2192 if (!isDownscaled) his2->Fill(rnd);
2193 }}
c683985a 2194
234cae2f 2195*/
c683985a 2196
2197}
2198
2199
2200
2201//_____________________________________________________________________________
2202Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
2203{
234cae2f 2204 // Constrain TPC inner params constrained
2205 //
2206 if(!tpcInnerC) return kFALSE;
2207 if(!vtx) return kFALSE;
c683985a 2208
234cae2f 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;
c683985a 2214
2215
234cae2f 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;
c683985a 2221
234cae2f 2222 if(!tpcInnerC->Update(p,c)) return kFALSE;
c683985a 2223
2224 return kTRUE;
2225}
2226
2227//_____________________________________________________________________________
2228Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
2229{
234cae2f 2230 // Constrain TPC inner params constrained
2231 //
2232 if(!trackInnerC) return kFALSE;
2233 if(!vtx) return kFALSE;
c683985a 2234
234cae2f 2235 const Double_t kRadius = 2.8;
2236 const Double_t kMaxStep = 1.0;
c683985a 2237
234cae2f 2238 Double_t dz[2],cov[3];
c683985a 2239
234cae2f 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;
c683985a 2243
234cae2f 2244 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
2245 if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
c683985a 2246
234cae2f 2247 //
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;
c683985a 2253
234cae2f 2254 if(!trackInnerC->Update(p,c)) return kFALSE;
c683985a 2255
2256 return kTRUE;
2257}
2258
2259
2260//_____________________________________________________________________________
2261TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack)
2262{
2263 if(!particle) return NULL;
2264 if(!stack) return NULL;
2265
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);
2271
234cae2f 2272 return mother;
c683985a 2273}
2274
2275//_____________________________________________________________________________
34e8fe1d 2276Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(Int_t label, AliStack *const stack)
c683985a 2277{
2278 Bool_t isFromConversion = kFALSE;
2279
2280 if(stack) {
2281 Int_t mcStackSize=stack->GetNtrack();
2282 if (label>=mcStackSize) return kFALSE;
2283 TParticle* particle = stack->Particle(label);
2284 if (!particle) return kFALSE;
2285
2286 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2287 {
234cae2f 2288 Int_t mech = particle->GetUniqueID(); // production mechanism
2289 Bool_t isPrim = stack->IsPhysicalPrimary(label);
c683985a 2290
234cae2f 2291 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2292 if (motherLabel>=mcStackSize) return kFALSE;
2293 TParticle* mother = stack->Particle(motherLabel);
2294 if(mother) {
2295 Int_t motherPdg = mother->GetPdgCode();
c683985a 2296
234cae2f 2297 if(!isPrim && mech==5 && motherPdg==kGamma) {
2298 isFromConversion=kTRUE;
2299 }
2300 }
c683985a 2301 }
2302 }
2303
2304 return isFromConversion;
2305}
2306
2307//_____________________________________________________________________________
34e8fe1d 2308Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(Int_t label, AliStack *const stack)
c683985a 2309{
2310 Bool_t isFromMaterial = kFALSE;
2311
2312 if(stack) {
2313 Int_t mcStackSize=stack->GetNtrack();
2314 if (label>=mcStackSize) return kFALSE;
2315 TParticle* particle = stack->Particle(label);
2316 if (!particle) return kFALSE;
2317
2318 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2319 {
234cae2f 2320 Int_t mech = particle->GetUniqueID(); // production mechanism
2321 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2322
2323 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2324 if (motherLabel>=mcStackSize) return kFALSE;
2325 TParticle* mother = stack->Particle(motherLabel);
2326 if(mother) {
2327 if(!isPrim && mech==13) {
2328 isFromMaterial=kTRUE;
2329 }
2330 }
2331 }
c683985a 2332 }
2333
2334 return isFromMaterial;
2335}
2336
2337//_____________________________________________________________________________
34e8fe1d 2338Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(Int_t label, AliStack *const stack)
c683985a 2339{
2340 Bool_t isFromStrangeness = kFALSE;
2341
2342 if(stack) {
2343 Int_t mcStackSize=stack->GetNtrack();
2344 if (label>=mcStackSize) return kFALSE;
2345 TParticle* particle = stack->Particle(label);
2346 if (!particle) return kFALSE;
2347
2348 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2349 {
234cae2f 2350 Int_t mech = particle->GetUniqueID(); // production mechanism
2351 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2352
2353 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2354 if (motherLabel>=mcStackSize) return kFALSE;
2355 TParticle* mother = stack->Particle(motherLabel);
2356 if(mother) {
2357 Int_t motherPdg = mother->GetPdgCode();
2358
2359 // K+-, lambda, antilambda, K0s decays
2360 if(!isPrim && mech==4 &&
2361 (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
2362 {
2363 isFromStrangeness = kTRUE;
2364 }
2365 }
c683985a 2366 }
2367 }
2368
2369 return isFromStrangeness;
2370}
2371
2372
2373//_____________________________________________________________________________
2374void AliAnalysisTaskFilteredTree::FinishTaskOutput()
2375{
2376 //
2377 // Called one at the end
2378 // locally on working node
2379 //
c683985a 2380 Bool_t deleteTrees=kTRUE;
2381 if ((AliAnalysisManager::GetAnalysisManager()))
2382 {
2383 if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
234cae2f 2384 AliAnalysisManager::kProofAnalysis)
c683985a 2385 deleteTrees=kFALSE;
2386 }
2387 if (deleteTrees) delete fTreeSRedirector;
2388 fTreeSRedirector=NULL;
c683985a 2389}
2390
2391//_____________________________________________________________________________
2392void AliAnalysisTaskFilteredTree::Terminate(Option_t *)
2393{
1a30582a 2394 //
c683985a 2395 // Called one at the end
1a30582a 2396 //
c683985a 2397}
2398
2399//_____________________________________________________________________________
2400Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
2401{
2402 //
2403 // calculate mc event true track multiplicity
2404 //
2405 if(!mcEvent) return 0;
2406
2407 AliStack* stack = 0;
2408 Int_t mult = 0;
2409
2410 // MC particle stack
2411 stack = mcEvent->Stack();
2412 if (!stack) return 0;
2413
2414 //
2415 //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
2416 //
2417
2418 Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
2419 if(!isEventOK) return 0;
2420
2421 Int_t nPart = stack->GetNtrack();
2422 for (Int_t iMc = 0; iMc < nPart; ++iMc)
2423 {
234cae2f 2424 TParticle* particle = stack->Particle(iMc);
2425 if (!particle)
2426 continue;
2427
2428 // only charged particles
2429 if(!particle->GetPDG()) continue;
2430 Double_t charge = particle->GetPDG()->Charge()/3.;
2431 if (TMath::Abs(charge) < 0.001)
2432 continue;
2433
2434 // physical primary
2435 Bool_t prim = stack->IsPhysicalPrimary(iMc);
2436 if(!prim) continue;
2437
2438 // checked accepted without pt cut
2439 //if(accCuts->AcceptTrack(particle))
2440 if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
2441 {
2442 mult++;
2443 }
c683985a 2444 }
2445
234cae2f 2446 return mult;
c683985a 2447}
2448
2449//_____________________________________________________________________________
34e8fe1d 2450void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC)
c683985a 2451{
234cae2f 2452 //
2453 // Fill pT relative resolution histograms for
2454 // TPC only, TPC only constrained to vertex and TPC+ITS tracking
2455 //
2456 if(!ptrack) return;
2457 if(!ptpcInnerC) return;
2458
2459 const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();
2460 if(!innerParam) return;
2461
2462 Float_t dxy, dz;
2463 ptrack->GetImpactParameters(dxy,dz);
2464
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)) &&
2480 abs(dz)<2.0 &&
2481 abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )
2482 {
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()));
2486 }
c683985a 2487
234cae2f 2488 // TPC primary tracks
2489 // and TPC constrained primary tracks
c683985a 2490
234cae2f 2491 AliExternalTrackParam *ptpcInner = (AliExternalTrackParam *) ptrack->GetTPCInnerParam();
2492 if(!ptpcInner) return;
c683985a 2493
2494
234cae2f 2495 Float_t dxyTPC, dzTPC;
2496 ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);
c683985a 2497
234cae2f 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 &&
2506 abs(dzTPC)<3.2 &&
2507 abs(dxyTPC)<2.4 )
2508 {
2509 // TPC only
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()));
2513
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()));
2518 }
c683985a 2519}
4a00e62d 2520
2521
2522void AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend){
2523 //
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
2529 //
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
2535 //
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
2542 //
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();
2547 //
2548 // 0.) Get parameters in reference radius TPC Inner wall
2549 //
2550 //
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];
2566 }
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();;
2571 }
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];
2581 }
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();;
2586 }
2587 }
2588 //
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
2593 //
2594 //
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++;
2609 //
2610 // 2.) find clostest TPCtrack
2611 // a.) all tracks
2612 Double_t minChi2All=10000000;
2613 Double_t minChi2TPC=10000000;
2614 Double_t minChi2TPCITS=10000000;
2615 Int_t indexAll=-1;
2616 Int_t indexTPC=-1;
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;
2626 // fast checks
2627 //
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;
2633 ncandidates0++;
2634 //
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;
2642 ncandidates1++;
2643 if (chi2<minChi2All){
2644 minChi2All=chi2;
2645 indexAll=iTrack1;
2646 }
2647 if (chi2<minChi2TPC && track1->IsOn(AliVTrack::kITSin)==0){
2648 minChi2TPC=chi2;
2649 indexTPC=iTrack1;
2650 itsAtTPC=outerITS;
2651 }
2652 if (chi2<minChi2TPCITS && track1->IsOn(AliVTrack::kITSin)){
2653 minChi2TPCITS=chi2;
2654 indexTPCITS=iTrack1;
2655 itsAtITSTPC=outerITS;
2656 }
2657 }
2658 //
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<<
2668 //
2669 "chi2All="<<minChi2All<< // chi2 of closest tracks
2670 "chi2TPC="<<minChi2TPC<<
2671 "chi2TPCITS="<<minChi2TPCITS<<
2672 //
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
2677 //
2678 "itsAtTPC.="<<&itsAtTPC<< // ITS track parameters at the TPC alone track frame
2679 "itsAtITSTPC.="<<&itsAtITSTPC<< // ITS track parameters at the TPC combeined track frame
2680 "\n";
2681 }
2682}
2683
b4990577 2684void AliAnalysisTaskFilteredTree::ProcessTrackMatch(AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/){
4a00e62d 2685/*
b4990577 2686 Use TPC standalone, ITS standalone and combined tracks to categorize/resp. recover track information.
2687
2688 Track categories:
2689 -TPC+ITS
2690 -TPC only
2691 -ITS only
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)
2697 Info:
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
2702 Information matrix:
2703 -matrix closest tracks from each categories
2704 -characterization - chi2, index,chi2, overlap ratio
2705
2706 Decissions:
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
4a00e62d 2712
b4990577 2713 1.) if (TPConly && !(TPC+ITS) && ITSonly match ) {
2714 if (!kink) TPCOnly.addoptITS
2715 if (kink) TPConly sign
2716 }
4a00e62d 2717
b4990577 2718 2.) Overlap tracks - Refit with doUnfold
4a00e62d 2719
2720*/
b4990577 2721}