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