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