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