]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtCutAnalysis.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Enlarging window for DCS DPs retrieval for short runs for GRP + Keeping connection...
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / AlidNdPtCutAnalysis.cxx
CommitLineData
a65a7e70 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15//------------------------------------------------------------------------------
16// AlidNdPtCutAnalysis class.
17//
18// a. functionality:
19// - fills generic cut histograms
20// - generates cuts (selection criteria)
21//
22// b. data members:
23// - generic cut histograms
24// - control histograms
25//
26// Author: J.Otwinowski 04/11/2008
27//------------------------------------------------------------------------------
28#include "TH1.h"
29#include "TH2.h"
30
31#include "AliHeader.h"
32#include "AliGenEventHeader.h"
33#include "AliInputEventHandler.h"
34#include "AliStack.h"
35#include "AliESDEvent.h"
36#include "AliMCEvent.h"
37#include "AliESDtrackCuts.h"
38#include "AliLog.h"
39#include "AliTracker.h"
40
41#include "AlidNdPtEventCuts.h"
42#include "AlidNdPtAcceptanceCuts.h"
43#include "AlidNdPtBackgroundCuts.h"
44#include "AlidNdPtAnalysis.h"
45#include "AliPhysicsSelection.h"
46
47#include "AliPWG0Helper.h"
48#include "AlidNdPtHelper.h"
49#include "AlidNdPtCutAnalysis.h"
50
51using namespace std;
52
53ClassImp(AlidNdPtCutAnalysis)
54
55//_____________________________________________________________________________
56 AlidNdPtCutAnalysis::AlidNdPtCutAnalysis(): AlidNdPt(),
57 fAnalysisFolder(0),
58 fEventCount(0),
59 fRecEventHist(0),
60 fMCEventHist(0),
61 fRecMCEventHist(0),
62 fRecMCTrackHist(0)
63{
64 // default constructor
65 Init();
66}
67
68//_____________________________________________________________________________
69AlidNdPtCutAnalysis::AlidNdPtCutAnalysis(Char_t* name, Char_t* title): AlidNdPt(name,title),
70 fAnalysisFolder(0),
71 fEventCount(0),
72 fRecEventHist(0),
73 fMCEventHist(0),
74 fRecMCEventHist(0),
75 fRecMCTrackHist(0)
76{
77 // constructor
78 Init();
79}
80
81//_____________________________________________________________________________
82AlidNdPtCutAnalysis::~AlidNdPtCutAnalysis() {
83 //
84 if(fEventCount) delete fEventCount; fEventCount=0;
85 if(fRecEventHist) delete fRecEventHist; fRecEventHist=0;
86 if(fMCEventHist) delete fMCEventHist; fMCEventHist=0;
87 if(fRecMCEventHist) delete fRecMCEventHist; fRecMCEventHist=0;
88 if(fRecMCTrackHist) delete fRecMCTrackHist; fRecMCTrackHist=0;
89
90 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
91}
92
93//_____________________________________________________________________________
94void AlidNdPtCutAnalysis::Init(){
95 //
96 // Init histograms
97 //
98 /*
99 const Int_t ptNbins = 56;
100 const Double_t ptMin = 0.;
101 const Double_t ptMax = 16.;
102 Double_t binsPt[ptNbins+1] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0};
103 */
104 // set pt bins
105 const Int_t ptNbins = 50;
106 const Double_t ptMin = 1.e-2, ptMax = 50.;
107 Double_t *binsPt = CreateLogAxis(ptNbins,ptMin,ptMax);
108
109 //
110 Int_t binsEventCount[2]={2,2};
111 Double_t minEventCount[2]={0,0};
112 Double_t maxEventCount[2]={2,2};
113 fEventCount = new THnSparseF("fEventCount","trig vs trig+vertex",2,binsEventCount,minEventCount,maxEventCount);
114 fEventCount->GetAxis(0)->SetTitle("trig");
115 fEventCount->GetAxis(1)->SetTitle("trig+vert");
116 fEventCount->Sumw2();
117
118 //Xv:Yv:Zv:ResZv:Mult
119 Double_t kFact = 1.0;
120
121 Int_t binsRecEventHist[5]={80,80,100,80,150};
122 Double_t minRecEventHist[5]={-3.*kFact,-3.*kFact,-35.,0.,0.};
123 Double_t maxRecEventHist[5]={3.*kFact,3.*kFact,35.,10.,150.};
124 fRecEventHist = new THnSparseF("fRecEventHist","Xv:Yv:Zv:ResZv:Mult",5,binsRecEventHist,minRecEventHist,maxRecEventHist);
125 fRecEventHist->GetAxis(0)->SetTitle("Xv (cm)");
126 fRecEventHist->GetAxis(1)->SetTitle("Yv (cm)");
127 fRecEventHist->GetAxis(2)->SetTitle("Zv (cm)");
128 fRecEventHist->GetAxis(3)->SetTitle("ResZv (cm)");
129 fRecEventHist->GetAxis(4)->SetTitle("Mult");
130 fRecEventHist->Sumw2();
131
132 //Xv:Yv:Zv
133 Int_t binsMCEventHist[3]={80,80,100};
134 Double_t minMCEventHist[3]={-0.1,-0.1,-35.};
135 Double_t maxMCEventHist[3]={0.1,0.1,35.};
136 fMCEventHist = new THnSparseF("fMCEventHist","mcXv:mcYv:mcZv",3,binsMCEventHist,minMCEventHist,maxMCEventHist);
137 fMCEventHist->GetAxis(0)->SetTitle("mcXv (cm)");
138 fMCEventHist->GetAxis(1)->SetTitle("mcYv (cm)");
139 fMCEventHist->GetAxis(2)->SetTitle("mcZv (cm)");
140 fMCEventHist->Sumw2();
141
142 //Xv-mcXv:Yv-mcYv:Zv-mcZv:Mult
143 Int_t binsRecMCEventHist[4]={100,100,100,150};
144 Double_t minRecMCEventHist[4]={-1.0*kFact,-1.0*kFact,-1.0*kFact,0.};
145 Double_t maxRecMCEventHist[4]={1.0*kFact,1.0*kFact,1.0*kFact,150.};
146 fRecMCEventHist = new THnSparseF("fRecMCEventHist","Xv-mcXv:Yv-mcYv:Zv-mcZv:Mult",4,binsRecMCEventHist,minRecMCEventHist,maxRecMCEventHist);
147 fRecMCEventHist->GetAxis(0)->SetTitle("Xv-mcXv (cm)");
148 fRecMCEventHist->GetAxis(1)->SetTitle("Yv-mcYv (cm)");
149 fRecMCEventHist->GetAxis(2)->SetTitle("Zv-mcZv (cm)");
150 fRecMCEventHist->GetAxis(3)->SetTitle("Mult");
151 fRecMCEventHist->Sumw2();
152
153 //
154 // THnSparse track histograms
155 //
156
157 //nCrossRows:chi2PerClust:nCrossRows/nFindableClust:fracSharedClust:DCAy:DCAz:eta:phi:pt:hasStrangeMother:isFromMaterial:isPrim:charge
158 Int_t binsRecMCTrackHist[13]= { 160, 10, 20, 20, 50, 50, 20, 90, ptNbins, 2, 2, 2, 3 };
159 Double_t minRecMCTrackHist[13]={ 0., 0., 0., 0., -0.5,-0.5,-1.0, 0., ptMin, 0., 0., 0.,-1. };
160 Double_t maxRecMCTrackHist[13]={ 160., 10., 1., 1., 0.5, 0.5, 1.0, 2.*TMath::Pi(), ptMax, 2., 2., 2., 2. };
161
162 fRecMCTrackHist = new THnSparseF("fRecMCTrackHist","nCrossRows:chi2PerClust:nCrossRows/nFindableClust:fracSharedClust:DCAy:DCAz:eta:phi:pt:hasStrangeMother:isFromMaterial:isPrim:charge",13,binsRecMCTrackHist,minRecMCTrackHist,maxRecMCTrackHist);
163 fRecMCTrackHist->SetBinEdges(8,binsPt);
164 fRecMCTrackHist->GetAxis(0)->SetTitle("nCrossRows");
165 fRecMCTrackHist->GetAxis(1)->SetTitle("chi2PerClust");
166 fRecMCTrackHist->GetAxis(2)->SetTitle("nCrossRows/nFindableClust");
167 fRecMCTrackHist->GetAxis(3)->SetTitle("fracSharedClust");
168 fRecMCTrackHist->GetAxis(4)->SetTitle("DCAy (cm)");
169 fRecMCTrackHist->GetAxis(5)->SetTitle("DCAz (cm)");
170 fRecMCTrackHist->GetAxis(6)->SetTitle("#eta");
171 fRecMCTrackHist->GetAxis(7)->SetTitle("#phi (rad)");
172 fRecMCTrackHist->GetAxis(8)->SetTitle("p_{T} (GeV/c)");
173 fRecMCTrackHist->GetAxis(9)->SetTitle("hasStrangeMother");
174 fRecMCTrackHist->GetAxis(10)->SetTitle("isFromMaterial");
175 fRecMCTrackHist->GetAxis(11)->SetTitle("isPrim");
176 fRecMCTrackHist->GetAxis(12)->SetTitle("charge");
177 fRecMCTrackHist->Sumw2();
178
179 //nClust:chi2PerClust:nClust/nFindableClust:DCAy:DCAz:eta:phi:pt:kinkIdx:isPrim:polarity
180 /*
181 Int_t binsRecMCTrackHist[11]={160,80,80,100,100,90,90,ptNbins, 3, 2, 2};
182 Double_t minRecMCTrackHist[11]={0., 0., 0., -1.,-1.,-1.5, 0., ptMin, -1., 0., 0.};
183 Double_t maxRecMCTrackHist[11]={160.,10.,1.2, 1.,1.,1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2.};
184
185 fRecMCTrackHist = new THnSparseF("fRecMCTrackHist","nClust:chi2PerClust:nClust/nFindableClust:DCAy:DCAz:eta:phi:pt:kinkIdx:isPrim:polarity",11,binsRecMCTrackHist,minRecMCTrackHist,maxRecMCTrackHist);
186 fRecMCTrackHist->SetBinEdges(7,binsPt);
187
188 fRecMCTrackHist->GetAxis(0)->SetTitle("nClust");
189 fRecMCTrackHist->GetAxis(1)->SetTitle("chi2PerClust");
190 fRecMCTrackHist->GetAxis(2)->SetTitle("nClust/nFindableClust");
191 fRecMCTrackHist->GetAxis(3)->SetTitle("DCAy (cm)");
192 fRecMCTrackHist->GetAxis(4)->SetTitle("DCAz (cm)");
193 fRecMCTrackHist->GetAxis(5)->SetTitle("#eta");
194 fRecMCTrackHist->GetAxis(6)->SetTitle("#phi (rad)");
195 fRecMCTrackHist->GetAxis(7)->SetTitle("p_{T} (GeV/c)");
196 fRecMCTrackHist->GetAxis(8)->SetTitle("kinkIdx"); // 0 - no kink, -1 - kink mother, 1 - kink daugther
197 fRecMCTrackHist->GetAxis(9)->SetTitle("isPrim");
198 fRecMCTrackHist->GetAxis(10)->SetTitle("polarity");
199 fRecMCTrackHist->Sumw2();
200 */
201
202 // init output folder
203 fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
204
205}
206
207//_____________________________________________________________________________
208void AlidNdPtCutAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent)
209{
210 //
211 // Process real and/or simulated events
212 //
213 if(!esdEvent) {
214 AliDebug(AliLog::kError, "esdEvent not available");
215 return;
216 }
217
218 // get selection cuts
219 AlidNdPtEventCuts *evtCuts = GetEventCuts();
220 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
221 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
222
223 if(!evtCuts || !accCuts || !esdTrackCuts) {
224 AliDebug(AliLog::kError, "cuts not available");
225 return;
226 }
227
228 // trigger selection
229
230 Bool_t isEventTriggered = kTRUE;
231 AliPhysicsSelection *physicsSelection = NULL;
232 AliTriggerAnalysis* triggerAnalysis = NULL;
233
234 //
235 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
236 if (!inputHandler)
237 {
238 Printf("ERROR: Could not receive input handler");
239 return;
240 }
241
242 if(evtCuts->IsTriggerRequired())
243 {
244 // always MB
245 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
246
247 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
248 if(!physicsSelection) return;
249 //SetPhysicsTriggerSelection(physicsSelection);
250
251 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
252 // set trigger (V0AND)
253 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
254 if(!triggerAnalysis) return;
255 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
256 }
257 }
258
259 // use MC information
260 AliHeader* header = 0;
261 AliGenEventHeader* genHeader = 0;
262 AliStack* stack = 0;
263 TArrayF vtxMC(3);
264 AliPWG0Helper::MCProcessType evtType = AliPWG0Helper::kInvalidProcess;
265
266 if(IsUseMCInfo())
267 {
268 if(!mcEvent) {
269 AliDebug(AliLog::kError, "mcEvent not available");
270 return;
271 }
272
273 // get MC event header
274 header = mcEvent->Header();
275 if (!header) {
276 AliDebug(AliLog::kError, "Header not available");
277 return;
278 }
279
280 // MC particle stack
281 stack = mcEvent->Stack();
282 if (!stack) {
283 AliDebug(AliLog::kError, "Stack not available");
284 return;
285 }
286
287 // get event type (ND=0x1, DD=0x2, SD=0x4)
288 evtType = AliPWG0Helper::GetEventProcessType(header);
289 AliDebug(AliLog::kDebug+1, Form("Found process type %d", evtType));
290
291 // get MC vertex
292 genHeader = header->GenEventHeader();
293 if (!genHeader) {
294 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
295 return;
296 }
297 genHeader->PrimaryVertex(vtxMC);
298
299 // Fill MC event histogram
300 Double_t vMCEventHist[3]={vtxMC[0],vtxMC[1],vtxMC[2]};
301 fMCEventHist->Fill(vMCEventHist);
302
303 } // end bUseMC
304
305 // get reconstructed vertex
306 Bool_t bRedoTPCVertex = evtCuts->IsRedoTPCVertex();
307 Bool_t bUseConstraints = evtCuts->IsUseBeamSpotConstraint();
308 const AliESDVertex* vtxESD = AlidNdPtHelper::GetVertex(esdEvent,evtCuts,accCuts,esdTrackCuts,GetAnalysisMode(),kFALSE,bRedoTPCVertex,bUseConstraints);
309 Bool_t isRecVertex = AlidNdPtHelper::TestRecVertex(vtxESD, esdEvent->GetPrimaryVertexSPD(), GetAnalysisMode(), kFALSE);
310 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD) && isRecVertex;
311
312 TObjArray *allChargedTracks=0;
313 Int_t multAll=0;
314
315 //
316 // event counter
317 //
318 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK,isEventTriggered);
319
320 Bool_t isTrigAndVertex = isEventTriggered && isEventOK;
2942f542 321 Double_t vEventCount[2] = { static_cast<Double_t>(isEventTriggered), static_cast<Double_t>(isTrigAndVertex)};
a65a7e70 322 fEventCount->Fill(vEventCount);
323
324 //
325 // cosmic background and splitted tracks
326 //
327 if(GetParticleMode() == AlidNdPtHelper::kBackgroundTrack)
328 {
329 AlidNdPtBackgroundCuts *backCuts = GetBackgroundCuts();
330 if(!backCuts) return;
331
332 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
333 {
334 AliESDtrack *track1 = esdEvent->GetTrack(iTrack);
335 if(!track1) continue;
336 if(track1->Charge()==0) continue;
337
338 for (Int_t jTrack = iTrack+1; jTrack < esdEvent->GetNumberOfTracks(); jTrack++)
339 {
340 AliESDtrack *track2 = esdEvent->GetTrack(jTrack);
341 if(!track2) continue;
342 if(track2->Charge()==0) continue;
343
344 //printf("track2->Charge() %d\n",track2->Charge());
345
346 backCuts->IsBackgroundTrack(track1, track2);
347 }
348 }
349 }
350
351 // check event cuts
352 if(isEventOK && isEventTriggered)
353 {
354 // get all charged tracks
355 allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
356 if(!allChargedTracks) return;
357
358 Int_t entries = allChargedTracks->GetEntries();
359 for(Int_t i=0; i<entries;++i)
360 {
361 AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
362 if(!track) continue;
363
364 if(!esdTrackCuts->AcceptTrack(track)) continue;
365
366 //
367 Bool_t isOK = kFALSE;
368 Double_t x[3]; track->GetXYZ(x);
369 Double_t b[3]; AliTracker::GetBxByBz(x,b);
370
371 //
372 // if TPC-ITS hybrid tracking (kTPCITSHybrid)
373 // replace track parameters with TPC-ony track parameters
374 //
375 if( GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybrid || GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtx || GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtxDCArPt)
376 {
377 // Relate TPC-only tracks to SPD vertex
378 isOK = track->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
379 if(!isOK) continue;
380
381 // replace esd track parameters with TPCinner
382 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
383 if (!tpcTrack) return;
384 track->Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),tpcTrack->GetParameter(),tpcTrack->GetCovariance());
385
386 if(tpcTrack) delete tpcTrack;
387 }
388
389 //
390 if (GetAnalysisMode()==AlidNdPtHelper::kTPCSPDvtxUpdate || GetAnalysisMode() == AlidNdPtHelper::kTPCTrackSPDvtxUpdate)
391 {
392 //
393 // update track parameters
394 //
395 AliExternalTrackParam cParam;
396 isOK = track->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig, &cParam);
397 if(!isOK) continue;
398
399 track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
400 }
401
402 FillHistograms(track, stack);
403 multAll++;
404 }
405
e690d4d0 406 Double_t vRecEventHist[5] = {vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(),vtxESD->GetZRes(),static_cast<Double_t>(multAll)};
a65a7e70 407 fRecEventHist->Fill(vRecEventHist);
408
409 if(IsUseMCInfo()) {
e690d4d0 410 Double_t vRecMCEventHist[5] = {vtxESD->GetX()-vtxMC[0],vtxESD->GetY()-vtxMC[1],vtxESD->GetZ()-vtxMC[2],static_cast<Double_t>(multAll)};
a65a7e70 411 fRecMCEventHist->Fill(vRecMCEventHist);
412 }
413 }
414
415 if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
416
417}
418
419//_____________________________________________________________________________
420void AlidNdPtCutAnalysis::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack) const
421{
422 //
423 // Fill ESD track and MC histograms
424 //
425 if(!esdTrack) return;
426 if(esdTrack->Charge() == 0.) return;
427
428 Float_t pt = esdTrack->Pt();
429 Float_t eta = esdTrack->Eta();
430 Float_t phi = esdTrack->Phi();
431
432
433 Int_t nClust = 0;
434 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
435 nClust = esdTrack->GetTPCNclsIter1();
436 } else {
437 nClust = esdTrack->GetTPCclusters(0);
438 }
439
440 Float_t chi2PerCluster = 0.;
441 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
442 if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2Iter1()/Float_t(nClust);
443 } else {
444 chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
445 }
446
447 Int_t nFindableClust = esdTrack->GetTPCNclsF();
448
449
450 Float_t clustPerFindClust = 0.;
451 if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
452
453 Float_t b[2], bCov[3];
454 esdTrack->GetImpactParameters(b,bCov);
455
456 //
457 Float_t nCrossedRowsTPC = esdTrack->GetTPCClusterInfo(2,1);
458
459 Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
460 if (esdTrack->GetTPCNclsF()>0) {
461 ratioCrossedRowsOverFindableClustersTPC = esdTrack->GetTPCClusterInfo(2,1)/esdTrack->GetTPCNclsF();
462 }
463
464 //
465 Int_t nClustersTPCShared = esdTrack->GetTPCnclsS();
466 Float_t fracClustersTPCShared = -1.;
467 fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClust);
468
469
470 // kink idx
471 Int_t kinkIdx = 0;
472 //if(esdTrack->GetKinkIndex(0) > 0.) isKink = kTRUE;
473 if(esdTrack->GetKinkIndex(0) > 0) kinkIdx = 1; // kink daughter
474 else if(esdTrack->GetKinkIndex(0) < 0) kinkIdx = -1; // kink mother
475 else kinkIdx = 0; // not kink
476
477 //printf("esdTrack->GetKinkIndex(0) %d \n", esdTrack->GetKinkIndex(0));
478 //printf("esdTrack->GetKinkIndex(1) %d \n", esdTrack->GetKinkIndex(1));
479 //printf("esdTrack->GetKinkIndex(2) %d \n", esdTrack->GetKinkIndex(2));
480 //printf("kinkIdx %d \n", kinkIdx);
481
482 //
483 // Fill rec vs MC information
484 //
485 Bool_t isPrim = kTRUE;
486 Bool_t hasStrangeMother = kFALSE;
487 Bool_t isFromMaterial = kFALSE;
488
489 if(IsUseMCInfo()) {
490 if(!stack) return;
491 Int_t label = TMath::Abs(esdTrack->GetLabel());
492 TParticle* particle = stack->Particle(label);
493 if(!particle) return;
494 if(particle->GetPDG() && particle->GetPDG()->Charge()==0.) return;
495 isPrim = stack->IsPhysicalPrimary(label);
496
497 // check whether has stange mother
498 //
499 Int_t motherPdg = -1;
500 TParticle* mother = 0;
501
502 Int_t motherLabel = particle->GetMother(0);
503 if(motherLabel>0) mother = stack->Particle(motherLabel);
504 if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
505 Int_t mech = particle->GetUniqueID(); // production mechanism
506
507 if( (motherPdg == 3122) || (motherPdg == -3122) || (motherPdg == 310)) // lambda, antilambda, k0s
508 {
509 if( (mech == 4) || (mech == 5) ) hasStrangeMother = kTRUE;
510 }
511 else {
512 //if(isPrim==0 && mech == 13)
513 //printf("mech %d \n", mech);
514 if(!isPrim) isFromMaterial = kTRUE;
515 }
516
517 //if(isPrim && pt > 1.5 && kinkIdx == -1) printf("nClust %d \n", nClust);
518 }
519
520 // fill histo
521 Int_t charge = esdTrack->Charge();
522
523 //Double_t vRecMCTrackHist[11] = { nClust,chi2PerCluster,clustPerFindClust,b[0],b[1],eta,phi,pt,kinkIdx,isPrim, polarity };
524 //fRecMCTrackHist->Fill(vRecMCTrackHist);
525
2942f542 526 Double_t vRecMCTrackHist[13] = { static_cast<Double_t>(nCrossedRowsTPC), chi2PerCluster, ratioCrossedRowsOverFindableClustersTPC, fracClustersTPCShared , b[0], b[1], eta, phi, pt, static_cast<Double_t>(hasStrangeMother), static_cast<Double_t>(isFromMaterial), static_cast<Double_t>(isPrim), static_cast<Double_t>(charge) };
a65a7e70 527 fRecMCTrackHist->Fill(vRecMCTrackHist);
528
529}
530
531//_____________________________________________________________________________
532Long64_t AlidNdPtCutAnalysis::Merge(TCollection* const list)
533{
534 // Merge list of objects (needed by PROOF)
535
536 if (!list)
537 return 0;
538
539 if (list->IsEmpty())
540 return 1;
541
542 TIterator* iter = list->MakeIterator();
543 TObject* obj = 0;
544
545 //TList *collPhysSelection = new TList;
546
547 // collection of generated histograms
548 Int_t count=0;
549 while((obj = iter->Next()) != 0) {
550 AlidNdPtCutAnalysis* entry = dynamic_cast<AlidNdPtCutAnalysis*>(obj);
551 if (entry == 0) continue;
552
553 // event histo
554 fEventCount->Add(entry->fEventCount);
555 fRecEventHist->Add(entry->fRecEventHist);
556 fRecMCEventHist->Add(entry->fRecMCEventHist);
557 fMCEventHist->Add(entry->fMCEventHist);
558
559 // track histo
560 fRecMCTrackHist->Add(entry->fRecMCTrackHist);
561
562 // physics selection
563 //collPhysSelection->Add(entry->GetPhysicsTriggerSelection());
564
565 count++;
566 }
567
568 //AliPhysicsSelection *trigSelection = GetPhysicsTriggerSelection();
569 //trigSelection->Merge(collPhysSelection);
570
571 //if(collPhysSelection) delete collPhysSelection;
572
573return count;
574}
575
576//_____________________________________________________________________________
577void AlidNdPtCutAnalysis::Analyse()
578{
579 //
580 // Analyse histograms
581 //
582 TH1::AddDirectory(kFALSE);
583 TObjArray *aFolderObj = new TObjArray;
584 if(!aFolderObj) return;
585
586 TH1D *h1D = 0;
587 TH2D *h2D = 0;
588
589
590 //
591 // get cuts
592 //
593 AlidNdPtEventCuts *evtCuts = GetEventCuts();
594 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
595 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
596
597 if(!evtCuts || !accCuts || !esdTrackCuts) {
598 Error("AlidNdPtCutAnalysis::Analyse()", "cuts not available");
599 return;
600 }
601
602 //
603 // set min and max values
604 //
605 Double_t minPt = accCuts->GetMinPt();
606 Double_t maxPt = accCuts->GetMaxPt();
607 Double_t minEta = accCuts->GetMinEta();
608 Double_t maxEta = accCuts->GetMaxEta()-0.00001;
609
610 Double_t maxDCAr = accCuts->GetMaxDCAr();
611
612 //
613 // Event counters
614 //
615 h2D = (TH2D*)fEventCount->Projection(0,1);
616 if(!h2D) return;
617 h2D->SetName("trig_vs_trigANDvertex");
618 aFolderObj->Add(h2D);
619
620 fEventCount->GetAxis(0)->SetRange(2,2); // triggered
621 h1D = (TH1D*)fEventCount->Projection(1);
622 if(!h1D) return;
623 h1D->SetTitle("rec. vertex for triggered events");
624 h1D->SetName("trigANDvertex");
625 aFolderObj->Add(h1D);
626
627 //
628 // Create rec. event histograms
629 //
630 h1D = (TH1D *)fRecEventHist->Projection(0);
631 if(!h1D) return;
632 h1D->SetName("rec_xv");
633 aFolderObj->Add(h1D);
634
635 h1D = (TH1D *)fRecEventHist->Projection(1);
636 if(!h1D) return;
637 h1D->SetName("rec_yv");
638 aFolderObj->Add(h1D);
639
640 h1D = (TH1D *)fRecEventHist->Projection(2);
641 if(!h1D) return;
642 h1D->SetName("rec_zv");
643 aFolderObj->Add(h1D);
644
645 h2D = (TH2D *)fRecEventHist->Projection(3,4);
646 if(!h2D) return;
647 h2D->SetName("rec_resZv_vs_Mult");
648 aFolderObj->Add(h2D);
649
650 h2D = (TH2D *)fRecEventHist->Projection(0,1);
651 if(!h2D) return;
652 h2D->SetName("rec_xv_vs_yv");
653 aFolderObj->Add(h2D);
654
655 h2D = (TH2D *)fRecEventHist->Projection(0,2);
656 if(!h2D) return;
657 h2D->SetName("rec_xv_vs_zv");
658 aFolderObj->Add(h2D);
659
660 h2D = (TH2D *)fRecEventHist->Projection(3,4);
661 if(!h2D) return;
662 h2D->SetName("rec_resZv_vs_Mult");
663 aFolderObj->Add(h2D);
664
665 //
666 // MC available
667 //
668 if(IsUseMCInfo()) {
669
670 //
671 // Create mc event histograms
672 //
673 h2D = (TH2D *)fMCEventHist->Projection(0,1);
674 if(!h2D) return;
675 h2D->SetName("mc_xv_vs_yv");
676 aFolderObj->Add(h2D);
677
678 h2D = (TH2D *)fMCEventHist->Projection(0,2);
679 if(!h2D) return;
680 h2D->SetName("mc_xv_vs_zv");
681 aFolderObj->Add(h2D);
682
683 //
684 // Create rec-mc event histograms
685 //
686 h2D = (TH2D *)fRecMCEventHist->Projection(0,3);
687 if(!h2D) return;
688 h2D->SetName("rec_mc_deltaXv_vs_mult");
689 aFolderObj->Add(h2D);
690
691 h2D = (TH2D *)fRecMCEventHist->Projection(1,3);
692 if(!h2D) return;
693 h2D->SetName("rec_mc_deltaYv_vs_mult");
694 aFolderObj->Add(h2D);
695
696 h2D = (TH2D *)fRecMCEventHist->Projection(2,3);
697 if(!h2D) return;
698 h2D->SetName("rec_mc_deltaZv_vs_mult");
699 aFolderObj->Add(h2D);
700
701 } // end use MC info
702
703
704
705 //
706 // Create rec-mc track track histograms
707 //
708
709 // DCA cuts
710 fRecMCTrackHist->GetAxis(3)->SetRangeUser(-maxDCAr,maxDCAr);
711 fRecMCTrackHist->GetAxis(4)->SetRangeUser(-maxDCAr,maxDCAr);
712
713 h2D = (TH2D *)fRecMCTrackHist->Projection(7,5);
714 if(!h2D) return;
715 h2D->SetName("pt_vs_eta");
716 aFolderObj->Add(h2D);
717
718 fRecMCTrackHist->GetAxis(7)->SetRangeUser(minPt,maxPt);
719
720 h2D = (TH2D *)fRecMCTrackHist->Projection(0,5);
721 if(!h2D) return;
722 h2D->SetName("nClust_vs_eta");
723 aFolderObj->Add(h2D);
724
725 h2D = (TH2D *)fRecMCTrackHist->Projection(1,5);
726 if(!h2D) return;
727 h2D->SetName("chi2PerClust_vs_eta");
728 aFolderObj->Add(h2D);
729
730 h2D = (TH2D *)fRecMCTrackHist->Projection(2,5);
731 if(!h2D) return;
732 h2D->SetName("ratio_nClust_nFindableClust_vs_eta");
733 aFolderObj->Add(h2D);
734
735 h2D = (TH2D *)fRecMCTrackHist->Projection(5,6);
736 if(!h2D) return;
737 h2D->SetName("eta_vs_phi");
738 aFolderObj->Add(h2D);
739
740 //
741 fRecMCTrackHist->GetAxis(5)->SetRangeUser(minEta,maxEta);
742
743 h2D = (TH2D *)fRecMCTrackHist->Projection(0,6);
744 if(!h2D) return;
745 h2D->SetName("nClust_vs_phi");
746 aFolderObj->Add(h2D);
747
748 h2D = (TH2D *)fRecMCTrackHist->Projection(1,6);
749 if(!h2D) return;
750 h2D->SetName("chi2PerClust_vs_phi");
751 aFolderObj->Add(h2D);
752
753 h2D = (TH2D *)fRecMCTrackHist->Projection(2,6);
754 if(!h2D) return;
755 h2D->SetName("ratio_nClust_nFindableClust_vs_phi");
756 aFolderObj->Add(h2D);
757
758 //
759 fRecMCTrackHist->GetAxis(7)->SetRangeUser(0.0,maxPt);
760
761 h2D = (TH2D *)fRecMCTrackHist->Projection(0,7);
762 if(!h2D) return;
763 h2D->SetName("nClust_vs_pt");
764 aFolderObj->Add(h2D);
765
766 h2D = (TH2D *)fRecMCTrackHist->Projection(1,7);
767 if(!h2D) return;
768 h2D->SetName("chi2PerClust_vs_pt");
769 aFolderObj->Add(h2D);
770
771 h2D = (TH2D *)fRecMCTrackHist->Projection(2,7);
772 if(!h2D) return;
773 h2D->SetName("ratio_nClust_nFindableClust_vs_pt");
774 aFolderObj->Add(h2D);
775
776 h2D = (TH2D *)fRecMCTrackHist->Projection(6,7);
777 if(!h2D) return;
778 h2D->SetName("phi_vs_pt");
779 aFolderObj->Add(h2D);
780
781
782 // fiducial volume
783 fRecMCTrackHist->GetAxis(5)->SetRangeUser(minEta,maxEta);
784 fRecMCTrackHist->GetAxis(7)->SetRangeUser(minPt,maxPt);
785
786 // DCA cuts
787 fRecMCTrackHist->GetAxis(3)->SetRangeUser(-maxDCAr,maxDCAr);
788 fRecMCTrackHist->GetAxis(4)->SetRangeUser(-maxDCAr,maxDCAr);
789
790 h2D = (TH2D *)fRecMCTrackHist->Projection(0,1);
791 if(!h2D) return;
792 h2D->SetName("nClust_vs_chi2PerClust");
793 aFolderObj->Add(h2D);
794
795 h2D = (TH2D *)fRecMCTrackHist->Projection(0,2);
796 if(!h2D) return;
797 h2D->SetName("nClust_vs_ratio_nClust_nFindableClust");
798 aFolderObj->Add(h2D);
799
800 //
801 // DCAy cuts
802 //
803 fRecMCTrackHist->GetAxis(0)->SetRange(50,160); // nClust/track > 50
804 fRecMCTrackHist->GetAxis(1)->SetRangeUser(0.,3.9999); // chi2/cluster < 4.0
805 fRecMCTrackHist->GetAxis(3)->SetRange(1,fRecMCTrackHist->GetAxis(3)->GetNbins());
806 //fRecMCTrackHist->GetAxis(4)->SetRangeUser(-1.0,1.0);
807 fRecMCTrackHist->GetAxis(4)->SetRange(1,fRecMCTrackHist->GetAxis(4)->GetNbins());
808
809 // sec
810 fRecMCTrackHist->GetAxis(9)->SetRange(1,1);
811 h1D = (TH1D *)fRecMCTrackHist->Projection(3);
812 if(!h1D) return;
813 h1D->SetName("dcay_sec");
814 aFolderObj->Add(h1D);
815
816 // prim
817 fRecMCTrackHist->GetAxis(9)->SetRange(2,2);
818 h1D = (TH1D *)fRecMCTrackHist->Projection(3);
819 if(!h1D) return;
820 h1D->SetName("dcay_prim");
821 aFolderObj->Add(h1D);
822
823 // DCAz cuts
824 //fRecMCTrackHist->GetAxis(3)->SetRangeUser(-1.0,1.0);
825 fRecMCTrackHist->GetAxis(4)->SetRange(1,fRecMCTrackHist->GetAxis(4)->GetNbins());
826
827 // sec
828 fRecMCTrackHist->GetAxis(9)->SetRange(1,1);
829 h1D = (TH1D *)fRecMCTrackHist->Projection(4);
830 if(!h1D) return;
831 h1D->SetName("dcaz_sec");
832 aFolderObj->Add(h1D);
833
834 // prim
835 fRecMCTrackHist->GetAxis(9)->SetRange(2,2);
836 h1D = (TH1D *)fRecMCTrackHist->Projection(4);
837 if(!h1D) return;
838 h1D->SetName("dcaz_prim");
839 aFolderObj->Add(h1D);
840
841
842 // export objects to analysis folder
843 fAnalysisFolder = ExportToFolder(aFolderObj);
844 if(!fAnalysisFolder) {
845 if(aFolderObj) delete aFolderObj;
846 return;
847 }
848
849 // delete only TObjArray
850 if(aFolderObj) delete aFolderObj;
851}
852
853//_____________________________________________________________________________
854TFolder* AlidNdPtCutAnalysis::ExportToFolder(TObjArray * const array)
855{
856 // recreate folder avery time and export objects to new one
857 //
858 AlidNdPtCutAnalysis * comp=this;
859 TFolder *folder = comp->GetAnalysisFolder();
860
861 TString name, title;
862 TFolder *newFolder = 0;
863 Int_t i = 0;
864 Int_t size = array->GetSize();
865
866 if(folder) {
867 // get name and title from old folder
868 name = folder->GetName();
869 title = folder->GetTitle();
870
871 // delete old one
872 delete folder;
873
874 // create new one
875 newFolder = CreateFolder(name.Data(),title.Data());
876 newFolder->SetOwner();
877
878 // add objects to folder
879 while(i < size) {
880 newFolder->Add(array->At(i));
881 i++;
882 }
883 }
884
885return newFolder;
886}
887
888//_____________________________________________________________________________
889TFolder* AlidNdPtCutAnalysis::CreateFolder(TString name,TString title) {
890// create folder for analysed histograms
891//
892TFolder *folder = 0;
893 folder = new TFolder(name.Data(),title.Data());
894
895 return folder;
896}