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