]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/dNdPt/AlidNdPtHelper.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
dNdPt analysis update for preliminary spectra
[u/mrichter/AliRoot.git] / PWG0 / dNdPt / AlidNdPtHelper.cxx
CommitLineData
0aaa8b91 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 **************************************************************************/
f2dec884 15//
16// static dNdPt helper functions
17//
18// basic functionality to select events and tracks
19// for dNdPt analysis
20//
21// Origin: Jan Fiete Grosse-Oetringhaus
22// Modified and Extended: Jacek Otwinowski 19/11/2009
23//
24//
0aaa8b91 25
26#include <TROOT.h>
f2dec884 27#include <TCanvas.h>
28#include <TF1.h>
0aaa8b91 29#include <TH1.h>
30#include <TH2.h>
31#include <TH3.h>
0aaa8b91 32
33#include <AliHeader.h>
34#include <AliStack.h>
35#include <AliLog.h>
0aaa8b91 36#include <AliESD.h>
37#include <AliESDEvent.h>
38#include <AliMCEvent.h>
39#include <AliESDVertex.h>
40#include <AliVertexerTracks.h>
0aaa8b91 41#include <AliMathBase.h>
42#include <AliESDtrackCuts.h>
774986ef 43#include <AliTracker.h>
0aaa8b91 44#include "dNdPt/AlidNdPtEventCuts.h"
45#include "dNdPt/AlidNdPtAcceptanceCuts.h"
46#include "dNdPt/AlidNdPtHelper.h"
47
48//____________________________________________________________________
49ClassImp(AlidNdPtHelper)
50
0aaa8b91 51//____________________________________________________________________
f2dec884 52const AliESDVertex* AlidNdPtHelper::GetVertex(AliESDEvent* const aEsd, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts, AliESDtrackCuts *const trackCuts, AnalysisMode analysisMode, Bool_t debug, Bool_t bRedoTPC, Bool_t bUseMeanVertex)
0aaa8b91 53{
54 // Get the vertex from the ESD and returns it if the vertex is valid
55 //
56 // Second argument decides which vertex is used (this selects
57 // also the quality criteria that are applied)
58
59 if(!aEsd)
60 {
61 ::Error("AlidNdPtHelper::GetVertex()","esd event is NULL");
62 return NULL;
63 }
64
65 if(!evtCuts || !accCuts || !trackCuts)
66 {
67 ::Error("AlidNdPtHelper::GetVertex()","cuts not available");
68 return NULL;
69 }
70
71 const AliESDVertex* vertex = 0;
72 AliESDVertex *initVertex = 0;
774986ef 73 if (analysisMode == kSPD || analysisMode == kTPCITS ||
74 analysisMode == kTPCSPDvtx || analysisMode == kTPCSPDvtxUpdate || analysisMode == kTPCITSHybrid)
0aaa8b91 75 {
76 vertex = aEsd->GetPrimaryVertexSPD();
77 if (debug)
78 Printf("AlidNdPtHelper::GetVertex: Returning SPD vertex");
79 }
791aaf54 80 else if (analysisMode == kTPCTrackSPDvtx || analysisMode == kTPCTrackSPDvtxUpdate || analysisMode == kTPCITSHybridTrackSPDvtx || analysisMode == kTPCITSHybridTrackSPDvtxDCArPt)
81 {
82 vertex = aEsd->GetPrimaryVertexTracks();
83 if(!vertex) return NULL;
84 if(vertex->GetNContributors()<1) {
85 // SPD vertex
86 vertex = aEsd->GetPrimaryVertexSPD();
87 }
88 }
847e74b2 89 else if (analysisMode == kTPC)
0aaa8b91 90 {
91 if(bRedoTPC) {
92
93 Double_t kBz = aEsd->GetMagneticField();
94 AliVertexerTracks vertexer(kBz);
95
96 if(bUseMeanVertex) {
97 Double_t pos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
98 Double_t err[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
0aaa8b91 99 initVertex = new AliESDVertex(pos,err);
100 vertexer.SetVtxStart(initVertex);
101 vertexer.SetConstraintOn();
102 }
103
0aaa8b91 104 Double_t maxDCAr = accCuts->GetMaxDCAr();
105 Double_t maxDCAz = accCuts->GetMaxDCAz();
106 Int_t minTPCClust = trackCuts->GetMinNClusterTPC();
107
bad4ba69 108 //vertexer.SetTPCMode(Double_t dcacut=0.1, Double_t dcacutIter0=1.0, Double_t maxd0z0=5.0, Int_t minCls=10, Int_t mintrks=1, Double_t nsigma=3., Double_t mindetfitter=0.1, Double_t maxtgl=1.5, Double_t fidR=3., Double_t fidZ=30., Int_t finderAlgo=1, Int_t finderAlgoIter0=4);
0aaa8b91 109 vertexer.SetTPCMode(0.1,1.0,5.0,minTPCClust,1,3.,0.1,2.0,maxDCAr,maxDCAz,1,4);
110
111 // TPC track preselection
112 Int_t ntracks = aEsd->GetNumberOfTracks();
113 TObjArray array(ntracks);
114 UShort_t *id = new UShort_t[ntracks];
115
116 Int_t count=0;
117 for (Int_t i=0;i <ntracks; i++) {
118 AliESDtrack *t = aEsd->GetTrack(i);
119 if (!t) continue;
120 if (t->Charge() == 0) continue;
121 if (!t->GetTPCInnerParam()) continue;
122 if (t->GetTPCNcls()<vertexer.GetMinClusters()) continue;
123 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
124 if(tpcTrack) {
125 array.AddLast(tpcTrack);
126 id[count] = (UShort_t)t->GetID();
127 count++;
128 }
129 }
130 AliESDVertex *vTPC = vertexer.VertexForSelectedTracks(&array,id, kTRUE, kTRUE, bUseMeanVertex);
bad4ba69 131
132 // set recreated TPC vertex
0aaa8b91 133 aEsd->SetPrimaryVertexTPC(vTPC);
134
135 for (Int_t i=0; i<aEsd->GetNumberOfTracks(); i++) {
136 AliESDtrack *t = aEsd->GetTrack(i);
774986ef 137 if(!t) continue;
138
139 Double_t x[3]; t->GetXYZ(x);
140 Double_t b[3]; AliTracker::GetBxByBz(x,b);
141 t->RelateToVertexTPCBxByBz(vTPC, b, kVeryBig);
0aaa8b91 142 }
143
144 delete vTPC;
145 array.Delete();
146 delete [] id; id=NULL;
147
148 }
149 vertex = aEsd->GetPrimaryVertexTPC();
150 if (debug)
774986ef 151 Printf("AlidNdPtHelper::GetVertex: Returning vertex from tracks");
152 }
153 else
154 Printf("AlidNdPtHelper::GetVertex: ERROR: Invalid second argument %d", analysisMode);
0aaa8b91 155
156 if (!vertex) {
157 if (debug)
158 Printf("AlidNdPtHelper::GetVertex: No vertex found in ESD");
159 return 0;
160 }
161
162 if (debug)
163 {
164 Printf("AlidNdPtHelper::GetVertex: Returning valid vertex: %s", vertex->GetTitle());
165 vertex->Print();
166 }
167
168 if(initVertex) delete initVertex; initVertex=NULL;
169 return vertex;
170}
171
bad4ba69 172//____________________________________________________________________
791aaf54 173Bool_t AlidNdPtHelper::TestRecVertex(const AliESDVertex* vertex, const AliESDVertex* vertexSPD, AnalysisMode analysisMode, Bool_t debug)
bad4ba69 174{
175 // Checks if a vertex meets the needed quality criteria
176 if(!vertex) return kFALSE;
774986ef 177 if(!vertex->GetStatus()) return kFALSE;
bad4ba69 178
179 Float_t requiredZResolution = -1;
774986ef 180 if (analysisMode == kSPD || analysisMode == kTPCITS ||
791aaf54 181 analysisMode == kTPCSPDvtx || analysisMode == kTPCSPDvtxUpdate || analysisMode == kTPCITSHybrid ||
182 analysisMode == kTPCTrackSPDvtx || analysisMode == kTPCTrackSPDvtxUpdate || analysisMode == kTPCITSHybridTrackSPDvtx || analysisMode == kTPCITSHybridTrackSPDvtxDCArPt
183 )
bad4ba69 184 {
774986ef 185 requiredZResolution = 1000;
bad4ba69 186 }
847e74b2 187 else if (analysisMode == kTPC)
bad4ba69 188 requiredZResolution = 10.;
189
774986ef 190 // check resolution
191 Double_t zRes = vertex->GetZRes();
192
193 if (zRes > requiredZResolution) {
194 if (debug)
195 Printf("AlidNdPtHelper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
196 return kFALSE;
197 }
791aaf54 198
199 // always check for SPD vertex
200 if(!vertexSPD) return kFALSE;
201 if(!vertexSPD->GetStatus()) return kFALSE;
202 if (vertexSPD->IsFromVertexerZ())
774986ef 203 {
791aaf54 204 if (vertexSPD->GetDispersion() > 0.02)
774986ef 205 {
206 if (debug)
207 Printf("AliPWG0Helper::TestVertex: Delta Phi too large in Vertexer Z: %f (required: %f", vertex->GetDispersion(), 0.02);
208 return kFALSE;
209 }
210 }
211
212 /*
bad4ba69 213 // check Ncontributors
214 if (vertex->GetNContributors() <= 0) {
215 if (debug){
216 Printf("AlidNdPtHelper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors());
217 Printf("AlidNdPtHelper::GetVertex: NIndices(): %d",vertex->GetNIndices());
218 vertex->Print();
219 }
220 return kFALSE;
221 }
774986ef 222 */
bad4ba69 223
224 return kTRUE;
225}
226
791aaf54 227//____________________________________________________________________
228Bool_t AlidNdPtHelper::IsGoodImpPar(AliESDtrack *const track)
229{
230//
231// check whether particle has good DCAr(Pt) impact
232// parameter. Only for TPC+ITS tracks (7*sigma cut)
233// Origin: Andrea Dainese
234//
235
236Float_t d0z0[2],covd0z0[3];
237track->GetImpactParameters(d0z0,covd0z0);
238Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
239Float_t d0max = 7.*sigma;
240if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
241
242return kFALSE;
243}
244
0aaa8b91 245//____________________________________________________________________
f2dec884 246Bool_t AlidNdPtHelper::IsPrimaryParticle(AliStack* const stack, Int_t idx, ParticleMode particleMode)
0aaa8b91 247{
774986ef 248//
0aaa8b91 249// check primary particles
847e74b2 250// depending on the particle mode
0aaa8b91 251//
252 if(!stack) return kFALSE;
253
254 TParticle* particle = stack->Particle(idx);
255 if (!particle) return kFALSE;
256
257 // only charged particles
258 Double_t charge = particle->GetPDG()->Charge()/3.;
f2dec884 259 if (TMath::Abs(charge) < 0.001) return kFALSE;
0aaa8b91 260
261 Int_t pdg = TMath::Abs(particle->GetPdgCode());
262
263 // physical primary
264 Bool_t prim = stack->IsPhysicalPrimary(idx);
265
847e74b2 266 if(particleMode==kMCPion) {
0aaa8b91 267 if(prim && pdg==kPiPlus) return kTRUE;
268 else return kFALSE;
269 }
270
847e74b2 271 if (particleMode==kMCKaon) {
0aaa8b91 272 if(prim && pdg==kKPlus) return kTRUE;
273 else return kFALSE;
274 }
275
847e74b2 276 if (particleMode==kMCProton) {
0aaa8b91 277 if(prim && pdg==kProton) return kTRUE;
278 else return kFALSE;
279 }
280
281return prim;
282}
283
d269b0e6 284//____________________________________________________________________
f2dec884 285Bool_t AlidNdPtHelper::IsCosmicTrack(AliESDtrack *const track1, AliESDtrack *const track2)
d269b0e6 286{
287//
288// check cosmic tracks
289//
290 if(!track1) return kFALSE;
291 if(!track2) return kFALSE;
d269b0e6 292
774986ef 293 //
294 // cosmic tracks in TPC
295 //
296 //if( TMath::Abs( track1->Theta() - track2->Theta() ) < 0.004 &&
297 // ((TMath::Abs(track1->Phi()-track2->Phi()) - TMath::Pi() )<0.004) &&
298 //if( track1->Pt() > 4.0 && (TMath::Abs(track1->Phi()-track2->Phi())-TMath::Pi())<0.1
299 // && (TMath::Abs(track1->Eta()+track2->Eta())-0.1) < 0.0 && (track1->Charge()+track2->Charge()) == 0)
300
301 //Float_t scaleF= 6.0;
302 if ( track1->Pt() > 4 && track2->Pt() > 4 &&
303 //(TMath::Abs(track1->GetSnp()-track2->GetSnp())-2.) < scaleF * TMath::Sqrt(track1->GetSigmaSnp2()+track2->GetSigmaSnp2()) &&
304 //TMath::Abs(track1->GetTgl()-track2->GetTgl()) < scaleF * TMath::Sqrt(track1->GetSigmaTgl2()+track2->GetSigmaTgl2()) &&
305 //TMath::Abs(track1->OneOverPt()-track2->OneOverPt()) < scaleF * TMath::Sqrt(track1->GetSigma1Pt2()+track2->GetSigma1Pt2()) &&
306 (track1->Charge()+track2->Charge()) == 0 &&
307 track1->Eta()*track2->Eta()<0.0 && TMath::Abs(track1->Eta()+track2->Eta())<0.03 &&
308 TMath::Abs(TMath::Abs(track1->Phi()-track2->Phi())-TMath::Pi())<0.1
309 )
d269b0e6 310 {
774986ef 311 printf("COSMIC candidate \n");
d269b0e6 312
774986ef 313 printf("track1->Pt() %f, track1->Theta() %f, track1->Eta() %f, track1->Phi() %f, track1->Charge() %d \n", track1->Pt(), track1->Theta(), track1->Eta(), track1->Phi(), track1->Charge());
314 printf("track2->Pt() %f, track2->Theta() %f, track2->Eta() %f, track2->Phi() %f, track2->Charge() %d \n", track2->Pt(), track2->Theta(), track2->Eta(), track2->Phi(), track2->Charge());
315 printf("dtheta %f, deta %f, dphi %f, dq %d \n", track1->Theta()-track2->Theta(), track1->Eta()-track2->Eta(), track1->Phi()-track2->Phi(), track1->Charge()+track2->Charge());
316 printf("dsphi %f, errsphi %f, dtanl %f, errtanl %f \n", TMath::Abs(track1->GetSnp()-track2->GetSnp()), TMath::Sqrt(track1->GetSigmaSnp2()+track2->GetSigmaSnp2()), TMath::Abs(track1->GetTgl()-track2->GetTgl()), TMath::Sqrt(track1->GetSigmaTgl2()+track2->GetSigmaTgl2()));
317 return kTRUE;
318 }
847e74b2 319
320return kFALSE;
321}
322
0aaa8b91 323//____________________________________________________________________
70fdd197 324void AlidNdPtHelper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger)
0aaa8b91 325{
326 //
327 // Prints the given configuration
328 //
329
330 TString str(">>>> Running with ");
331
332 switch (analysisMode)
333 {
334 case kInvalid: str += "invalid setting"; break;
335 case kSPD : str += "SPD-only"; break;
336 case kTPC : str += "TPC-only"; break;
337 case kTPCITS : str += "Global tracking"; break;
338 case kTPCSPDvtx : str += "TPC tracking + SPD event vertex"; break;
847e74b2 339 case kTPCSPDvtxUpdate : str += "TPC tracks updated with SPD event vertex point"; break;
791aaf54 340 case kTPCTrackSPDvtx : str += "TPC tracking + Tracks event vertex or SPD event vertex"; break;
341 case kTPCTrackSPDvtxUpdate : str += "TPC tracks updated with Track or SPD event vertex point"; break;
774986ef 342 case kTPCITSHybrid : str += "TPC tracking + ITS refit + >1 SPD cluster"; break;
791aaf54 343 case kTPCITSHybridTrackSPDvtx : str += "TPC tracking + ITS refit + >1 SPD cluster + Tracks event vertex or SPD event vertex"; break;
344 case kTPCITSHybridTrackSPDvtxDCArPt : str += "TPC tracking + ITS refit + >1 SPD cluster + Tracks event vertex or SPD event vertex + DCAr(pt)"; break;
0aaa8b91 345 case kMCRec : str += "TPC tracking + Replace rec. with MC values"; break;
0aaa8b91 346 }
347 str += " and trigger ";
348
70fdd197 349 str += AliTriggerAnalysis::GetTriggerName(trigger);
0aaa8b91 350
351 str += " <<<<";
352
353 Printf("%s", str.Data());
354}
355
356//____________________________________________________________________
357Int_t AlidNdPtHelper::ConvertPdgToPid(TParticle *particle) {
358//
359// Convert Abs(pdg) to pid
360// (0 - e, 1 - muons, 2 - pions, 3 - kaons, 4 - protons, 5 -all rest)
361//
362Int_t pid=-1;
363
364 if (TMath::Abs(particle->GetPdgCode()) == kElectron) { pid = 0; }
365 else if (TMath::Abs(particle->GetPdgCode()) == kMuonMinus) { pid = 1; }
366 else if (TMath::Abs(particle->GetPdgCode()) == kPiPlus) { pid = 2; }
367 else if (TMath::Abs(particle->GetPdgCode()) == kKPlus) { pid = 3; }
368 else if (TMath::Abs(particle->GetPdgCode()) == kProton) { pid = 4; }
369 else { pid = 5; }
370
371return pid;
372}
373
374//_____________________________________________________________________________
f2dec884 375TH1F* AlidNdPtHelper::CreateResHisto(TH2F* const hRes2, TH1F **phMean, Int_t integ, Bool_t drawBinFits, Int_t minHistEntries)
0aaa8b91 376{
bad4ba69 377//
378// Create mean and resolution
379// histograms
380//
0aaa8b91 381 TVirtualPad* currentPad = gPad;
382 TAxis* axis = hRes2->GetXaxis();
383 Int_t nBins = axis->GetNbins();
384 Bool_t overflowBinFits = kFALSE;
385 TH1F* hRes, *hMean;
386 if (axis->GetXbins()->GetSize()){
387 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
388 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
389 }
390 else{
391 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
392 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
393
394 }
395 hRes->SetStats(false);
396 hRes->SetOption("E");
397 hRes->SetMinimum(0.);
398 //
399 hMean->SetStats(false);
400 hMean->SetOption("E");
401
402 // create the fit function
403 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
404
405 fitFunc->SetLineWidth(2);
406 fitFunc->SetFillStyle(0);
407 // create canvas for fits
408 TCanvas* canBinFits = NULL;
409 Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
410 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
411 Int_t ny = (nPads-1) / nx + 1;
412 if (drawBinFits) {
413 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
414 if (canBinFits) delete canBinFits;
415 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
416 canBinFits->Divide(nx, ny);
417 }
418
419 // loop over x bins and fit projection
420 Int_t dBin = ((overflowBinFits) ? 1 : 0);
421 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
422 if (drawBinFits) canBinFits->cd(bin + dBin);
423 Int_t bin0=TMath::Max(bin-integ,0);
424 Int_t bin1=TMath::Min(bin+integ,nBins);
425 TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
426 //
427 if (hBin->GetEntries() > minHistEntries) {
428 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
429 hBin->Fit(fitFunc,"s");
430 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
431
432 if (sigma > 0.){
433 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
434 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
435 }
436 else{
437 hRes->SetBinContent(bin, 0.);
438 hMean->SetBinContent(bin,0);
439 }
440 hRes->SetBinError(bin, fitFunc->GetParError(2));
441 hMean->SetBinError(bin, fitFunc->GetParError(1));
442
443 //
444 //
445
446 } else {
447 hRes->SetBinContent(bin, 0.);
448 hRes->SetBinError(bin, 0.);
449 hMean->SetBinContent(bin, 0.);
450 hMean->SetBinError(bin, 0.);
451 }
452
453
454 if (drawBinFits) {
455 char name[256];
456 if (bin == 0) {
457 sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
458 } else if (bin == nBins+1) {
459 sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
460 } else {
461 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
462 axis->GetTitle(), axis->GetBinUpEdge(bin));
463 }
464 canBinFits->cd(bin + dBin);
465 hBin->SetTitle(name);
466 hBin->SetStats(kTRUE);
467 hBin->DrawCopy("E");
468 canBinFits->Update();
469 canBinFits->Modified();
470 canBinFits->Update();
471 }
472
473 delete hBin;
474 }
475
476 delete fitFunc;
477 currentPad->cd();
478 *phMean = hMean;
479 return hRes;
480}
481
482//_____________________________________________________________________________
483TH1F* AlidNdPtHelper::MakeResol(TH2F * his, Int_t integ, Bool_t type, Bool_t drawBins, Int_t minHistEntries){
484// Create resolution histograms
485
486 TH1F *hisr=0, *hism=0;
487 if (!gPad) new TCanvas;
0aaa8b91 488 hisr = CreateResHisto(his,&hism,integ,drawBins,minHistEntries);
489 if (type) return hism;
490 else return hisr;
491
492return hisr;
493}
494
495//_____________________________________________________________________________
bad4ba69 496TObjArray* AlidNdPtHelper::GetAllChargedTracks(AliESDEvent *esdEvent, AnalysisMode analysisMode)
0aaa8b91 497{
498 //
499 // all charged TPC particles
500 //
501 TObjArray *allTracks = new TObjArray();
502 if(!allTracks) return allTracks;
503
504 AliESDtrack *track=0;
505 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
506 {
847e74b2 507 if(analysisMode == AlidNdPtHelper::kTPC) {
774986ef 508 //
847e74b2 509 // track must be deleted by user
774986ef 510 // esd track parameters are replaced by TPCinner
511 //
bad4ba69 512 track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
847e74b2 513 if(!track) continue;
514 }
774986ef 515 else if (analysisMode == AlidNdPtHelper::kTPCSPDvtx || analysisMode == AlidNdPtHelper::kTPCSPDvtxUpdate)
847e74b2 516 {
774986ef 517 //
847e74b2 518 // track must be deleted by the user
774986ef 519 // esd track parameters are replaced by TPCinner
520 //
847e74b2 521 track = AlidNdPtHelper::GetTPCOnlyTrackSPDvtx(esdEvent,iTrack,kFALSE);
522 if(!track) continue;
523 }
791aaf54 524 else if (analysisMode == AlidNdPtHelper::kTPCTrackSPDvtx || analysisMode == AlidNdPtHelper::kTPCTrackSPDvtxUpdate)
525 {
526 //
527 // track must be deleted by the user
528 // esd track parameters are replaced by TPCinner
529 //
530 track = AlidNdPtHelper::GetTPCOnlyTrackTrackSPDvtx(esdEvent,iTrack,kFALSE);
531 if(!track) continue;
532 }
774986ef 533 else if( analysisMode == AlidNdPtHelper::kTPCITSHybrid )
534 {
535 track = AlidNdPtHelper::GetTrackSPDvtx(esdEvent,iTrack,kFALSE);
536 }
791aaf54 537 else if( analysisMode == AlidNdPtHelper::kTPCITSHybridTrackSPDvtx || analysisMode == AlidNdPtHelper::kTPCITSHybridTrackSPDvtxDCArPt)
538 {
539 track = AlidNdPtHelper::GetTrackTrackSPDvtx(esdEvent,iTrack,kFALSE);
540 }
774986ef 541 else
542 {
543 track = esdEvent->GetTrack(iTrack);
0aaa8b91 544 }
774986ef 545
0aaa8b91 546 if(!track) continue;
547
548 if(track->Charge()==0) {
791aaf54 549 if(analysisMode == AlidNdPtHelper::kTPC ||
550 analysisMode == AlidNdPtHelper::kTPCSPDvtx || analysisMode == AlidNdPtHelper::kTPCTrackSPDvtx ||
551 analysisMode == AlidNdPtHelper::kTPCSPDvtxUpdate || analysisMode == AlidNdPtHelper::kTPCTrackSPDvtxUpdate)
847e74b2 552 {
0aaa8b91 553 delete track; continue;
554 } else {
555 continue;
556 }
557 }
558
559 allTracks->Add(track);
560 }
bad4ba69 561
791aaf54 562 if(analysisMode == AlidNdPtHelper::kTPC ||
563 analysisMode == AlidNdPtHelper::kTPCSPDvtx || analysisMode == AlidNdPtHelper::kTPCTrackSPDvtx ||
564 analysisMode == AlidNdPtHelper::kTPCSPDvtxUpdate || analysisMode == AlidNdPtHelper::kTPCTrackSPDvtxUpdate) {
bad4ba69 565
566 allTracks->SetOwner(kTRUE);
567 }
0aaa8b91 568
569return allTracks;
570}
571
847e74b2 572//_____________________________________________________________________________
573AliESDtrack *AlidNdPtHelper::GetTPCOnlyTrackSPDvtx(AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
574{
575//
576// Create ESD tracks from TPCinner parameters.
577// Propagte to DCA to SPD vertex.
578// Update using SPD vertex point (parameter)
579//
580// It is user responsibility to delete these tracks
581//
582
583 if (!esdEvent) return NULL;
584 if (!esdEvent->GetPrimaryVertexSPD() ) { return NULL; }
585 if (!esdEvent->GetPrimaryVertexSPD()->GetStatus() ) { return NULL; }
586
587 //
588 AliESDtrack* track = esdEvent->GetTrack(iTrack);
589 if (!track)
590 return NULL;
591
774986ef 592 Bool_t isOK = kFALSE;
593 Double_t x[3]; track->GetXYZ(x);
594 Double_t b[3]; AliTracker::GetBxByBz(x,b);
595
847e74b2 596 // create new ESD track
597 AliESDtrack *tpcTrack = new AliESDtrack();
598
599 // relate TPC-only tracks (TPCinner) to SPD vertex
600 AliExternalTrackParam cParam;
601 if(bUpdate) {
774986ef 602 isOK = track->RelateToVertexTPCBxByBz(esdEvent->GetPrimaryVertexSPD(),b,kVeryBig,&cParam);
847e74b2 603 track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
604
605 // reject fake tracks
606 if(track->Pt() > 10000.) {
607 ::Error("Exclude no physical tracks","pt>10000. GeV");
608 delete tpcTrack;
609 return NULL;
610 }
611 }
612 else {
774986ef 613 isOK = track->RelateToVertexTPCBxByBz(esdEvent->GetPrimaryVertexSPD(), b, kVeryBig);
847e74b2 614 }
615
616 // only true if we have a tpc track
617 if (!track->FillTPCOnlyTrack(*tpcTrack))
618 {
619 delete tpcTrack;
620 return NULL;
621 }
774986ef 622
623 if(!isOK) return NULL;
847e74b2 624
625return tpcTrack;
626}
627
791aaf54 628//_____________________________________________________________________________
629AliESDtrack *AlidNdPtHelper::GetTPCOnlyTrackTrackSPDvtx(AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
630{
631//
632// Create ESD tracks from TPCinner parameters.
633// Propagte to DCA to Track or SPD vertex.
634// Update using SPD vertex point (parameter)
635//
636// It is user responsibility to delete these tracks
637//
638 if (!esdEvent) return NULL;
639 const AliESDVertex *vertex = esdEvent->GetPrimaryVertexTracks();
640 if(vertex->GetNContributors()<1) {
641 // SPD vertex
642 vertex = esdEvent->GetPrimaryVertexSPD();
643 }
644 if(!vertex) return NULL;
645
646 //
647 AliESDtrack* track = esdEvent->GetTrack(iTrack);
648 if (!track)
649 return NULL;
650
651 Bool_t isOK = kFALSE;
652 Double_t x[3]; track->GetXYZ(x);
653 Double_t b[3]; AliTracker::GetBxByBz(x,b);
654
655 // create new ESD track
656 AliESDtrack *tpcTrack = new AliESDtrack();
657
658 // relate TPC-only tracks (TPCinner) to SPD vertex
659 AliExternalTrackParam cParam;
660 if(bUpdate) {
661 isOK = track->RelateToVertexTPCBxByBz(vertex,b,kVeryBig,&cParam);
662 track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
663
664 // reject fake tracks
665 if(track->Pt() > 10000.) {
666 ::Error("Exclude no physical tracks","pt>10000. GeV");
667 delete tpcTrack;
668 return NULL;
669 }
670 }
671 else {
672 isOK = track->RelateToVertexTPCBxByBz(vertex, b, kVeryBig);
673 }
674
675 // only true if we have a tpc track
676 if (!track->FillTPCOnlyTrack(*tpcTrack))
677 {
678 delete tpcTrack;
679 return NULL;
680 }
681
682 if(!isOK) return NULL;
683
684return tpcTrack;
685}
686
774986ef 687//_____________________________________________________________________________
688AliESDtrack *AlidNdPtHelper::GetTrackSPDvtx(AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
689{
690//
691// Propagte track to DCA to SPD vertex.
692// Update using SPD vertex point (parameter)
693//
694 if (!esdEvent) return NULL;
695 if (!esdEvent->GetPrimaryVertexSPD() ) { return NULL; }
696 if (!esdEvent->GetPrimaryVertexSPD()->GetStatus() ) { return NULL; }
697
698 //
699 AliESDtrack* track = esdEvent->GetTrack(iTrack);
700 if (!track)
701 return NULL;
702
703 Bool_t isOK = kFALSE;
704 Double_t x[3]; track->GetXYZ(x);
705 Double_t b[3]; AliTracker::GetBxByBz(x,b);
706
707 // relate tracks to SPD vertex
708 AliExternalTrackParam cParam;
709 if(bUpdate) {
710 isOK = track->RelateToVertexBxByBz(esdEvent->GetPrimaryVertexSPD(),b,kVeryBig,&cParam);
711 track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
712
713 // reject fake tracks
714 if(track->Pt() > 10000.) {
715 ::Error("Exclude no physical tracks","pt>10000. GeV");
716 return NULL;
717 }
718 }
719 else {
720 isOK = track->RelateToVertexBxByBz(esdEvent->GetPrimaryVertexSPD(), b, kVeryBig);
721 }
722
723 if(!isOK) return NULL;
724
725return track;
726}
727
791aaf54 728//_____________________________________________________________________________
729AliESDtrack *AlidNdPtHelper::GetTrackTrackSPDvtx(AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
730{
731//
732// Propagte track to DCA to Track or SPD vertex.
733// Update using SPD vertex point (parameter)
734//
735 if (!esdEvent) return NULL;
736
737 const AliESDVertex *vertex = esdEvent->GetPrimaryVertexTracks();
738 if(vertex->GetNContributors()<1) {
739 // SPD vertex
740 vertex = esdEvent->GetPrimaryVertexSPD();
741 }
742 if(!vertex) return NULL;
743
744 //
745 AliESDtrack* track = esdEvent->GetTrack(iTrack);
746 if (!track)
747 return NULL;
748
749 Bool_t isOK = kFALSE;
750 Double_t x[3]; track->GetXYZ(x);
751 Double_t b[3]; AliTracker::GetBxByBz(x,b);
752
753 // relate tracks to SPD vertex
754 AliExternalTrackParam cParam;
755 if(bUpdate) {
756 isOK = track->RelateToVertexBxByBz(vertex,b,kVeryBig,&cParam);
757 track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
758
759 // reject fake tracks
760 if(track->Pt() > 10000.) {
761 ::Error("Exclude no physical tracks","pt>10000. GeV");
762 return NULL;
763 }
764 }
765 else {
766 isOK = track->RelateToVertexBxByBz(vertex, b, kVeryBig);
767 }
768
769 if(!isOK) return NULL;
770
771return track;
772}
773
0aaa8b91 774//_____________________________________________________________________________
f2dec884 775Int_t AlidNdPtHelper::GetTPCMBTrackMult(AliESDEvent *const esdEvent, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts, AliESDtrackCuts *const trackCuts)
0aaa8b91 776{
777 //
778 // get MB event track multiplicity
779 //
780 if(!esdEvent)
781 {
782 ::Error("AlidNdPtHelper::GetTPCMBTrackMult()","esd event is NULL");
783 return 0;
784 }
785
786 if(!evtCuts || !accCuts || !trackCuts)
787 {
788 ::Error("AlidNdPtHelper::GetTPCMBTrackMult()","cuts not available");
789 return 0;
790 }
791
792 //
793 Double_t pos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
794 Double_t err[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
795 AliESDVertex vtx0(pos,err);
796
797 //
798 Float_t maxDCAr = accCuts->GetMaxDCAr();
799 Float_t maxDCAz = accCuts->GetMaxDCAz();
800 Float_t minTPCClust = trackCuts->GetMinNClusterTPC();
801 //
802 Int_t ntracks = esdEvent->GetNumberOfTracks();
803 Double_t dca[2],cov[3];
804 Int_t mult=0;
805 for (Int_t i=0;i <ntracks; i++){
806 AliESDtrack *t = esdEvent->GetTrack(i);
807 if (!t) continue;
808 if (t->Charge() == 0) continue;
809 if (!t->GetTPCInnerParam()) continue;
810 if (t->GetTPCNcls()<minTPCClust) continue;
811 //
774986ef 812 Double_t x[3]; t->GetXYZ(x);
813 Double_t b[3]; AliTracker::GetBxByBz(x,b);
814 const Double_t kMaxStep = 1; //max step over the material
815
0aaa8b91 816 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
774986ef 817 if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov))
0aaa8b91 818 {
819 if(tpcTrack) delete tpcTrack;
820 continue;
821 }
822 //
823 if (TMath::Abs(dca[0])>maxDCAr || TMath::Abs(dca[1])>maxDCAz) {
824 if(tpcTrack) delete tpcTrack;
825 continue;
826 }
827
828 mult++;
829
830 if(tpcTrack) delete tpcTrack;
831 }
832
833return mult;
834}
835
836//_____________________________________________________________________________
f2dec884 837Int_t AlidNdPtHelper::GetTPCMBPrimTrackMult(AliESDEvent *const esdEvent, AliStack *const stack, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts, AliESDtrackCuts *const trackCuts)
0aaa8b91 838{
839 //
840 // get MB primary event track multiplicity
841 //
842 if(!esdEvent)
843 {
844 ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","esd event is NULL");
845 return 0;
846 }
847
848 if(!stack)
849 {
850 ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","esd event is NULL");
851 return 0;
852 }
853
854 if(!evtCuts || !accCuts || !trackCuts)
855 {
856 ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","cuts not available");
857 return 0;
858 }
859
860 //
861 Double_t pos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
862 Double_t err[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
863 AliESDVertex vtx0(pos,err);
864
865 //
866 Float_t maxDCAr = accCuts->GetMaxDCAr();
867 Float_t maxDCAz = accCuts->GetMaxDCAz();
868 Float_t minTPCClust = trackCuts->GetMinNClusterTPC();
869
870 //
871 Int_t ntracks = esdEvent->GetNumberOfTracks();
872 Double_t dca[2],cov[3];
873 Int_t mult=0;
874 for (Int_t i=0;i <ntracks; i++){
875 AliESDtrack *t = esdEvent->GetTrack(i);
876 if (!t) continue;
877 if (t->Charge() == 0) continue;
878 if (!t->GetTPCInnerParam()) continue;
879 if (t->GetTPCNcls()<minTPCClust) continue;
880 //
774986ef 881 Double_t x[3]; t->GetXYZ(x);
882 Double_t b[3]; AliTracker::GetBxByBz(x,b);
883 const Double_t kMaxStep = 1; //max step over the material
884
0aaa8b91 885 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
774986ef 886 if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov))
0aaa8b91 887 {
888 if(tpcTrack) delete tpcTrack;
889 continue;
890 }
891 //
892 if (TMath::Abs(dca[0])>maxDCAr || TMath::Abs(dca[1])>maxDCAz) {
893 if(tpcTrack) delete tpcTrack;
894 continue;
895 }
896
897 Int_t label = TMath::Abs(t->GetLabel());
898 TParticle *part = stack->Particle(label);
899 if(!part) {
900 if(tpcTrack) delete tpcTrack;
901 continue;
902 }
903 if(!stack->IsPhysicalPrimary(label))
904 {
905 if(tpcTrack) delete tpcTrack;
906 continue;
907 }
908
909 mult++;
910
911 if(tpcTrack) delete tpcTrack;
912 }
913
914return mult;
915}
916
917
791aaf54 918//_____________________________________________________________________________
919Double_t AlidNdPtHelper::GetStrangenessCorrFactor(const Double_t pt)
920{
921// data driven correction factor for secondaries
922// underestimated secondaries with strangeness in Pythia (A. Dainese)
0aaa8b91 923
791aaf54 924//
925// pt=0.17; fact=1
926// pt=0.4; fact=1.07
927// pt=0.6; fact=1.25
928// pt>=1.2; fact=1.5
929//
0aaa8b91 930
791aaf54 931if (pt <= 0.17) return 1.0;
932if (pt <= 0.4) return GetLinearInterpolationValue(0.17,1.0,0.4,1.07, pt);
933if (pt <= 0.6) return GetLinearInterpolationValue(0.4,1.07,0.6,1.25, pt);
934if (pt <= 1.2) return GetLinearInterpolationValue(0.6,1.25,1.2,1.5, pt);
935return 1.5;
936
937}
938
939//___________________________________________________________________________
940Double_t AlidNdPtHelper::GetLinearInterpolationValue(Double_t const x1, Double_t const y1, Double_t const x2, Double_t const y2, const Double_t pt)
941{
942//
943// linear interpolation
944//
945 return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
946}
0aaa8b91 947
948//_____________________________________________________________________________
f2dec884 949Int_t AlidNdPtHelper::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts)
0aaa8b91 950{
951 //
952 // calculate mc event true track multiplicity
953 //
954 if(!mcEvent) return 0;
955
956 AliStack* stack = 0;
957 Int_t mult = 0;
958
959 // MC particle stack
960 stack = mcEvent->Stack();
961 if (!stack) return 0;
962
963 Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
964 if(!isEventOK) return 0;
965
966 Int_t nPart = stack->GetNtrack();
967 for (Int_t iMc = 0; iMc < nPart; ++iMc)
968 {
969 TParticle* particle = stack->Particle(iMc);
970 if (!particle)
971 continue;
972
973 // only charged particles
974 Double_t charge = particle->GetPDG()->Charge()/3.;
f2dec884 975 if (TMath::Abs(charge) < 0.001)
0aaa8b91 976 continue;
977
978 // physical primary
979 Bool_t prim = stack->IsPhysicalPrimary(iMc);
980 if(!prim) continue;
981
774986ef 982 // checked accepted without pt cut
983 //if(accCuts->AcceptTrack(particle))
984 if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
0aaa8b91 985 {
986 mult++;
987 }
988 }
989
990return mult;
991}
992
993//_______________________________________________________________________
f2dec884 994void AlidNdPtHelper::PrintMCInfo(AliStack *const pStack,Int_t label)
0aaa8b91 995{
996// print information about particles in the stack
997
998 if(!pStack)return;
999 label = TMath::Abs(label);
1000 TParticle *part = pStack->Particle(label);
1001 Printf("########################");
1002 Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
1003 part->Print();
1004 TParticle* mother = part;
1005 Int_t imo = part->GetFirstMother();
1006 Int_t nprim = pStack->GetNprimary();
1007
1008 while((imo >= nprim)) {
1009 mother = pStack->Particle(imo);
1010 Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
1011 mother->Print();
1012 imo = mother->GetFirstMother();
1013 }
1014
1015 Printf("########################");
1016}
1017
1018
1019//_____________________________________________________________________________
f2dec884 1020TH1* AlidNdPtHelper::GetContCorrHisto(TH1 *const hist)
0aaa8b91 1021{
1022//
1023// get contamination histogram
1024//
1025 if(!hist) return 0;
1026
1027 Int_t nbins = hist->GetNbinsX();
f2dec884 1028 TH1 *hCont = (TH1D *)hist->Clone();
0aaa8b91 1029
1030 for(Int_t i=0; i<=nbins+1; i++) {
1031 Double_t binContent = hist->GetBinContent(i);
1032 Double_t binError = hist->GetBinError(i);
1033
f2dec884 1034 hCont->SetBinContent(i,1.-binContent);
1035 hCont->SetBinError(i,binError);
0aaa8b91 1036 }
1037
f2dec884 1038return hCont;
0aaa8b91 1039}
1040
1041
1042//_____________________________________________________________________________
f2dec884 1043TH1* AlidNdPtHelper::ScaleByBinWidth(TH1 *const hist)
0aaa8b91 1044{
1045//
1046// scale by bin width
1047//
1048 if(!hist) return 0;
1049
f2dec884 1050 TH1 *hScale = (TH1D *)hist->Clone();
1051 hScale->Scale(1.,"width");
0aaa8b91 1052
f2dec884 1053return hScale;
0aaa8b91 1054}
1055
1056//_____________________________________________________________________________
f2dec884 1057TH1* AlidNdPtHelper::CalcRelativeDifference(TH1 *const hist1, TH1 *const hist2)
0aaa8b91 1058{
1059//
1060// calculate rel. difference
1061//
1062
1063 if(!hist1) return 0;
1064 if(!hist2) return 0;
1065
f2dec884 1066 TH1 *h1Clone = (TH1D *)hist1->Clone();
1067 h1Clone->Sumw2();
0aaa8b91 1068
1069 // (rec-mc)/mc
f2dec884 1070 h1Clone->Add(hist2,-1);
1071 h1Clone->Divide(hist2);
0aaa8b91 1072
f2dec884 1073return h1Clone;
0aaa8b91 1074}
1075
1076//_____________________________________________________________________________
f2dec884 1077TH1* AlidNdPtHelper::CalcRelativeDifferenceFun(TH1 *const hist1, TF1 *const fun)
0aaa8b91 1078{
1079//
08a80bfe 1080// calculate rel. difference
0aaa8b91 1081// between histogram and function
1082//
1083 if(!hist1) return 0;
1084 if(!fun) return 0;
1085
f2dec884 1086 TH1 *h1Clone = (TH1D *)hist1->Clone();
1087 h1Clone->Sumw2();
0aaa8b91 1088
1089 //
f2dec884 1090 h1Clone->Add(fun,-1);
1091 h1Clone->Divide(hist1);
0aaa8b91 1092
f2dec884 1093return h1Clone;
0aaa8b91 1094}
1095
1096//_____________________________________________________________________________
f2dec884 1097TH1* AlidNdPtHelper::NormalizeToEvent(TH2 *const hist1, TH1 *const hist2)
0aaa8b91 1098{
1099// normalise to event for a given multiplicity bin
1100// return pt histogram
1101
1102 if(!hist1) return 0;
1103 if(!hist2) return 0;
1104 char name[256];
1105
1106 Int_t nbinsX = hist1->GetNbinsX();
1107 //Int_t nbinsY = hist1->GetNbinsY();
1108
f2dec884 1109 TH1D *histNorm = 0;
0aaa8b91 1110 for(Int_t i=0; i<=nbinsX+1; i++) {
1111 sprintf(name,"mom_%d",i);
1112 TH1D *hist = (TH1D*)hist1->ProjectionY(name,i+1,i+1);
1113
1114 sprintf(name,"mom_norm");
1115 if(i==0) {
f2dec884 1116 histNorm = (TH1D *)hist->Clone(name);
1117 histNorm->Reset();
0aaa8b91 1118 }
1119
1120 Double_t nbEvents = hist2->GetBinContent(i);
1121 if(!nbEvents) { nbEvents = 1.; };
1122
1123 hist->Scale(1./nbEvents);
f2dec884 1124 histNorm->Add(hist);
0aaa8b91 1125 }
1126
f2dec884 1127return histNorm;
0aaa8b91 1128}
1129
1130//_____________________________________________________________________________
f2dec884 1131THnSparse* AlidNdPtHelper::GenerateCorrMatrix(THnSparse *const hist1, THnSparse *const hist2, char *const name) {
0aaa8b91 1132// generate correction matrix
1133if(!hist1 || !hist2) return 0;
1134
1135THnSparse *h =(THnSparse*)hist1->Clone(name);;
1136h->Divide(hist1,hist2,1,1,"B");
1137
1138return h;
1139}
1140
1141//_____________________________________________________________________________
f2dec884 1142TH2* AlidNdPtHelper::GenerateCorrMatrix(TH2 *const hist1, TH2 *const hist2, char *const name) {
0aaa8b91 1143// generate correction matrix
1144if(!hist1 || !hist2) return 0;
1145
1146TH2D *h =(TH2D*)hist1->Clone(name);;
1147h->Divide(hist1,hist2,1,1,"B");
1148
1149return h;
1150}
1151
1152//_____________________________________________________________________________
f2dec884 1153TH1* AlidNdPtHelper::GenerateCorrMatrix(TH1 *const hist1, TH1 *const hist2, char *const name) {
0aaa8b91 1154// generate correction matrix
1155if(!hist1 || !hist2) return 0;
1156
1157TH1D *h =(TH1D*)hist1->Clone(name);;
1158h->Divide(hist1,hist2,1,1,"B");
1159
1160return h;
1161}
1162
1163//_____________________________________________________________________________
f2dec884 1164THnSparse* AlidNdPtHelper::GenerateContCorrMatrix(THnSparse *const hist1, THnSparse *const hist2, char *const name) {
0aaa8b91 1165// generate contamination correction matrix
1166if(!hist1 || !hist2) return 0;
1167
1168THnSparse *hist = GenerateCorrMatrix(hist1, hist2, name);
1169if(!hist) return 0;
1170
1171// only for non ZERO bins!!!!
1172
1173Int_t* coord = new Int_t[hist->GetNdimensions()];
1174memset(coord, 0, sizeof(Int_t) * hist->GetNdimensions());
1175
1176 for (Long64_t i = 0; i < hist->GetNbins(); ++i) {
1177 Double_t v = hist->GetBinContent(i, coord);
1178 hist->SetBinContent(coord, 1.0-v);
1179 //printf("v %f, hist->GetBinContent(i, coord) %f \n",v,hist->GetBinContent(i, coord));
1180 Double_t err = hist->GetBinError(coord);
1181 hist->SetBinError(coord, err);
1182 }
1183
1184delete [] coord;
1185
1186return hist;
1187}
1188
1189//_____________________________________________________________________________
f2dec884 1190TH2* AlidNdPtHelper::GenerateContCorrMatrix(TH2 *const hist1, TH2 *const hist2, char *const name) {
0aaa8b91 1191// generate contamination correction matrix
1192if(!hist1 || !hist2) return 0;
1193
1194TH2 *hist = GenerateCorrMatrix(hist1, hist2, name);
1195if(!hist) return 0;
1196
1197Int_t nBinsX = hist->GetNbinsX();
1198Int_t nBinsY = hist->GetNbinsY();
1199
1200 for (Int_t i = 0; i < nBinsX+1; i++) {
1201 for (Int_t j = 0; j < nBinsY+1; j++) {
1202 Double_t cont = hist->GetBinContent(i,j);
1203 hist->SetBinContent(i,j,1.-cont);
1204 Double_t err = hist->GetBinError(i,j);
1205 hist->SetBinError(i,j,err);
1206 }
1207 }
1208
1209return hist;
1210}
1211
1212//_____________________________________________________________________________
f2dec884 1213TH1* AlidNdPtHelper::GenerateContCorrMatrix(TH1 *const hist1, TH1 *const hist2, char *const name) {
0aaa8b91 1214// generate contamination correction matrix
1215if(!hist1 || !hist2) return 0;
1216
1217TH1 *hist = GenerateCorrMatrix(hist1, hist2, name);
1218if(!hist) return 0;
1219
1220Int_t nBinsX = hist->GetNbinsX();
1221
1222 for (Int_t i = 0; i < nBinsX+1; i++) {
1223 Double_t cont = hist->GetBinContent(i);
1224 hist->SetBinContent(i,1.-cont);
1225 Double_t err = hist->GetBinError(i);
1226 hist->SetBinError(i,err);
1227 }
1228
1229return hist;
1230}
1231
1232//_____________________________________________________________________________
f2dec884 1233const AliESDVertex* AlidNdPtHelper::GetTPCVertexZ(AliESDEvent* const esdEvent, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts, AliESDtrackCuts *const trackCuts, Float_t fraction, Int_t ntracksMin){
0aaa8b91 1234 //
1235 // TPC Z vertexer
1236 //
bad4ba69 1237 if(!esdEvent)
1238 {
1239 ::Error("AlidNdPtHelper::GetTPCVertexZ()","cuts not available");
1240 return NULL;
1241 }
1242
1243 if(!evtCuts || !accCuts || !trackCuts)
1244 {
1245 ::Error("AlidNdPtHelper::GetTPCVertexZ()","cuts not available");
1246 return NULL;
1247 }
1248
1249 Double_t vtxpos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
1250 Double_t vtxsigma[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
0aaa8b91 1251 AliESDVertex vtx0(vtxpos,vtxsigma);
bad4ba69 1252
1253 Double_t maxDCAr = accCuts->GetMaxDCAr();
1254 Double_t maxDCAz = accCuts->GetMaxDCAz();
1255 Int_t minTPCClust = trackCuts->GetMinNClusterTPC();
1256
0aaa8b91 1257 //
1258 Int_t ntracks = esdEvent->GetNumberOfTracks();
1259 TVectorD ztrack(ntracks);
0aaa8b91 1260 Double_t dca[2],cov[3];
1261 Int_t counter=0;
1262 for (Int_t i=0;i <ntracks; i++){
1263 AliESDtrack *t = esdEvent->GetTrack(i);
1264 if (!t) continue;
1265 if (!t->GetTPCInnerParam()) continue;
bad4ba69 1266 if (t->GetTPCNcls()<minTPCClust) continue;
0aaa8b91 1267 //
774986ef 1268
1269 Double_t x[3]; t->GetXYZ(x);
1270 Double_t b[3]; AliTracker::GetBxByBz(x,b);
1271 const Double_t kMaxStep = 1; //max step over the material
1272
0aaa8b91 1273 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
774986ef 1274 if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov)) continue;
0aaa8b91 1275
1276 //
bad4ba69 1277 if (TMath::Abs(dca[0])>maxDCAr) continue;
1278 //if (TMath::Sqrt(cov[0])>sigmaXYcut) continue;
1279 if (TMath::Abs(tpcTrack->GetZ())>maxDCAz) continue;
0aaa8b91 1280
0aaa8b91 1281 ztrack[counter]=tpcTrack->GetZ();
1282 counter++;
1283
1284 if(tpcTrack) delete tpcTrack;
1285 }
bad4ba69 1286
0aaa8b91 1287 //
1288 // Find LTM z position
1289 //
1290 Double_t mean=0, sigma=0;
1291 if (counter<ntracksMin) return 0;
1292 //
1293 Int_t nused = TMath::Nint(counter*fraction);
1294 if (nused==counter) nused=counter-1;
1295 if (nused>1){
1296 AliMathBase::EvaluateUni(counter, ztrack.GetMatrixArray(), mean,sigma, TMath::Nint(counter*fraction));
1297 sigma/=TMath::Sqrt(nused);
1298 }else{
1299 mean = TMath::Mean(counter, ztrack.GetMatrixArray());
1300 sigma = TMath::RMS(counter, ztrack.GetMatrixArray());
1301 sigma/=TMath::Sqrt(counter-1);
1302 }
1303 vtxpos[2]=mean;
1304 vtxsigma[2]=sigma;
1305 const AliESDVertex* vertex = new AliESDVertex(vtxpos, vtxsigma);
1306 return vertex;
1307}
1308
1309//_____________________________________________________________________________
f2dec884 1310Int_t AlidNdPtHelper::GetSPDMBTrackMult(AliESDEvent* const esdEvent, Float_t deltaThetaCut, Float_t deltaPhiCut)
0aaa8b91 1311{
1312 //
1313 // SPD track multiplicity
1314 //
1315
1316 // get tracklets
1317 const AliMultiplicity* mult = esdEvent->GetMultiplicity();
1318 if (!mult)
1319 return 0;
1320
1321 // get multiplicity from SPD tracklets
1322 Int_t inputCount = 0;
1323 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
1324 {
1325 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
1326
1327 Float_t phi = mult->GetPhi(i);
1328 if (phi < 0)
1329 phi += TMath::Pi() * 2;
1330 Float_t deltaPhi = mult->GetDeltaPhi(i);
1331 Float_t deltaTheta = mult->GetDeltaTheta(i);
1332
1333 if (TMath::Abs(deltaPhi) > 1)
1334 printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
1335
1336 if (deltaThetaCut > 0. && TMath::Abs(deltaTheta) > deltaThetaCut)
1337 continue;
1338
1339 if (deltaPhiCut > 0. && TMath::Abs(deltaPhi) > deltaPhiCut)
1340 continue;
1341
1342 ++inputCount;
1343 }
1344
1345return inputCount;
1346}
1347
1348//_____________________________________________________________________________
f2dec884 1349Int_t AlidNdPtHelper::GetSPDMBPrimTrackMult(AliESDEvent* const esdEvent, AliStack* const stack, Float_t deltaThetaCut, Float_t deltaPhiCut)
0aaa8b91 1350{
1351 //
1352 // SPD track multiplicity
1353 //
1354
1355 // get tracklets
1356 const AliMultiplicity* mult = esdEvent->GetMultiplicity();
1357 if (!mult)
1358 return 0;
1359
1360 // get multiplicity from SPD tracklets
1361 Int_t inputCount = 0;
1362 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
1363 {
1364 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
1365
1366 Float_t phi = mult->GetPhi(i);
1367 if (phi < 0)
1368 phi += TMath::Pi() * 2;
1369 Float_t deltaPhi = mult->GetDeltaPhi(i);
1370 Float_t deltaTheta = mult->GetDeltaTheta(i);
1371
1372 if (TMath::Abs(deltaPhi) > 1)
1373 printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
1374
1375 if (deltaThetaCut > 0. && TMath::Abs(deltaTheta) > deltaThetaCut)
1376 continue;
1377
1378 if (deltaPhiCut > 0. && TMath::Abs(deltaPhi) > deltaPhiCut)
1379 continue;
1380
1381
1382 if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) ||
1383 !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
1384 continue;
1385
1386
1387 ++inputCount;
1388 }
1389
1390return inputCount;
1391}
1392
1393