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