]> 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
Added macros and scripts for ZDC QA automatic trending (M.Leoncino)
[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;
371
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;
392 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack0 = esdFriend->GetTrack(itrack0);} //this guy can be NULL
234cae2f 393
394 for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
395 AliESDtrack *track1 = event->GetTrack(itrack1);
396 if (!track1) continue;
397 if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
398 if (track1->GetKinkIndex(0)>0) continue;
399 if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue;
400 if (track1->Pt() < kMinPt) continue;
401 if (track1->GetTPCncls()<kMinNcl) continue;
402 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
403 if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
404 //track1->GetImpactParametersTPC(dcaTPC,covTPC);
405 // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
c683985a 406 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
234cae2f 407 //
408 const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
409 //
410 Bool_t isPair=kTRUE;
411 for (Int_t ipar=0; ipar<5; ipar++){
412 if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
413 if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
414 }
415 if (!isPair) continue;
416 if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
417 //delta with correct sign
418 /*
82522b46 419 TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
420 TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
421 TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
422 */
234cae2f 423 if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
424 if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
425 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
426 if (!isPair) continue;
427 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
428 Int_t eventNumber = event->GetEventNumberInFile();
234cae2f 429 //
430 //
431 Int_t ntracksSPD = vertexSPD->GetNContributors();
432 Int_t ntracksTPC = vertexTPC->GetNContributors();
433 Int_t runNumber = event->GetRunNumber();
434 Int_t timeStamp = event->GetTimeStamp();
435 ULong64_t triggerMask = event->GetTriggerMask();
436 Float_t magField = event->GetMagneticField();
437 TObjString triggerClass = event->GetFiredTriggerClasses().Data();
438
b55a470b 439 // Global event id calculation using orbitID, bunchCrossingID and periodID
440 ULong64_t orbitID = (ULong64_t)event->GetOrbitNumber();
441 ULong64_t bunchCrossID = (ULong64_t)event->GetBunchCrossNumber();
442 ULong64_t periodID = (ULong64_t)event->GetPeriodNumber();
443 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
444
c683985a 445
234cae2f 446 AliESDfriendTrack* friendTrack1=NULL;
447 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack1 = esdFriend->GetTrack(itrack1);} //this guy can be NULL
82522b46 448 //
449 AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
450 AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
451 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
452 if (gRandom->Rndm()>1./fFriendDownscaling){
453 friendTrackStore0 = 0;
454 friendTrackStore1 = 0;
455 }
456 }
457 if (fFriendDownscaling<=0){
458 if (((*fTreeSRedirector)<<"CosmicPairs").GetTree()){
459 TTree * tree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
460 if (tree){
461 Double_t sizeAll=tree->GetZipBytes();
462 TBranch * br= tree->GetBranch("friendTrack0.fPoints");
463 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
464 br= tree->GetBranch("friendTrack0.fCalibContainer");
465 if (br) sizeFriend+=br->GetZipBytes();
466 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
467 friendTrackStore0=0;
468 friendTrackStore1=0;
469 }
470 }
471 }
472 }
c683985a 473 if(!fFillTree) return;
474 if(!fTreeSRedirector) return;
234cae2f 475 (*fTreeSRedirector)<<"CosmicPairs"<<
1a30582a 476 "gid="<<gid<< // global id of track
477 "fileName.="<<&fCurrentFileName<< // file name
478 "runNumber="<<runNumber<< // run number
479 "evtTimeStamp="<<timeStamp<< // time stamp of event
480 "evtNumberInFile="<<eventNumber<< // event number
481 "trigger="<<triggerMask<< // trigger mask
482 "triggerClass="<<&triggerClass<< // trigger class
483 "Bz="<<magField<< // magnetic field
234cae2f 484 //
1a30582a 485 "multSPD="<<ntracksSPD<< // event ultiplicity
486 "multTPC="<<ntracksTPC<< //
487 "vertSPD.="<<vertexSPD<< // primary vertex -SPD
488 "vertTPC.="<<vertexTPC<< // primary vertex -TPC
489 "t0.="<<track0<< // first half of comsic trak
490 "t1.="<<track1<< // second half of cosmic track
491 "friendTrack0.="<<friendTrackStore0<< // friend information first track + points
492 "friendTrack1.="<<friendTrackStore1<< // frined information first track + points
234cae2f 493 "\n";
234cae2f 494 }
234cae2f 495 }
c683985a 496}
497
498
499//_____________________________________________________________________________
500void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
501{
502 //
503 // Select real events with high-pT tracks
504 //
c683985a 505
506 // get selection cuts
507 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
508 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
509 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
510
511 if(!evtCuts || !accCuts || !esdTrackCuts) {
512 Printf("ERROR cuts not available");
513 return;
514 }
515
516 // trigger selection
517 Bool_t isEventTriggered = kTRUE;
518 AliPhysicsSelection *physicsSelection = NULL;
519 AliTriggerAnalysis* triggerAnalysis = NULL;
520
521 //
522 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
c683985a 523 // trigger
524 if(evtCuts->IsTriggerRequired())
525 {
526 // always MB
527 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
528
529 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
530 if(!physicsSelection) return;
531 //SetPhysicsTriggerSelection(physicsSelection);
532
533 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
534 // set trigger (V0AND)
535 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
536 if(!triggerAnalysis) return;
537 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
538 }
539 }
540
541 // centrality determination
542 Float_t centralityF = -1;
543 AliCentrality *esdCentrality = esdEvent->GetCentrality();
544 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
545
546 // use MC information
547 AliHeader* header = 0;
548 AliGenEventHeader* genHeader = 0;
549 AliStack* stack = 0;
550 TArrayF vtxMC(3);
551
552 Int_t multMCTrueTracks = 0;
553 if(mcEvent)
554 {
555 // get MC event header
556 header = mcEvent->Header();
557 if (!header) {
558 AliDebug(AliLog::kError, "Header not available");
559 return;
560 }
561 // MC particle stack
562 stack = mcEvent->Stack();
563 if (!stack) {
564 AliDebug(AliLog::kError, "Stack not available");
565 return;
566 }
567
568 // get MC vertex
569 genHeader = header->GenEventHeader();
570 if (!genHeader) {
571 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
572 return;
573 }
574 genHeader->PrimaryVertex(vtxMC);
575
576 // multipliticy of all MC primary tracks
577 // in Zv, pt and eta ranges)
578 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
579 } // end bUseMC
580
581 // get reconstructed vertex
582 //const AliESDVertex* vtxESD = 0;
583 AliESDVertex* vtxESD = 0;
584 if(GetAnalysisMode() == kTPCAnalysisMode) {
234cae2f 585 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
c683985a 586 }
587 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
234cae2f 588 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
c683985a 589 }
590 else {
234cae2f 591 return;
c683985a 592 }
593
594 AliESDVertex* vtxTPC = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
595 AliESDVertex* vtxSPD = (AliESDVertex*)esdEvent->GetPrimaryVertexSPD();
596
597 if(!vtxESD) return;
598 if(!vtxTPC) return;
599 if(!vtxSPD) return;
600
601 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
602 //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZv());
603 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
604 Int_t ntracks = esdEvent->GetNumberOfTracks();
605
606 // check event cuts
607 if(isEventOK && isEventTriggered)
608 {
609
610 //
611 // get IR information
612 //
613 AliESDHeader *esdHeader = 0;
614 esdHeader = esdEvent->GetHeader();
615 if(!esdHeader) return;
616 //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
617 //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
618
619 // Use when Peter commit the changes in the header
620 Int_t ir1 = 0;
621 Int_t ir2 = 0;
622
623 //
07aef7d1 624 //Double_t vert[3] = {0};
625 //vert[0] = vtxESD->GetXv();
626 //vert[1] = vtxESD->GetYv();
627 //vert[2] = vtxESD->GetZv();
c683985a 628 Int_t mult = vtxESD->GetNContributors();
629 Int_t multSPD = vtxSPD->GetNContributors();
630 Int_t multTPC = vtxTPC->GetNContributors();
631
632 Float_t bz = esdEvent->GetMagneticField();
633 Int_t runNumber = esdEvent->GetRunNumber();
634 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
635 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
b55a470b 636
637 // Global event id calculation using orbitID, bunchCrossingID and periodID
638 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
639 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
640 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
641 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
642
c683985a 643
644 // high pT tracks
645 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
646 {
647 AliESDtrack *track = esdEvent->GetTrack(iTrack);
648 if(!track) continue;
649 if(track->Charge()==0) continue;
650 if(!esdTrackCuts->AcceptTrack(track)) continue;
651 if(!accCuts->AcceptTrack(track)) continue;
234cae2f 652
c683985a 653 // downscale low-pT tracks
654 Double_t scalempt= TMath::Min(track->Pt(),10.);
655 Double_t downscaleF = gRandom->Rndm();
656 downscaleF *= fLowPtTrackDownscaligF;
657 if(TMath::Exp(2*scalempt)<downscaleF) continue;
658 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
659
660 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
661 if (!tpcInner) continue;
662 // transform to the track reference frame
663 Bool_t isOK = kFALSE;
664 isOK = tpcInner->Rotate(track->GetAlpha());
665 isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
666 if(!isOK) continue;
667
668 // Dump to the tree
669 // vertex
670 // TPC-ITS tracks
671 //
672 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
c683985a 673 if(!fFillTree) return;
674 if(!fTreeSRedirector) return;
675 (*fTreeSRedirector)<<"highPt"<<
b55a470b 676 "gid="<<gid<<
1a30582a 677 "fileName.="<<&fCurrentFileName<<
c683985a 678 "runNumber="<<runNumber<<
679 "evtTimeStamp="<<evtTimeStamp<<
680 "evtNumberInFile="<<evtNumberInFile<<
234cae2f 681 "triggerClass="<<&triggerClass<< // trigger
c683985a 682 "Bz="<<bz<< // magnetic field
683 "vtxESD.="<<vtxESD<<
234cae2f 684 "ntracksESD="<<ntracks<< // number of tracks in the ESD
234cae2f 685 "IRtot="<<ir1<< // interaction record history info
686 "IRint2="<<ir2<<
c683985a 687 "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
688 "multSPD="<<multSPD<< // multiplicity of tracks pointing to the SPD primary vertex
689 "multTPC="<<multTPC<< // multiplicity of tracks pointing to the TPC primary vertex
690 "esdTrack.="<<track<<
691 "centralityF="<<centralityF<<
692 "\n";
693 }
694 }
234cae2f 695
c683985a 696}
697
698
699//_____________________________________________________________________________
234cae2f 700void AliAnalysisTaskFilteredTree::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const esdFriend)
c683985a 701{
702 //
82522b46 703 // Process laser events -> dump tracks and clusters to the special tree
c683985a 704 //
1236ed2e 705 const Double_t kMinPt = 5;
234cae2f 706 if(!fFillTree) return;
707 if(!fTreeSRedirector) return;
c683985a 708 const AliESDHeader* esdHeader = esdEvent->GetHeader();
1a30582a 709 if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) {
c683985a 710 Int_t countLaserTracks = 0;
234cae2f 711 Int_t runNumber = esdEvent->GetRunNumber();
712 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
713 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
714 Float_t bz = esdEvent->GetMagneticField();
82522b46 715 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
b55a470b 716 // Global event id calculation using orbitID, bunchCrossingID and periodID
717 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
718 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
719 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
720 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1a30582a 721 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
c683985a 722 AliESDtrack *track = esdEvent->GetTrack(iTrack);
723 if(!track) continue;
1a30582a 724 if(track->GetTPCInnerParam()) countLaserTracks++;
234cae2f 725 AliESDfriendTrack* friendTrack=NULL;
1236ed2e 726 // suppress beam background and CE random reacks
15120578 727 if (track->GetInnerParam()->Pt()<kMinPt) continue;
1236ed2e 728 Bool_t skipTrack=gRandom->Rndm()>1/(1+TMath::Abs(fFriendDownscaling));
729 if (skipTrack) continue;
1a30582a 730 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
c683985a 731 (*fTreeSRedirector)<<"Laser"<<
82522b46 732 "gid="<<gid<< // global identifier of event
1a30582a 733 "fileName.="<<&fCurrentFileName<< //
c683985a 734 "runNumber="<<runNumber<<
735 "evtTimeStamp="<<evtTimeStamp<<
736 "evtNumberInFile="<<evtNumberInFile<<
82522b46 737 "triggerClass="<<&triggerClass<< // trigger
738 "Bz="<<bz<< // magnetic field
739 "multTPCtracks="<<countLaserTracks<< // multiplicity of tracks
1236ed2e 740 "track.="<<track<< // track parameters
82522b46 741 "friendTrack.="<<friendTrack<< // friend track information
c683985a 742 "\n";
743 }
744 }
745}
746
747//_____________________________________________________________________________
748void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
749{
750 //
751 // Select real events with high-pT tracks
752 // Calculate and stor in the output tree:
753 // TPC constrained tracks
754 // InnerParams constrained tracks
755 // TPC-ITS tracks
756 // ITSout-InnerParams tracks
757 // chi2 distance between TPC constrained and TPC-ITS tracks
758 // chi2 distance between TPC refitted constrained and TPC-ITS tracks
759 // chi2 distance between ITSout and InnerParams tracks
760 // MC information:
761 // track references at ITSin, TPCin; InnerParam at first TPC track reference,
762 // particle ID, mother ID, production mechanism ...
763 //
c683985a 764 // get selection cuts
765 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
766 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
767 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
768
769 if(!evtCuts || !accCuts || !esdTrackCuts) {
770 AliDebug(AliLog::kError, "cuts not available");
771 return;
772 }
773
774 // trigger selection
775 Bool_t isEventTriggered = kTRUE;
776 AliPhysicsSelection *physicsSelection = NULL;
777 AliTriggerAnalysis* triggerAnalysis = NULL;
778
779 //
780 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
234cae2f 781
c683985a 782
783 // trigger
784 if(evtCuts->IsTriggerRequired())
785 {
786 // always MB
787 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
788
789 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
790 if(!physicsSelection) {AliInfo("no physics selection"); return;}
791 //SetPhysicsTriggerSelection(physicsSelection);
792
793 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
794 // set trigger (V0AND)
795 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
796 if(!triggerAnalysis) {AliInfo("no trigger analysis");return;}
797 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
798 }
799 }
800
801 // centrality determination
802 Float_t centralityF = -1;
803 AliCentrality *esdCentrality = esdEvent->GetCentrality();
804 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
805
806 // use MC information
807 AliHeader* header = 0;
808 AliGenEventHeader* genHeader = 0;
809 AliStack* stack = 0;
810 TArrayF vtxMC(3);
811 Int_t mcStackSize=0;
812
813 Int_t multMCTrueTracks = 0;
814 if(mcEvent)
815 {
816 // get MC event header
817 header = mcEvent->Header();
818 if (!header) {
819 AliDebug(AliLog::kError, "Header not available");
820 return;
821 }
822 // MC particle stack
823 stack = mcEvent->Stack();
824 if (!stack) {
825 AliDebug(AliLog::kError, "Stack not available");
826 return;
827 }
828 mcStackSize=stack->GetNtrack();
829
830 // get MC vertex
831 genHeader = header->GenEventHeader();
832 if (!genHeader) {
833 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
834 return;
835 }
836 genHeader->PrimaryVertex(vtxMC);
837
838 // multipliticy of all MC primary tracks
839 // in Zv, pt and eta ranges)
840 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
841
842 } // end bUseMC
843
844 // get reconstructed vertex
845 //const AliESDVertex* vtxESD = 0;
846 AliESDVertex* vtxESD = 0;
847 if(GetAnalysisMode() == kTPCAnalysisMode) {
234cae2f 848 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
c683985a 849 }
850 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
234cae2f 851 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
c683985a 852 }
853 else {
854 AliInfo("no ESD vertex");
234cae2f 855 return;
c683985a 856 }
857
858 if(!vtxESD) return;
859
860 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
234cae2f 861
c683985a 862 //
863 // Vertex info comparison and track multiplicity
864 //
865 AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
866 AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
867 Int_t contSPD = vertexSPD->GetNContributors();
868 Int_t contTPC = vertexTPC->GetNContributors();
869 TVectorD vertexPosTPC(3), vertexPosSPD(3);
870 vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
871 vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
872 Int_t ntracksTPC=0;
873 Int_t ntracksITS=0;
874 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
875 AliESDtrack *track = esdEvent->GetTrack(iTrack);
876 if(!track) continue;
877 if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
878 if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
879 }
234cae2f 880
c683985a 881 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
882 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
883
884
885 // check event cuts
886 if(isEventOK && isEventTriggered)
887 {
888 //
889 // get IR information
890 //
891 AliESDHeader *esdHeader = 0;
892 esdHeader = esdEvent->GetHeader();
893 if(!esdHeader) {AliInfo("no esdHeader");return;}
894 //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
895 //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
896 //
897 Int_t ir1 = 0;
898 Int_t ir2 = 0;
899
900 //
901 Double_t vert[3] = {0};
902 vert[0] = vtxESD->GetXv();
903 vert[1] = vtxESD->GetYv();
904 vert[2] = vtxESD->GetZv();
905 Int_t mult = vtxESD->GetNContributors();
906 Float_t bz = esdEvent->GetMagneticField();
907 Int_t runNumber = esdEvent->GetRunNumber();
908 Int_t evtTimeStamp = esdEvent->GetTimeStamp();
909 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
910
2e16763c 911 Int_t numberOfTracks=esdEvent->GetNumberOfTracks();
c683985a 912 // high pT tracks
2e16763c 913 for (Int_t iTrack = 0; iTrack < numberOfTracks; iTrack++)
c683985a 914 {
915 AliESDtrack *track = esdEvent->GetTrack(iTrack);
234cae2f 916 AliESDfriendTrack* friendTrack=NULL;
2e16763c 917 Int_t numberOfFriendTracks=0;
918 if (esdFriend) numberOfFriendTracks=esdFriend->GetNumberOfTracks();
919 if (esdFriend && iTrack<numberOfFriendTracks) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
c683985a 920 if(!track) continue;
921 if(track->Charge()==0) continue;
922 if(!esdTrackCuts->AcceptTrack(track)) continue;
923 if(!accCuts->AcceptTrack(track)) continue;
234cae2f 924
c683985a 925 // downscale low-pT tracks
926 Double_t scalempt= TMath::Min(track->Pt(),10.);
927 Double_t downscaleF = gRandom->Rndm();
928 downscaleF *= fLowPtTrackDownscaligF;
929 if(TMath::Exp(2*scalempt)<downscaleF) continue;
930 //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
931
932 // Dump to the tree
933 // vertex
934 // TPC constrained tracks
935 // InnerParams constrained tracks
936 // TPC-ITS tracks
937 // ITSout-InnerParams tracks
938 // chi2 distance between TPC constrained and TPC-ITS tracks
939 // chi2 distance between TPC refitted constrained and TPC-ITS tracks
940 // chi2 distance between ITSout and InnerParams tracks
941 // MC information
234cae2f 942
c683985a 943 Double_t x[3]; track->GetXYZ(x);
944 Double_t b[3]; AliTracker::GetBxByBz(x,b);
945
946 //
947 // Transform TPC inner params to track reference frame
948 //
949 Bool_t isOKtpcInner = kFALSE;
950 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
951 if (tpcInner) {
952 // transform to the track reference frame
953 isOKtpcInner = tpcInner->Rotate(track->GetAlpha());
954 isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
955 }
956
957 //
958 // Constrain TPC inner to vertex
959 // clone TPCinner has to be deleted
960 //
961 Bool_t isOKtpcInnerC = kFALSE;
962 AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
963 if (tpcInnerC) {
964 isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);
965 isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());
966 isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
967 }
968
969 //
970 // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex
971 // Clone track InnerParams has to be deleted
972 //
973 Bool_t isOKtrackInnerC = kFALSE;
974 AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));
975 if (trackInnerC) {
976 isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);
977 isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());
978 isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
979 }
234cae2f 980
c683985a 981 //
982 // calculate chi2 between vi and vj vectors
983 // with covi and covj covariance matrices
984 // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
985 //
986 TMatrixD deltaT(5,1), deltaTtrackC(5,1);
987 TMatrixD delta(1,5), deltatrackC(1,5);
988 TMatrixD covarM(5,5), covarMtrackC(5,5);
989 TMatrixD chi2(1,1);
990 TMatrixD chi2trackC(1,1);
991
992 if(isOKtpcInnerC && isOKtrackInnerC)
993 {
994 for (Int_t ipar=0; ipar<5; ipar++) {
995 deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
234cae2f 996 delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
c683985a 997
998 deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
234cae2f 999 deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
c683985a 1000
1001 for (Int_t jpar=0; jpar<5; jpar++) {
234cae2f 1002 Int_t index=track->GetIndex(ipar,jpar);
1003 covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
1004 covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];
c683985a 1005 }
1006 }
1007
1008 // chi2 distance TPC constrained and TPC+ITS
1009 TMatrixD covarMInv = covarM.Invert();
1010 TMatrixD mat2 = covarMInv*deltaT;
1011 chi2 = delta*mat2;
1012 //chi2.Print();
1013
1014 // chi2 distance TPC refitted constrained and TPC+ITS
1015 TMatrixD covarMInvtrackC = covarMtrackC.Invert();
1016 TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;
1017 chi2trackC = deltatrackC*mat2trackC;
1018 //chi2trackC.Print();
1019 }
1020
1021
1022 //
1023 // Propagate ITSout to TPC inner wall
1024 // and calculate chi2 distance to track (InnerParams)
1025 //
1026 const Double_t kTPCRadius=85;
1027 const Double_t kStep=3;
1028
1029 // clone track InnerParams has to be deleted
1030 Bool_t isOKtrackInnerC2 = kFALSE;
1031 AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));
1032 if (trackInnerC2) {
1033 isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);
1034 }
1035
1036 Bool_t isOKouterITSc = kFALSE;
1037 AliExternalTrackParam *outerITSc = NULL;
1038 TMatrixD chi2OuterITS(1,1);
1039
234cae2f 1040 if(esdFriend && !esdFriend->TestSkipBit())
c683985a 1041 {
1042 // propagate ITSout to TPC inner wall
c683985a 1043 if(friendTrack)
234cae2f 1044 {
1045
1046 outerITSc = NULL;
1047 if (friendTrack->GetITSOut()) outerITSc = new AliExternalTrackParam(*(friendTrack->GetITSOut()));
c683985a 1048 if(outerITSc)
234cae2f 1049 {
c683985a 1050 isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);
1051 isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());
1052 isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());
1053
234cae2f 1054 //
c683985a 1055 // calculate chi2 between outerITS and innerParams
234cae2f 1056 // cov without z-coordinate at the moment
1057 //
c683985a 1058 TMatrixD deltaTouterITS(4,1);
1059 TMatrixD deltaouterITS(1,4);
1060 TMatrixD covarMouterITS(4,4);
1061
1062 if(isOKtrackInnerC2 && isOKouterITSc) {
234cae2f 1063 Int_t kipar = 0;
1064 Int_t kjpar = 0;
c683985a 1065 for (Int_t ipar=0; ipar<5; ipar++) {
234cae2f 1066 if(ipar!=1) {
c683985a 1067 deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
234cae2f 1068 deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
1069 }
c683985a 1070
1071 kjpar=0;
1072 for (Int_t jpar=0; jpar<5; jpar++) {
234cae2f 1073 Int_t index=outerITSc->GetIndex(ipar,jpar);
1074 if(ipar !=1 || jpar!=1) {
1075 covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];
1076 }
c683985a 1077 if(jpar!=1) kjpar++;
234cae2f 1078 }
1079 if(ipar!=1) kipar++;
1080 }
c683985a 1081
1082 // chi2 distance ITSout and InnerParams
1083 TMatrixD covarMInvouterITS = covarMouterITS.Invert();
1084 TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;
1085 chi2OuterITS = deltaouterITS*mat2outerITS;
1086 //chi2OuterITS.Print();
234cae2f 1087 }
c683985a 1088 }
1089 }
1090 }
1091
1092 //
1093 // MC info
1094 //
1095 TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;
1096 TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;
1097 Int_t mech=-1, mechTPC=-1, mechITS=-1;
1098 Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;
1099 Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;
1100 Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;
1101 Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;
1102
1103 AliTrackReference *refTPCIn = NULL;
1104 AliTrackReference *refTPCOut = NULL;
1105 AliTrackReference *refITS = NULL;
1106 AliTrackReference *refTRD = NULL;
1107 AliTrackReference *refTOF = NULL;
1108 AliTrackReference *refEMCAL = NULL;
1109 AliTrackReference *refPHOS = NULL;
1110 Int_t nrefTPC=0, nrefTRD=0, nrefTOF=0, nrefITS=0, nrefEMCAL=0, nrefPHOS=0;
1111
1112 Bool_t isOKtrackInnerC3 = kFALSE;
1113 AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
2e16763c 1114 if(mcEvent && stack)
c683985a 1115 {
2e16763c 1116 do //artificial loop (once) to make the continue statements jump out of the MC part
234cae2f 1117 {
2e16763c 1118 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1119 //
1120 // global track
1121 //
1122 Int_t label = TMath::Abs(track->GetLabel());
1123 if (label >= mcStackSize) continue;
1124 particle = stack->Particle(label);
1125 if (!particle) continue;
1126 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)
1127 {
1128 particleMother = GetMother(particle,stack);
1129 mech = particle->GetUniqueID();
1130 isPrim = stack->IsPhysicalPrimary(label);
1131 isFromStrangess = IsFromStrangeness(label,stack);
1132 isFromConversion = IsFromConversion(label,stack);
1133 isFromMaterial = IsFromMaterial(label,stack);
1134 }
c683985a 1135
2e16763c 1136 //
1137 // TPC track
1138 //
1139 Int_t labelTPC = TMath::Abs(track->GetTPCLabel());
1140 if (labelTPC >= mcStackSize) continue;
1141 particleTPC = stack->Particle(labelTPC);
1142 if (!particleTPC) continue;
1143 if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)
1144 {
1145 particleMotherTPC = GetMother(particleTPC,stack);
1146 mechTPC = particleTPC->GetUniqueID();
1147 isPrimTPC = stack->IsPhysicalPrimary(labelTPC);
1148 isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);
1149 isFromConversionTPC = IsFromConversion(labelTPC,stack);
1150 isFromMaterialTPC = IsFromMaterial(labelTPC,stack);
1151 }
c683985a 1152
2e16763c 1153 //
1154 // store first track reference
1155 // for TPC track
1156 //
1157 TParticle *part=0;
1158 TClonesArray *trefs=0;
1159 Int_t status = mcEvent->GetParticleAndTR(TMath::Abs(labelTPC), part, trefs);
c683985a 1160
2e16763c 1161 if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.)
c683985a 1162 {
2e16763c 1163 Int_t nTrackRef = trefs->GetEntries();
1164 //printf("nTrackRef %d \n",nTrackRef);
c683985a 1165
2e16763c 1166 Int_t countITS = 0;
1167 for (Int_t iref = 0; iref < nTrackRef; iref++)
c683985a 1168 {
2e16763c 1169 AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);
1170
1171 // Ref. in the middle ITS
1172 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kITS)
1173 {
1174 nrefITS++;
1175 if(!refITS && countITS==2) {
1176 refITS = ref;
1177 //printf("refITS %p \n",refITS);
1178 }
1179 countITS++;
234cae2f 1180 }
c683985a 1181
2e16763c 1182 // TPC
1183 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC)
1184 {
1185 nrefTPC++;
1186 refTPCOut=ref;
1187 if(!refTPCIn) {
1188 refTPCIn = ref;
1189 //printf("refTPCIn %p \n",refTPCIn);
1190 //break;
1191 }
234cae2f 1192 }
2e16763c 1193 // TRD
1194 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD)
1195 {
1196 nrefTRD++;
1197 if(!refTRD) {
1198 refTRD = ref;
1199 }
234cae2f 1200 }
2e16763c 1201 // TOF
1202 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF)
1203 {
1204 nrefTOF++;
1205 if(!refTOF) {
1206 refTOF = ref;
1207 }
234cae2f 1208 }
2e16763c 1209 // EMCAL
1210 if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL)
1211 {
1212 nrefEMCAL++;
1213 if(!refEMCAL) {
1214 refEMCAL = ref;
1215 }
1216 }
1217 // PHOS
1218 // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS)
1219 // {
1220 // nrefPHOS++;
1221 // if(!refPHOS) {
1222 // refPHOS = ref;
1223 // }
1224 // }
c683985a 1225 }
2e16763c 1226
1227 // transform inner params to TrackRef
1228 // reference frame
1229 if(refTPCIn && trackInnerC3)
c683985a 1230 {
2e16763c 1231 Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());
1232 isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);
1233 isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);
c683985a 1234 }
234cae2f 1235 }
c683985a 1236
2e16763c 1237 //
1238 // ITS track
1239 //
1240 Int_t labelITS = TMath::Abs(track->GetITSLabel());
1241 if (labelITS >= mcStackSize) continue;
1242 particleITS = stack->Particle(labelITS);
1243 if (!particleITS) continue;
1244 if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)
234cae2f 1245 {
2e16763c 1246 particleMotherITS = GetMother(particleITS,stack);
1247 mechITS = particleITS->GetUniqueID();
1248 isPrimITS = stack->IsPhysicalPrimary(labelITS);
1249 isFromStrangessITS = IsFromStrangeness(labelITS,stack);
1250 isFromConversionITS = IsFromConversion(labelITS,stack);
1251 isFromMaterialITS = IsFromMaterial(labelITS,stack);
234cae2f 1252 }
c683985a 1253 }
2e16763c 1254 while (0);
c683985a 1255 }
1256
234cae2f 1257 //
2e16763c 1258 Bool_t dumpToTree=kFALSE;
1259
1260 if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE;
1261 //if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1262 if(isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
1263 if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE;
1264 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1265 if (fReducePileUp){
1266 //
1267 // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks
1268 // Done only in case no MC info
1269 //
1270 Float_t dcaTPC[2];
1271 track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]);
1272 Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10;
1273 Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin));
1274 Bool_t keepPileUp=gRandom->Rndm()<0.05;
1275 if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){
1276 dumpToTree=kFALSE;
1277 }
234cae2f 1278 }
c683985a 1279
2e16763c 1280 //init dummy objects
b3982e10 1281 static AliESDVertex dummyvtxESD;
1282 //if (!dummyvtxESD)
1283 //{
1284 // dummyvtxESD=new AliESDVertex();
1285 // //dummyvtxESD->SetNContributors(1);
1286 // //UShort_t pindices[1]; pindices[0]=0;
1287 // //dummyvtxESD->SetIndices(1,pindices);
1288 // //dummyvtxESD->SetNContributors(0);
1289 //}
1290 static AliExternalTrackParam dummyexternaltrackparam;
1291 //if (!dummyexternaltrackparam) dummyexternaltrackparam=new AliExternalTrackParam();
1292 static AliTrackReference dummytrackreference;
1293 //if (!dummytrackreference) dummytrackreference=new AliTrackReference();
1294 static TParticle dummyparticle;
1295 //if (!dummyparticle) dummyparticle=new TParticle();
2e16763c 1296
1297 //assign the dummy objects if needed
b3982e10 1298 if (!track) {track=fDummyTrack;}
b4990577 1299 AliESDfriendTrack *friendTrackStore=friendTrack; // store friend track for later processing
1300 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
1301 friendTrackStore = (gRandom->Rndm()<1./fFriendDownscaling)? friendTrack:0;
1302 }
1303 if (fFriendDownscaling<=0){
1304 if (((*fTreeSRedirector)<<"highPt").GetTree()){
1305 TTree * tree = ((*fTreeSRedirector)<<"highPt").GetTree();
1306 if (tree){
1307 Double_t sizeAll=tree->GetZipBytes();
1308 TBranch * br= tree->GetBranch("friendTrack.fPoints");
1309 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
82522b46 1310 br= tree->GetBranch("friendTrack.fCalibContainer");
b4990577 1311 if (br) sizeFriend+=br->GetZipBytes();
1312 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) friendTrackStore=0;
1313 }
1314 }
1315 }
1316
1317
1318 // if (!friendTrackStore && fFriendDownscaling<=1) {friendTrack=fDummyFriendTrack;}
b3982e10 1319 if (!vtxESD) {vtxESD=&dummyvtxESD;}
c683985a 1320 if (mcEvent)
1321 {
b3982e10 1322 if (!refTPCIn) {refTPCIn=&dummytrackreference;}
1323 if (!refITS) {refITS=&dummytrackreference;}
1324 if (!particle) {particle=&dummyparticle;}
1325 if (!particleMother) {particleMother=&dummyparticle;}
1326 if (!particleTPC) {particleTPC=&dummyparticle;}
1327 if (!particleMotherTPC) {particleMotherTPC=&dummyparticle;}
1328 if (!particleITS) {particleITS=&dummyparticle;}
1329 if (!particleMotherITS) {particleMotherITS=&dummyparticle;}
c683985a 1330 }
2e16763c 1331
1332 Int_t ntracks = esdEvent->GetNumberOfTracks();
1333
1334 // Global event id calculation using orbitID, bunchCrossingID and periodID
1335 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1336 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1337 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1338 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1339
1340 // fill histograms
1341 FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0));
1342
1a30582a 1343 if(fTreeSRedirector && dumpToTree && fFillTree) {
2e16763c 1344 (*fTreeSRedirector)<<"highPt"<<
1345 "gid="<<gid<<
1a30582a 1346 "fileName.="<<&fCurrentFileName<< // name of the chunk file (hopefully full)
2e16763c 1347 "runNumber="<<runNumber<< // runNumber
1348 "evtTimeStamp="<<evtTimeStamp<< // time stamp of event (in seconds)
1349 "evtNumberInFile="<<evtNumberInFile<< // event number
1350 "triggerClass="<<&triggerClass<< // trigger class as a string
1351 "Bz="<<bz<< // solenoid magnetic field in the z direction (in kGaus)
1352 "vtxESD.="<<vtxESD<< // vertexer ESD tracks (can be biased by TPC pileup tracks)
2e16763c 1353 "IRtot="<<ir1<< // interaction record (trigger) counters - coutner 1
1354 "IRint2="<<ir2<< // interaction record (trigger) coutners - counter 2
1355 "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
1356 "ntracks="<<ntracks<< // number of the esd tracks (to take into account the pileup in the TPC)
1357 // important variables for the pile-up studies
1358 "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
1359 "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
1360 "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
1361 "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
1362 "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
1363 "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
1364 //
1365 "esdTrack.="<<track<< // esdTrack as used in the physical analysis
b4990577 1366 // "friendTrack.="<<friendTrack<< // esdFriendTrack associated to the esdTrack
1367 "friendTrack.="<<friendTrackStore<< // esdFriendTrack associated to the esdTrack
2e16763c 1368 "extTPCInnerC.="<<tpcInnerC<< // ???
1369 "extInnerParamC.="<<trackInnerC<< // ???
1370 "extInnerParam.="<<trackInnerC2<< // ???
1371 "extOuterITS.="<<outerITSc<< // ???
1372 "extInnerParamRef.="<<trackInnerC3<< // ???
1373 "chi2TPCInnerC="<<chi2(0,0)<< // chi2 of tracks ???
1374 "chi2InnerC="<<chi2trackC(0,0)<< // chi2s of tracks TPCinner to the combined
1375 "chi2OuterITS="<<chi2OuterITS(0,0)<< // chi2s of tracks TPC at inner wall to the ITSout
1376 "centralityF="<<centralityF;
1a30582a 1377 if (mcEvent){
b3982e10 1378 static AliTrackReference refDummy;
2e16763c 1379 if (!refITS) refITS = &refDummy;
1380 if (!refTRD) refTRD = &refDummy;
1381 if (!refTOF) refTOF = &refDummy;
1382 if (!refEMCAL) refEMCAL = &refDummy;
1383 if (!refPHOS) refPHOS = &refDummy;
1384 (*fTreeSRedirector)<<"highPt"<<
1385 "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
1386 "nrefITS="<<nrefITS<< // number of track references in the ITS
1387 "nrefTPC="<<nrefTPC<< // number of track references in the TPC
1388 "nrefTRD="<<nrefTRD<< // number of track references in the TRD
1389 "nrefTOF="<<nrefTOF<< // number of track references in the TOF
1390 "nrefEMCAL="<<nrefEMCAL<< // number of track references in the TOF
1391 "nrefPHOS="<<nrefPHOS<< // number of track references in the TOF
1392 "refTPCIn.="<<refTPCIn<<
1393 "refTPCOut.="<<refTPCOut<<
1394 "refITS.="<<refITS<<
1395 "refTRD.="<<refTRD<<
1396 "refTOF.="<<refTOF<<
1397 "refEMCAL.="<<refEMCAL<<
1398 "refPHOS.="<<refPHOS<<
1399 "particle.="<<particle<<
1400 "particleMother.="<<particleMother<<
1401 "mech="<<mech<<
1402 "isPrim="<<isPrim<<
1403 "isFromStrangess="<<isFromStrangess<<
1404 "isFromConversion="<<isFromConversion<<
1405 "isFromMaterial="<<isFromMaterial<<
1406 "particleTPC.="<<particleTPC<<
1407 "particleMotherTPC.="<<particleMotherTPC<<
1408 "mechTPC="<<mechTPC<<
1409 "isPrimTPC="<<isPrimTPC<<
1410 "isFromStrangessTPC="<<isFromStrangessTPC<<
1411 "isFromConversionTPC="<<isFromConversionTPC<<
1412 "isFromMaterialTPC="<<isFromMaterialTPC<<
1413 "particleITS.="<<particleITS<<
1414 "particleMotherITS.="<<particleMotherITS<<
1415 "mechITS="<<mechITS<<
1416 "isPrimITS="<<isPrimITS<<
1417 "isFromStrangessITS="<<isFromStrangessITS<<
1418 "isFromConversionITS="<<isFromConversionITS<<
1419 "isFromMaterialITS="<<isFromMaterialITS;
1420 }
1421 //finish writing the entry
1422 AliInfo("writing tree highPt");
1423 (*fTreeSRedirector)<<"highPt"<<"\n";
1424 }
82522b46 1425 AliSysInfo::AddStamp("filteringTask",iTrack,numberOfTracks,numberOfFriendTracks,(friendTrackStore)?0:1);
2e16763c 1426 delete tpcInnerC;
1427 delete trackInnerC;
1428 delete trackInnerC2;
1429 delete outerITSc;
1430 delete trackInnerC3;
1431 }
c683985a 1432 }
c683985a 1433}
1434
1435
1436//_____________________________________________________________________________
1437void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
1438{
1439 //
1440 // Fill tree for efficiency studies MC only
1441 AliInfo("we start!");
c683985a 1442 if(!mcEvent) {
1443 AliDebug(AliLog::kError, "mcEvent not available");
1444 return;
1445 }
1446
1447 // get selection cuts
1448 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1449 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1450 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1451
1452 if(!evtCuts || !accCuts || !esdTrackCuts) {
1453 AliDebug(AliLog::kError, "cuts not available");
1454 return;
1455 }
1456
1457 // trigger selection
1458 Bool_t isEventTriggered = kTRUE;
1459 AliPhysicsSelection *physicsSelection = NULL;
1460 AliTriggerAnalysis* triggerAnalysis = NULL;
1461
1462 //
1463 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
c683985a 1464
c683985a 1465
1466 // trigger
1467 if(evtCuts->IsTriggerRequired())
1468 {
1469 // always MB
1470 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
c683985a 1471 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1472 if(!physicsSelection) return;
1473
1474 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1475 // set trigger (V0AND)
1476 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1477 if(!triggerAnalysis) return;
1478 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1479 }
1480 }
1481
1482 // centrality determination
1483 Float_t centralityF = -1;
1484 AliCentrality *esdCentrality = esdEvent->GetCentrality();
1485 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1486
1487 // use MC information
1488 AliHeader* header = 0;
1489 AliGenEventHeader* genHeader = 0;
1490 AliStack* stack = 0;
1491 Int_t mcStackSize=0;
1492 TArrayF vtxMC(3);
1493
1494 Int_t multMCTrueTracks = 0;
1495 //
1496 if(!mcEvent) {
1497 AliDebug(AliLog::kError, "mcEvent not available");
1498 return;
1499 }
1500 // get MC event header
1501 header = mcEvent->Header();
1502 if (!header) {
1503 AliDebug(AliLog::kError, "Header not available");
1504 return;
1505 }
1506 // MC particle stack
1507 stack = mcEvent->Stack();
1508 if (!stack) {
1509 AliDebug(AliLog::kError, "Stack not available");
1510 return;
1511 }
1512 mcStackSize=stack->GetNtrack();
1513
1514 // get MC vertex
1515 genHeader = header->GenEventHeader();
1516 if (!genHeader) {
1517 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1518 return;
1519 }
1520 genHeader->PrimaryVertex(vtxMC);
1521
1522 // multipliticy of all MC primary tracks
1523 // in Zv, pt and eta ranges)
1524 multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
1525
1526
1527 // get reconstructed vertex
1528 //const AliESDVertex* vtxESD = 0;
1529 AliESDVertex* vtxESD = 0;
1530 if(GetAnalysisMode() == kTPCAnalysisMode) {
1531 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
1532 }
1533 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
1534 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
1535 }
1536 else {
1537 return;
1538 }
1539
1540 if(!vtxESD) return;
1541 //
1542 // Vertex info comparison and track multiplicity
1543 //
1544 AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
1545 AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
1546 Int_t contSPD = vertexSPD->GetNContributors();
1547 Int_t contTPC = vertexTPC->GetNContributors();
1548 TVectorD vertexPosTPC(3), vertexPosSPD(3);
1549 vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
1550 vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
1551 Int_t ntracksTPC=0;
1552 Int_t ntracksITS=0;
1553 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
1554 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1555 if(!track) continue;
1556 if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
1557 if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
1558 }
1559
1560 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1561 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1562 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1563
1564 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1565
1566 // check event cuts
1567 if(isEventOK && isEventTriggered)
1568 {
1569 //if(!stack) return;
1570
1571 //
1572 // MC info
1573 //
1574 TParticle *particle=NULL;
1575 TParticle *particleMother=NULL;
1576 Int_t mech=-1;
1577
1578 // reco event info
1579 Double_t vert[3] = {0};
1580 vert[0] = vtxESD->GetXv();
1581 vert[1] = vtxESD->GetYv();
1582 vert[2] = vtxESD->GetZv();
1583 Int_t mult = vtxESD->GetNContributors();
1584 Double_t bz = esdEvent->GetMagneticField();
1585 Double_t runNumber = esdEvent->GetRunNumber();
1586 Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1587 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
234cae2f 1588 // loop over MC stack
c683985a 1589 for (Int_t iMc = 0; iMc < mcStackSize; ++iMc)
1590 {
1591 particle = stack->Particle(iMc);
1592 if (!particle)
1593 continue;
1594
1595 // only charged particles
1596 if(!particle->GetPDG()) continue;
1597 Double_t charge = particle->GetPDG()->Charge()/3.;
1598 if (TMath::Abs(charge) < 0.001)
1599 continue;
1600
1601 // only primary particles
1602 Bool_t prim = stack->IsPhysicalPrimary(iMc);
1603 if(!prim) continue;
1604
1605 // downscale low-pT particles
1606 Double_t scalempt= TMath::Min(particle->Pt(),10.);
1607 Double_t downscaleF = gRandom->Rndm();
1608 downscaleF *= fLowPtTrackDownscaligF;
1609 if(TMath::Exp(2*scalempt)<downscaleF) continue;
1610
1611 // is particle in acceptance
1612 if(!accCuts->AcceptTrack(particle)) continue;
1613
1614 // check if particle reconstructed
1615 Bool_t isRec = kFALSE;
1616 Int_t trackIndex = -1;
1617 Int_t trackLoopIndex = -1;
1618 Int_t isESDtrackCut= 0;
1619 Int_t isAccCuts = 0;
1620 Int_t nRec = 0; // how many times reconstructed
1621 Int_t nFakes = 0; // how many times reconstructed as a fake track
1622 AliESDtrack *recTrack = NULL;
1623
1624 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1625 {
1626 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1627 if(!track) continue;
1628 if(track->Charge()==0) continue;
234cae2f 1629 //
1630 Int_t label = TMath::Abs(track->GetLabel());
1631 if (label >= mcStackSize) continue;
1632 if(label == iMc) {
1633 Bool_t isAcc=esdTrackCuts->AcceptTrack(track);
1634 if (isAcc) isESDtrackCut=1;
1635 if (accCuts->AcceptTrack(track)) isAccCuts=1;
1636 isRec = kTRUE;
c683985a 1637 trackIndex = iTrack;
1638
234cae2f 1639 if (recTrack){
1640 if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track
1641 if (!isAcc) continue;
1642 trackLoopIndex = iTrack;
1643 }
1644 recTrack = esdEvent->GetTrack(trackIndex);
c683985a 1645 nRec++;
1646 if(track->GetLabel()<0) nFakes++;
1647
234cae2f 1648 continue;
1649 }
c683985a 1650 }
234cae2f 1651
c683985a 1652 // Store information in the output tree
1653 if (trackLoopIndex>-1) {
1654 recTrack = esdEvent->GetTrack(trackLoopIndex);
1655 } else if (trackIndex >-1) {
1656 recTrack = esdEvent->GetTrack(trackIndex);
1657 } else {
2e16763c 1658 recTrack = fDummyTrack;
c683985a 1659 }
1660
1661 particleMother = GetMother(particle,stack);
1662 mech = particle->GetUniqueID();
1663
1664 //MC particle track length
1665 Double_t tpcTrackLength = 0.;
1666 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
1667 if(mcParticle) {
1668 Int_t counter;
1669 tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
1670 }
1671
1672
1673 //
1674 if(fTreeSRedirector && fFillTree) {
1675 (*fTreeSRedirector)<<"MCEffTree"<<
1a30582a 1676 "fileName.="<<&fCurrentFileName<<
c683985a 1677 "triggerClass.="<<&triggerClass<<
1678 "runNumber="<<runNumber<<
1679 "evtTimeStamp="<<evtTimeStamp<<
1680 "evtNumberInFile="<<evtNumberInFile<< //
1681 "Bz="<<bz<< // magnetic field
1682 "vtxESD.="<<vtxESD<< // vertex info
234cae2f 1683 //
c683985a 1684 "mult="<<mult<< // primary vertex 9whatewe found) multiplicity
234cae2f 1685 "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
1686 // important variables for the pile-up studies
1687 "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
1688 "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
1689 "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
1690 "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
1691 "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
1692 "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
1693 //
1694 //
1695 "isAcc0="<<isESDtrackCut<< // track accepted by ESD track cuts
1696 "isAcc1="<<isAccCuts<< // track accepted by acceptance cuts flag
c683985a 1697 "esdTrack.="<<recTrack<< // reconstructed track (only the longest from the loopers)
1698 "isRec="<<isRec<< // track was reconstructed
1699 "tpcTrackLength="<<tpcTrackLength<< // track length in the TPC r projection
1700 "particle.="<<particle<< // particle properties
1701 "particleMother.="<<particleMother<< // particle mother
1702 "mech="<<mech<< // production mechanizm
1703 "nRec="<<nRec<< // how many times reconstruted
1704 "nFakes="<<nFakes<< // how many times reconstructed as a fake track
1705 "\n";
1706 }
1707
2e16763c 1708 //if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
c683985a 1709 }
1710 }
1711
1712}
1713
1714//_____________________________________________________________________________
1715Bool_t AliAnalysisTaskFilteredTree::IsHighDeDxParticle(AliESDtrack * track) {
1716 //
1717 // check if particle is Z > 1
1718 //
1719 if (track->GetTPCNcls() < 60) return kFALSE;
1720 Double_t mom = track->GetInnerParam()->GetP();
1721 if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
1722 Float_t dca[2], bCov[3];
1723 track->GetImpactParameters(dca,bCov);
1724 //
1725
1726 Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
1727
1728 if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
1729
1730 return kFALSE;
1731}
1732
1733//_____________________________________________________________________________
234cae2f 1734void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
c683985a 1735{
1736 //
1737 // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates
1738 //
1739 if(!esdEvent) {
1740 AliDebug(AliLog::kError, "esdEvent not available");
1741 return;
1742 }
1743
1744 // get selection cuts
1745 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1746 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1747 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1748
1749 if(!evtCuts || !accCuts || !esdTrackCuts) {
1750 AliDebug(AliLog::kError, "cuts not available");
1751 return;
1752 }
1753
1754 // trigger selection
1755 Bool_t isEventTriggered = kTRUE;
1756 AliPhysicsSelection *physicsSelection = NULL;
1757 AliTriggerAnalysis* triggerAnalysis = NULL;
1758
1759 //
1760 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
c683985a 1761
1762 // trigger
1763 if(evtCuts->IsTriggerRequired())
1764 {
1765 // always MB
1766 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1767
1768 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1769 if(!physicsSelection) return;
1770 //SetPhysicsTriggerSelection(physicsSelection);
1771
1772 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1773 // set trigger (V0AND)
1774 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1775 if(!triggerAnalysis) return;
1776 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1777 }
1778 }
1779
1780 // centrality determination
1781 Float_t centralityF = -1;
1782 AliCentrality *esdCentrality = esdEvent->GetCentrality();
1783 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
1784
1785
1786 // get reconstructed vertex
1787 //const AliESDVertex* vtxESD = 0;
1788 AliESDVertex* vtxESD = 0;
1789 if(GetAnalysisMode() == kTPCAnalysisMode) {
234cae2f 1790 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
c683985a 1791 }
1792 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
234cae2f 1793 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
c683985a 1794 }
1795 else {
234cae2f 1796 return;
c683985a 1797 }
1798
1799 if(!vtxESD) return;
1800
1801 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1802 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1803 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1804
1805 // check event cuts
1806 if(isEventOK && isEventTriggered) {
c683985a 1807 //
234cae2f 1808 // Dump the pt downscaled V0 into the tree
1809 //
1810 Int_t ntracks = esdEvent->GetNumberOfTracks();
1811 Int_t nV0s = esdEvent->GetNumberOfV0s();
1812 Int_t run = esdEvent->GetRunNumber();
1813 Int_t time= esdEvent->GetTimeStamp();
1814 Int_t evNr=esdEvent->GetEventNumberInFile();
1815 Double_t bz = esdEvent->GetMagneticField();
c683985a 1816
b55a470b 1817 // Global event id calculation using orbitID, bunchCrossingID and periodID
1818 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1819 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1820 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1821 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1822
234cae2f 1823
1824 for (Int_t iv0=0; iv0<nV0s; iv0++){
1825
1826 AliESDv0 * v0 = esdEvent->GetV0(iv0);
1827 if (!v0) continue;
1828 AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
1829 AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
1830 if (!track0) continue;
1831 if (!track1) continue;
1832 AliESDfriendTrack* friendTrack0=NULL;
1833 AliESDfriendTrack* friendTrack1=NULL;
82522b46 1834 if (esdFriend) {
1835 if (!esdFriend->TestSkipBit()){
234cae2f 1836 friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
1837 friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
1838 }
1839 }
1a30582a 1840 if (track0->GetSign()<0) {
1841 track1 = esdEvent->GetTrack(v0->GetIndex(0));
1842 track0 = esdEvent->GetTrack(v0->GetIndex(1));
1843 }
1844
82522b46 1845 //
1846 AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
1847 AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
1848 if (fFriendDownscaling>=1){ // downscaling number of friend tracks
1849 if (gRandom->Rndm()>1./fFriendDownscaling){
1850 friendTrackStore0 = 0;
1851 friendTrackStore1 = 0;
1852 }
1853 }
1854 if (fFriendDownscaling<=0){
1855 if (((*fTreeSRedirector)<<"V0s").GetTree()){
1856 TTree * tree = ((*fTreeSRedirector)<<"V0s").GetTree();
1857 if (tree){
1858 Double_t sizeAll=tree->GetZipBytes();
1859 TBranch * br= tree->GetBranch("friendTrack0.fPoints");
1860 Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
1861 br= tree->GetBranch("friendTrack0.fCalibContainer");
1862 if (br) sizeFriend+=br->GetZipBytes();
1863 if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
1864 friendTrackStore0=0;
1865 friendTrackStore1=0;
1866 }
1867 }
1868 }
1869 }
1870
234cae2f 1871 //
1872 Bool_t isDownscaled = IsV0Downscaled(v0);
1873 if (isDownscaled) continue;
1874 AliKFParticle kfparticle; //
1875 Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
1876 if (type==0) continue;
1877 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
1878
1879 if(!fFillTree) return;
1880 if(!fTreeSRedirector) return;
1881 (*fTreeSRedirector)<<"V0s"<<
1a30582a 1882 "gid="<<gid<< // global id of event
1883 "isDownscaled="<<isDownscaled<< //
234cae2f 1884 "triggerClass="<<&triggerClass<< // trigger
1a30582a 1885 "Bz="<<bz<< //
1886 "fileName.="<<&fCurrentFileName<< // full path - file name with ESD
1887 "runNumber="<<run<< //
1888 "evtTimeStamp="<<time<< // time stamp of event in secons
1889 "evtNumberInFile="<<evNr<< //
1890 "type="<<type<< // type of V0-
234cae2f 1891 "ntracks="<<ntracks<<
1892 "v0.="<<v0<<
1893 "kf.="<<&kfparticle<<
1a30582a 1894 "track0.="<<track0<< // track
234cae2f 1895 "track1.="<<track1<<
82522b46 1896 "friendTrack0.="<<friendTrackStore0<<
1897 "friendTrack1.="<<friendTrackStore1<<
234cae2f 1898 "centralityF="<<centralityF<<
1899 "\n";
234cae2f 1900 }
c683985a 1901 }
1902}
1903
1904//_____________________________________________________________________________
234cae2f 1905void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
c683985a 1906{
1907 //
1908 // Select real events with large TPC dEdx signal
1909 //
1910 if(!esdEvent) {
1911 AliDebug(AliLog::kError, "esdEvent not available");
1912 return;
1913 }
1914
1915 // get selection cuts
1916 AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
1917 AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
1918 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1919
1920 if(!evtCuts || !accCuts || !esdTrackCuts) {
1921 AliDebug(AliLog::kError, "cuts not available");
1922 return;
1923 }
1924
c683985a 1925
1926 // trigger
1927 Bool_t isEventTriggered = kTRUE;
1928 AliPhysicsSelection *physicsSelection = NULL;
1929 AliTriggerAnalysis* triggerAnalysis = NULL;
1930
1931 //
1932 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
c683985a 1933
1934 if(evtCuts->IsTriggerRequired())
1935 {
1936 // always MB
1937 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1938
1939 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1940 if(!physicsSelection) return;
1941
1942 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1943 // set trigger (V0AND)
1944 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1945 if(!triggerAnalysis) return;
1946 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1947 }
1948 }
1949
1950 // get reconstructed vertex
1951 AliESDVertex* vtxESD = 0;
1952 if(GetAnalysisMode() == kTPCAnalysisMode) {
234cae2f 1953 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
c683985a 1954 }
1955 else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
234cae2f 1956 vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
c683985a 1957 }
1958 else {
234cae2f 1959 return;
c683985a 1960 }
1961 if(!vtxESD) return;
1962
1963 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
1964 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1965 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1966
1967
1968 // check event cuts
1969 if(isEventOK && isEventTriggered)
1970 {
1971 Double_t vert[3] = {0};
1972 vert[0] = vtxESD->GetXv();
1973 vert[1] = vtxESD->GetYv();
1974 vert[2] = vtxESD->GetZv();
1975 Int_t mult = vtxESD->GetNContributors();
1976 Double_t bz = esdEvent->GetMagneticField();
1977 Double_t runNumber = esdEvent->GetRunNumber();
1978 Double_t evtTimeStamp = esdEvent->GetTimeStamp();
1979 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
b55a470b 1980
1981 // Global event id calculation using orbitID, bunchCrossingID and periodID
1982 ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
1983 ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
1984 ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
1985 ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
1986
c683985a 1987 // large dEdx
1988 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
1989 {
1990 AliESDtrack *track = esdEvent->GetTrack(iTrack);
1991 if(!track) continue;
234cae2f 1992 AliESDfriendTrack* friendTrack=NULL;
1993 if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
c683985a 1994 if(track->Charge()==0) continue;
1995 if(!esdTrackCuts->AcceptTrack(track)) continue;
1996 if(!accCuts->AcceptTrack(track)) continue;
1997
1998 if(!IsHighDeDxParticle(track)) continue;
1999 TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
2000
2001 if(!fFillTree) return;
2002 if(!fTreeSRedirector) return;
1a30582a 2003 (*fTreeSRedirector)<<"dEdx"<< // high dEdx tree
2004 "gid="<<gid<< // global id
2005 "fileName.="<<&fCurrentFileName<< // file name
234cae2f 2006 "runNumber="<<runNumber<<
2007 "evtTimeStamp="<<evtTimeStamp<<
2008 "evtNumberInFile="<<evtNumberInFile<<
2009 "triggerClass="<<&triggerClass<< // trigger
2010 "Bz="<<bz<<
1a30582a 2011 "vtxESD.="<<vtxESD<< //
234cae2f 2012 "mult="<<mult<<
2013 "esdTrack.="<<track<<
2014 "friendTrack.="<<friendTrack<<
2015 "\n";
c683985a 2016 }
2017 }
2018}
2019
2020//_____________________________________________________________________________
2021Int_t AliAnalysisTaskFilteredTree::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
2022{
2023 //
2024 // Create KF particle in case the V0 fullfill selection criteria
2025 //
2026 // Selection criteria
2027 // 0. algorithm cut
2028 // 1. track cut
2029 // 3. chi2 cut
2030 // 4. rough mass cut
2031 // 5. Normalized pointing angle cut
2032 //
2033 const Double_t cutMass=0.2;
2034 const Double_t kSigmaDCACut=3;
2035 //
2036 // 0.) algo cut - accept only on the fly
2037 //
2038 if (v0->GetOnFlyStatus() ==kFALSE) return 0;
2039 //
2040 // 1.) track cut
2041 //
2042 AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
2043 AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
2044 /*
234cae2f 2045 TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
2046 TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
2047 TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
2048 */
c683985a 2049 if (TMath::Abs(track0->GetTgl())>1) return 0;
2050 if (TMath::Abs(track1->GetTgl())>1) return 0;
2051 if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
2052 if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
c683985a 2053 Float_t pos0[2]={0}, cov0[3]={0};
2054 Float_t pos1[2]={0}, cov1[3]={0};
2055 track0->GetImpactParameters(pos0,cov0);
2056 track0->GetImpactParameters(pos1,cov1);
2057 //
2058 if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
2059 if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
2060 //
2061 //
2062 // 3.) Chi2 cut
2063 //
2064 Double_t chi2KF = v0->GetKFInfo(2,2,2);
2065 if (chi2KF>25) return 0;
2066 //
2067 // 4.) Rough mass cut - 0.200 GeV
2068 //
2069 static Double_t masses[2]={-1};
2070 if (masses[0]<0){
2071 masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
2072 masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
2073 }
2074 Double_t mass00= v0->GetEffMass(0,0);
2075 Double_t mass22= v0->GetEffMass(2,2);
2076 Double_t mass42= v0->GetEffMass(4,2);
2077 Double_t mass24= v0->GetEffMass(2,4);
2078 Bool_t massOK=kFALSE;
2079 Int_t type=0;
2080 Int_t ptype=0;
2081 Double_t dmass=1;
2082 Int_t p1=0, p2=0;
2083 if (TMath::Abs(mass00-0)<cutMass) {
2084 massOK=kTRUE; type+=1;
2085 if (TMath::Abs(mass00-0)<dmass) {
2086 ptype=1;
2087 dmass=TMath::Abs(mass00-0);
2088 p1=0; p2=0;
2089 }
2090 }
2091 if (TMath::Abs(mass24-masses[1])<cutMass) {
2092 massOK=kTRUE; type+=2;
2093 if (TMath::Abs(mass24-masses[1])<dmass){
2094 dmass = TMath::Abs(mass24-masses[1]);
2095 ptype=2;
2096 p1=2; p2=4;
2097 }
2098 }
2099 if (TMath::Abs(mass42-masses[1])<cutMass) {
2100 massOK=kTRUE; type+=4;
2101 if (TMath::Abs(mass42-masses[1])<dmass){
2102 dmass = TMath::Abs(mass42-masses[1]);
2103 ptype=4;
2104 p1=4; p2=2;
2105 }
2106 }
2107 if (TMath::Abs(mass22-masses[0])<cutMass) {
2108 massOK=kTRUE; type+=8;
2109 if (TMath::Abs(mass22-masses[0])<dmass){
2110 dmass = TMath::Abs(mass22-masses[0]);
2111 ptype=8;
2112 p1=2; p2=2;
2113 }
2114 }
2115 if (type==0) return 0;
2116 //
2117 const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
2118 const AliExternalTrackParam *paramP = v0->GetParamP();
2119 const AliExternalTrackParam *paramN = v0->GetParamN();
2120 if (paramP->GetSign()<0){
2121 paramP=v0->GetParamP();
2122 paramN=v0->GetParamN();
2123 }
2124 //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
2125 //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
2126 //
2127 AliKFParticle kfp1( *paramP, spdg[p1] );
2128 AliKFParticle kfp2( *paramN, -1 *spdg[p2] );
2129 AliKFParticle V0KF;
2130 (V0KF)+=kfp1;
2131 (V0KF)+=kfp2;
2132 kfparticle=V0KF;
2133 //
2134 // Pointing angle
2135 //
2136 Double_t errPhi = V0KF.GetErrPhi();
2137 Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
2138 if (pointAngle/errPhi>10) return 0;
2139 //
2140 return ptype;
2141}
2142
2143//_____________________________________________________________________________
2144Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0)
2145{
2146 //
2147 // Downscale randomly low pt V0
2148 //
2149 //return kFALSE;
2150 Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
2151 Double_t scalempt= TMath::Min(maxPt,10.);
2152 Double_t downscaleF = gRandom->Rndm();
2153 downscaleF *= fLowPtV0DownscaligF;
2154 //
2155 // Special treatment of the gamma conversion pt spectra is softer -
2156 Double_t mass00= v0->GetEffMass(0,0);
2157 const Double_t cutMass=0.2;
2158 if (TMath::Abs(mass00-0)<cutMass){
2159 downscaleF/=10.; // 10 times smaller downscaling for the gamma concersion candidate
2160 }
2161 //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
2162 if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;
2163 return kFALSE;
2164
2165 /*
234cae2f 2166 TH1F his1("his1","his1",100,0,10);
2167 TH1F his2("his2","his2",100,0,10);
2168 {for (Int_t i=0; i<10000; i++){
2169 Double_t rnd=gRandom->Exp(1);
2170 Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
2171 his1->Fill(rnd);
2172 if (!isDownscaled) his2->Fill(rnd);
2173 }}
c683985a 2174
234cae2f 2175*/
c683985a 2176
2177}
2178
2179
2180
2181//_____________________________________________________________________________
2182Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
2183{
234cae2f 2184 // Constrain TPC inner params constrained
2185 //
2186 if(!tpcInnerC) return kFALSE;
2187 if(!vtx) return kFALSE;
c683985a 2188
234cae2f 2189 Double_t dz[2],cov[3];
2190 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2191 //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2192 //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
2193 if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
c683985a 2194
2195
234cae2f 2196 Double_t covar[6]; vtx->GetCovMatrix(covar);
2197 Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
2198 Double_t c[3]={covar[2],0.,covar[5]};
2199 Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
2200 if (chi2C>kVeryBig) return kFALSE;
c683985a 2201
234cae2f 2202 if(!tpcInnerC->Update(p,c)) return kFALSE;
c683985a 2203
2204 return kTRUE;
2205}
2206
2207//_____________________________________________________________________________
2208Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
2209{
234cae2f 2210 // Constrain TPC inner params constrained
2211 //
2212 if(!trackInnerC) return kFALSE;
2213 if(!vtx) return kFALSE;
c683985a 2214
234cae2f 2215 const Double_t kRadius = 2.8;
2216 const Double_t kMaxStep = 1.0;
c683985a 2217
234cae2f 2218 Double_t dz[2],cov[3];
c683985a 2219
234cae2f 2220 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
2221 //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
2222 //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
c683985a 2223
234cae2f 2224 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
2225 if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
c683985a 2226
234cae2f 2227 //
2228 Double_t covar[6]; vtx->GetCovMatrix(covar);
2229 Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
2230 Double_t c[3]={covar[2],0.,covar[5]};
2231 Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
2232 if (chi2C>kVeryBig) return kFALSE;
c683985a 2233
234cae2f 2234 if(!trackInnerC->Update(p,c)) return kFALSE;
c683985a 2235
2236 return kTRUE;
2237}
2238
2239
2240//_____________________________________________________________________________
2241TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack)
2242{
2243 if(!particle) return NULL;
2244 if(!stack) return NULL;
2245
2246 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2247 TParticle* mother = NULL;
2248 Int_t mcStackSize=stack->GetNtrack();
2249 if (motherLabel>=mcStackSize) return NULL;
2250 mother = stack->Particle(motherLabel);
2251
234cae2f 2252 return mother;
c683985a 2253}
2254
2255//_____________________________________________________________________________
34e8fe1d 2256Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(Int_t label, AliStack *const stack)
c683985a 2257{
2258 Bool_t isFromConversion = kFALSE;
2259
2260 if(stack) {
2261 Int_t mcStackSize=stack->GetNtrack();
2262 if (label>=mcStackSize) return kFALSE;
2263 TParticle* particle = stack->Particle(label);
2264 if (!particle) return kFALSE;
2265
2266 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2267 {
234cae2f 2268 Int_t mech = particle->GetUniqueID(); // production mechanism
2269 Bool_t isPrim = stack->IsPhysicalPrimary(label);
c683985a 2270
234cae2f 2271 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2272 if (motherLabel>=mcStackSize) return kFALSE;
2273 TParticle* mother = stack->Particle(motherLabel);
2274 if(mother) {
2275 Int_t motherPdg = mother->GetPdgCode();
c683985a 2276
234cae2f 2277 if(!isPrim && mech==5 && motherPdg==kGamma) {
2278 isFromConversion=kTRUE;
2279 }
2280 }
c683985a 2281 }
2282 }
2283
2284 return isFromConversion;
2285}
2286
2287//_____________________________________________________________________________
34e8fe1d 2288Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(Int_t label, AliStack *const stack)
c683985a 2289{
2290 Bool_t isFromMaterial = kFALSE;
2291
2292 if(stack) {
2293 Int_t mcStackSize=stack->GetNtrack();
2294 if (label>=mcStackSize) return kFALSE;
2295 TParticle* particle = stack->Particle(label);
2296 if (!particle) return kFALSE;
2297
2298 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2299 {
234cae2f 2300 Int_t mech = particle->GetUniqueID(); // production mechanism
2301 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2302
2303 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2304 if (motherLabel>=mcStackSize) return kFALSE;
2305 TParticle* mother = stack->Particle(motherLabel);
2306 if(mother) {
2307 if(!isPrim && mech==13) {
2308 isFromMaterial=kTRUE;
2309 }
2310 }
2311 }
c683985a 2312 }
2313
2314 return isFromMaterial;
2315}
2316
2317//_____________________________________________________________________________
34e8fe1d 2318Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(Int_t label, AliStack *const stack)
c683985a 2319{
2320 Bool_t isFromStrangeness = kFALSE;
2321
2322 if(stack) {
2323 Int_t mcStackSize=stack->GetNtrack();
2324 if (label>=mcStackSize) return kFALSE;
2325 TParticle* particle = stack->Particle(label);
2326 if (!particle) return kFALSE;
2327
2328 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
2329 {
234cae2f 2330 Int_t mech = particle->GetUniqueID(); // production mechanism
2331 Bool_t isPrim = stack->IsPhysicalPrimary(label);
2332
2333 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
2334 if (motherLabel>=mcStackSize) return kFALSE;
2335 TParticle* mother = stack->Particle(motherLabel);
2336 if(mother) {
2337 Int_t motherPdg = mother->GetPdgCode();
2338
2339 // K+-, lambda, antilambda, K0s decays
2340 if(!isPrim && mech==4 &&
2341 (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
2342 {
2343 isFromStrangeness = kTRUE;
2344 }
2345 }
c683985a 2346 }
2347 }
2348
2349 return isFromStrangeness;
2350}
2351
2352
2353//_____________________________________________________________________________
2354void AliAnalysisTaskFilteredTree::FinishTaskOutput()
2355{
2356 //
2357 // Called one at the end
2358 // locally on working node
2359 //
c683985a 2360 Bool_t deleteTrees=kTRUE;
2361 if ((AliAnalysisManager::GetAnalysisManager()))
2362 {
2363 if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
234cae2f 2364 AliAnalysisManager::kProofAnalysis)
c683985a 2365 deleteTrees=kFALSE;
2366 }
2367 if (deleteTrees) delete fTreeSRedirector;
2368 fTreeSRedirector=NULL;
c683985a 2369}
2370
2371//_____________________________________________________________________________
2372void AliAnalysisTaskFilteredTree::Terminate(Option_t *)
2373{
1a30582a 2374 //
c683985a 2375 // Called one at the end
1a30582a 2376 //
c683985a 2377}
2378
2379//_____________________________________________________________________________
2380Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
2381{
2382 //
2383 // calculate mc event true track multiplicity
2384 //
2385 if(!mcEvent) return 0;
2386
2387 AliStack* stack = 0;
2388 Int_t mult = 0;
2389
2390 // MC particle stack
2391 stack = mcEvent->Stack();
2392 if (!stack) return 0;
2393
2394 //
2395 //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
2396 //
2397
2398 Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
2399 if(!isEventOK) return 0;
2400
2401 Int_t nPart = stack->GetNtrack();
2402 for (Int_t iMc = 0; iMc < nPart; ++iMc)
2403 {
234cae2f 2404 TParticle* particle = stack->Particle(iMc);
2405 if (!particle)
2406 continue;
2407
2408 // only charged particles
2409 if(!particle->GetPDG()) continue;
2410 Double_t charge = particle->GetPDG()->Charge()/3.;
2411 if (TMath::Abs(charge) < 0.001)
2412 continue;
2413
2414 // physical primary
2415 Bool_t prim = stack->IsPhysicalPrimary(iMc);
2416 if(!prim) continue;
2417
2418 // checked accepted without pt cut
2419 //if(accCuts->AcceptTrack(particle))
2420 if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
2421 {
2422 mult++;
2423 }
c683985a 2424 }
2425
234cae2f 2426 return mult;
c683985a 2427}
2428
2429//_____________________________________________________________________________
34e8fe1d 2430void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC)
c683985a 2431{
234cae2f 2432 //
2433 // Fill pT relative resolution histograms for
2434 // TPC only, TPC only constrained to vertex and TPC+ITS tracking
2435 //
2436 if(!ptrack) return;
2437 if(!ptpcInnerC) return;
2438
2439 const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();
2440 if(!innerParam) return;
2441
2442 Float_t dxy, dz;
2443 ptrack->GetImpactParameters(dxy,dz);
2444
2445 // TPC+ITS primary tracks
2446 if( abs(ptrack->Eta())<0.8 &&
2447 ptrack->GetTPCClusterInfo(3,1)>120 &&
2448 ptrack->IsOn(0x40) &&
2449 ptrack->GetTPCclusters(0)>0.0 &&
2450 ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2451 //abs(innerParam->GetX())>0.0 &&
2452 //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2453 //abs(innerParam->GetTgl())<0.85 &&
2454 ptrack->IsOn(0x0004) &&
2455 ptrack->GetNcls(0)>0 &&
2456 ptrack->GetITSchi2()>0 &&
2457 sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 &&
2458 sqrt(chi2TPCInnerC)<6 &&
2459 (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) &&
2460 abs(dz)<2.0 &&
2461 abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )
2462 {
2463 fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2464 fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2465 fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
2466 }
c683985a 2467
234cae2f 2468 // TPC primary tracks
2469 // and TPC constrained primary tracks
c683985a 2470
234cae2f 2471 AliExternalTrackParam *ptpcInner = (AliExternalTrackParam *) ptrack->GetTPCInnerParam();
2472 if(!ptpcInner) return;
c683985a 2473
2474
234cae2f 2475 Float_t dxyTPC, dzTPC;
2476 ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);
c683985a 2477
234cae2f 2478 if( abs(ptrack->Eta())<0.8 &&
2479 ptrack->GetTPCClusterInfo(3,1)>120 &&
2480 ptrack->IsOn(0x40)&&
2481 ptrack->GetTPCclusters(0)>0.0 &&
2482 ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
2483 //abs(innerParam->GetX())>0.0 &&
2484 //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
2485 //abs(innerParam->GetTgl())<0.85 &&
2486 abs(dzTPC)<3.2 &&
2487 abs(dxyTPC)<2.4 )
2488 {
2489 // TPC only
2490 fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2491 fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2492 fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
2493
2494 // TPC constrained to vertex
2495 fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2496 fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2497 fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
2498 }
c683985a 2499}
4a00e62d 2500
2501
2502void AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend){
2503 //
2504 // Process ITS standalone tracks find match with closest TPC(or combined tracks) tracks
2505 // marian.ivanov@cern.ch
2506 // 0.) Init variables
2507 // 1.) GetTrack parameters at TPC inner wall
2508 // 2.) Match closest TPC track (STANDALONE/global) - chi2 match criteria
2509 //
2510 // Logic to be used in reco:
2511 // 1.) Find matching ITSalone->TPCalone
2512 // 2.) if (!TPCalone.FindClose(TPCother)) TPCalone.Addopt(ITSalone)
2513 // 3.) ff ((ITSalone.FindClose(Global)==0) CreateGlobaltrack
2514 const Double_t radiusMatch=84.; // redius to propagate
2515 //
2516 const Double_t dFastPhiCut=0.2; // 6 sigma (200 MeV) fast angular cut
2517 const Double_t dFastThetaCut=0.12; // 6 sigma (200 MeV) fast angular cut
2518 const Double_t dFastPosPhiCut=0.06; // 6 sigma (200 MeV) fast angular cut
2519 const Double_t dFastZCut=6; // 6 sigma (200 MeV) fast z difference cut
2520 const Double_t dFastPtCut=2.; // 6 sigma (200 MeV) fast 1/pt cut
2521 const Double_t chi2Cut=100; // chi2 matching cut
2522 //
2523 if (!esdFriend) return; // not ITS standalone track
2524 if (esdFriend->TestSkipBit()) return; // friends tracks not stored
2525 Int_t ntracks=esdEvent->GetNumberOfTracks();
2526 Float_t bz = esdEvent->GetMagneticField();
2527 //
2528 // 0.) Get parameters in reference radius TPC Inner wall
2529 //
2530 //
2531 TMatrixD vecPosR0(ntracks,6); // possition and momentum estimate at reference radius
2532 TMatrixD vecMomR0(ntracks,6); //
2533 TMatrixD vecPosR1(ntracks,6); // possition and momentum estimate at reference radius TPC track
2534 TMatrixD vecMomR1(ntracks,6); //
2535 Double_t xyz[3], pxyz[3]; //
2536 for (Int_t iTrack=0; iTrack<ntracks; iTrack++){
2537 AliESDtrack *track = esdEvent->GetTrack(iTrack);
2538 if(!track) continue;
2539 if (track->GetInnerParam()){
2540 const AliExternalTrackParam *trackTPC=track->GetInnerParam();
2541 trackTPC->GetXYZAt(radiusMatch,bz,xyz);
2542 trackTPC->GetPxPyPzAt(radiusMatch,bz,pxyz);
2543 for (Int_t i=0; i<3; i++){
2544 vecPosR1(iTrack,i)=xyz[i];
2545 vecMomR1(iTrack,i)=pxyz[i];
2546 }
2547 vecPosR1(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); // phi pos angle
2548 vecMomR1(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); // phi mom angle
2549 vecMomR1(iTrack,4)= trackTPC->GetSigned1Pt();;
2550 vecMomR1(iTrack,5)= trackTPC->GetTgl();;
2551 }
2552 AliESDfriendTrack* friendTrack=esdFriend->GetTrack(iTrack);
2553 if(!friendTrack) continue;
2554 if (friendTrack->GetITSOut()){
2555 const AliExternalTrackParam *trackITS=friendTrack->GetITSOut();
2556 trackITS->GetXYZAt(radiusMatch,bz,xyz);
2557 trackITS->GetPxPyPzAt(radiusMatch,bz,pxyz);
2558 for (Int_t i=0; i<3; i++){
2559 vecPosR0(iTrack,i)=xyz[i];
2560 vecMomR0(iTrack,i)=pxyz[i];
2561 }
2562 vecPosR0(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]);
2563 vecMomR0(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]);
2564 vecMomR0(iTrack,4)= trackITS->GetSigned1Pt();;
2565 vecMomR0(iTrack,5)= trackITS->GetTgl();;
2566 }
2567 }
2568 //
2569 // 1.) Find closest matching tracks, between the ITS standalone track
2570 // and the all other tracks
2571 // a.) caltegory - All
2572 // b.) category - without ITS
2573 //
2574 //
2575 Int_t ntracksPropagated=0;
2576 AliExternalTrackParam extTrackDummy;
2577 AliESDtrack esdTrackDummy;
2578 AliExternalTrackParam itsAtTPC;
2579 AliExternalTrackParam itsAtITSTPC;
2580 for (Int_t iTrack0=0; iTrack0<ntracks; iTrack0++){
2581 AliESDtrack *track0 = esdEvent->GetTrack(iTrack0);
2582 if(!track0) continue;
2583 if (track0->IsOn(AliVTrack::kTPCin)) continue;
2584 AliESDfriendTrack* friendTrack0=esdFriend->GetTrack(iTrack0);
2585 if (!friendTrack0) continue;
2586 //if (!track0->IsOn(AliVTrack::kITSpureSA)) continue;
2587 //if (!friendTrack0->GetITSOut()) continue; // is there flag for ITS standalone?
2588 ntracksPropagated++;
2589 //
2590 // 2.) find clostest TPCtrack
2591 // a.) all tracks
2592 Double_t minChi2All=10000000;
2593 Double_t minChi2TPC=10000000;
2594 Double_t minChi2TPCITS=10000000;
2595 Int_t indexAll=-1;
2596 Int_t indexTPC=-1;
2597 Int_t indexTPCITS=-1;
2598 Int_t ncandidates0=0; // n candidates - rough cut
2599 Int_t ncandidates1=0; // n candidates - rough + chi2 cut
2600 itsAtTPC=*(friendTrack0->GetITSOut());
2601 itsAtITSTPC=*(friendTrack0->GetITSOut());
2602 for (Int_t iTrack1=0; iTrack1<ntracks; iTrack1++){
2603 AliESDtrack *track1 = esdEvent->GetTrack(iTrack1);
2604 if(!track1) continue;
2605 if (!track1->IsOn(AliVTrack::kTPCin)) continue;
2606 // fast checks
2607 //
2608 if (TMath::Abs(vecPosR1(iTrack1,2)-vecPosR0(iTrack0,2))>dFastZCut) continue;
2609 if (TMath::Abs(vecPosR1(iTrack1,3)-vecPosR0(iTrack0,3))>dFastPosPhiCut) continue;
2610 if (TMath::Abs(vecMomR1(iTrack1,3)-vecMomR0(iTrack0,3))>dFastPhiCut) continue;
2611 if (TMath::Abs(vecMomR1(iTrack1,5)-vecMomR0(iTrack0,5))>dFastThetaCut) continue;
2612 if (TMath::Abs(vecMomR1(iTrack1,4)-vecMomR0(iTrack0,4))>dFastPtCut) continue;
2613 ncandidates0++;
2614 //
2615 const AliExternalTrackParam * param1= track1->GetInnerParam();
2616 if (!friendTrack0->GetITSOut()) continue;
2617 AliExternalTrackParam outerITS = *(friendTrack0->GetITSOut());
2618 if (!outerITS.Rotate(param1->GetAlpha())) continue;
2619 if (!outerITS.PropagateTo(param1->GetX(),bz)) continue; // assume track close to the TPC inner wall
2620 Double_t chi2 = outerITS.GetPredictedChi2(param1);
2621 if (chi2>chi2Cut) continue;
2622 ncandidates1++;
2623 if (chi2<minChi2All){
2624 minChi2All=chi2;
2625 indexAll=iTrack1;
2626 }
2627 if (chi2<minChi2TPC && track1->IsOn(AliVTrack::kITSin)==0){
2628 minChi2TPC=chi2;
2629 indexTPC=iTrack1;
2630 itsAtTPC=outerITS;
2631 }
2632 if (chi2<minChi2TPCITS && track1->IsOn(AliVTrack::kITSin)){
2633 minChi2TPCITS=chi2;
2634 indexTPCITS=iTrack1;
2635 itsAtITSTPC=outerITS;
2636 }
2637 }
2638 //
2639 AliESDtrack * trackAll= (indexAll>=0)? esdEvent->GetTrack(indexAll):&esdTrackDummy;
2640 AliESDtrack * trackTPC= (indexTPC>=0)? esdEvent->GetTrack(indexTPC):&esdTrackDummy;
2641 AliESDtrack * trackTPCITS= (indexTPCITS>=0)? esdEvent->GetTrack(indexTPCITS):&esdTrackDummy;
2642 (*fTreeSRedirector)<<"itsTPC"<<
2643 "indexAll="<<indexAll<< // index of closest track (chi2)
2644 "indexTPC="<<indexTPC<< // index of closest TPCalone tracks
2645 "indexTPCITS="<<indexTPCITS<< // index of closest cobined tracks
2646 "ncandidates0="<<ncandidates0<< // number of candidates
2647 "ncandidates1="<<ncandidates1<<
2648 //
2649 "chi2All="<<minChi2All<< // chi2 of closest tracks
2650 "chi2TPC="<<minChi2TPC<<
2651 "chi2TPCITS="<<minChi2TPCITS<<
2652 //
2653 "track0.="<<track0<< // ITS standalone tracks
2654 "trackAll.="<<trackAll<< // Closets other track
2655 "trackTPC.="<<trackTPC<< // Closest TPC only track
2656 "trackTPCITS.="<<trackTPCITS<< // closest combined track
2657 //
2658 "itsAtTPC.="<<&itsAtTPC<< // ITS track parameters at the TPC alone track frame
2659 "itsAtITSTPC.="<<&itsAtITSTPC<< // ITS track parameters at the TPC combeined track frame
2660 "\n";
2661 }
2662}
2663
b4990577 2664void AliAnalysisTaskFilteredTree::ProcessTrackMatch(AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/){
4a00e62d 2665/*
b4990577 2666 Use TPC standalone, ITS standalone and combined tracks to categorize/resp. recover track information.
2667
2668 Track categories:
2669 -TPC+ITS
2670 -TPC only
2671 -ITS only
2672 Topology categories:
2673 -Nice isolated tracks with full information
2674 -Overlapped tracks - Refit and sign them
2675 -Multiple found (check overlap factor) - Merge and sign
2676 -Charge particle (kink) decays - Change of direction - Sign them)
2677 Info:
2678 -Array of indexes of closest tracks in each individual category
2679 -Chi2 distance to the closest tracks at reference radius of TPCin
2680 -Overlap factors - fraction of shared clusters or missing region
2681 -Chi2 distance at DCA
2682 Information matrix:
2683 -matrix closest tracks from each categories
2684 -characterization - chi2, index,chi2, overlap ratio
2685
2686 Decissions:
2687 0.) Kink decay or catastophic multiple scaterring
2688 (combining all track categories)
2689 - small chi2 at the DCA
2690 - significantly large deflection angle
2691 - Combinatorial algorithm - to decrease CPU time restriction of investigation to tracks with small DCA at refernce radius
4a00e62d 2692
b4990577 2693 1.) if (TPConly && !(TPC+ITS) && ITSonly match ) {
2694 if (!kink) TPCOnly.addoptITS
2695 if (kink) TPConly sign
2696 }
4a00e62d 2697
b4990577 2698 2.) Overlap tracks - Refit with doUnfold
4a00e62d 2699
2700*/
b4990577 2701}