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