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