]> 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
Warning fixed
[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//____________________________________________________________________
e9f7cb15 52const AliESDVertex* AlidNdPtHelper::GetVertex(AliESDEvent* const aEsd, const AlidNdPtEventCuts *const evtCuts, const AlidNdPtAcceptanceCuts *const accCuts, const 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;
24c88fc4 73 if (analysisMode == kSPD ||
774986ef 74 analysisMode == kTPCSPDvtx || analysisMode == kTPCSPDvtxUpdate || analysisMode == kTPCITSHybrid)
0aaa8b91 75 {
76 vertex = aEsd->GetPrimaryVertexSPD();
77 if (debug)
78 Printf("AlidNdPtHelper::GetVertex: Returning SPD vertex");
24c88fc4 79 }
80 else if (analysisMode == kTPCITS || analysisMode == kTPCTrackSPDvtx || analysisMode == kTPCTrackSPDvtxUpdate ||
695facdf 81 analysisMode == kTPCITSHybridTrackSPDvtx || analysisMode == kTPCITSHybridTrackSPDvtxDCArPt ||
82 analysisMode == kITSStandAloneTrackSPDvtx || analysisMode == kITSStandAloneTPCTrackSPDvtx)
791aaf54 83 {
84 vertex = aEsd->GetPrimaryVertexTracks();
85 if(!vertex) return NULL;
86 if(vertex->GetNContributors()<1) {
87 // SPD vertex
88 vertex = aEsd->GetPrimaryVertexSPD();
89 }
90 }
847e74b2 91 else if (analysisMode == kTPC)
0aaa8b91 92 {
93 if(bRedoTPC) {
94
95 Double_t kBz = aEsd->GetMagneticField();
96 AliVertexerTracks vertexer(kBz);
97
98 if(bUseMeanVertex) {
99 Double_t pos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
100 Double_t err[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
0aaa8b91 101 initVertex = new AliESDVertex(pos,err);
102 vertexer.SetVtxStart(initVertex);
103 vertexer.SetConstraintOn();
104 }
105
0aaa8b91 106 Double_t maxDCAr = accCuts->GetMaxDCAr();
107 Double_t maxDCAz = accCuts->GetMaxDCAz();
108 Int_t minTPCClust = trackCuts->GetMinNClusterTPC();
109
bad4ba69 110 //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 111 vertexer.SetTPCMode(0.1,1.0,5.0,minTPCClust,1,3.,0.1,2.0,maxDCAr,maxDCAz,1,4);
112
113 // TPC track preselection
114 Int_t ntracks = aEsd->GetNumberOfTracks();
115 TObjArray array(ntracks);
116 UShort_t *id = new UShort_t[ntracks];
117
118 Int_t count=0;
119 for (Int_t i=0;i <ntracks; i++) {
120 AliESDtrack *t = aEsd->GetTrack(i);
121 if (!t) continue;
122 if (t->Charge() == 0) continue;
123 if (!t->GetTPCInnerParam()) continue;
124 if (t->GetTPCNcls()<vertexer.GetMinClusters()) continue;
125 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
126 if(tpcTrack) {
127 array.AddLast(tpcTrack);
128 id[count] = (UShort_t)t->GetID();
129 count++;
130 }
131 }
132 AliESDVertex *vTPC = vertexer.VertexForSelectedTracks(&array,id, kTRUE, kTRUE, bUseMeanVertex);
bad4ba69 133
134 // set recreated TPC vertex
0aaa8b91 135 aEsd->SetPrimaryVertexTPC(vTPC);
136
137 for (Int_t i=0; i<aEsd->GetNumberOfTracks(); i++) {
138 AliESDtrack *t = aEsd->GetTrack(i);
774986ef 139 if(!t) continue;
140
141 Double_t x[3]; t->GetXYZ(x);
142 Double_t b[3]; AliTracker::GetBxByBz(x,b);
143 t->RelateToVertexTPCBxByBz(vTPC, b, kVeryBig);
0aaa8b91 144 }
145
146 delete vTPC;
147 array.Delete();
148 delete [] id; id=NULL;
149
150 }
151 vertex = aEsd->GetPrimaryVertexTPC();
152 if (debug)
774986ef 153 Printf("AlidNdPtHelper::GetVertex: Returning vertex from tracks");
154 }
155 else
156 Printf("AlidNdPtHelper::GetVertex: ERROR: Invalid second argument %d", analysisMode);
0aaa8b91 157
158 if (!vertex) {
159 if (debug)
160 Printf("AlidNdPtHelper::GetVertex: No vertex found in ESD");
161 return 0;
162 }
163
164 if (debug)
165 {
166 Printf("AlidNdPtHelper::GetVertex: Returning valid vertex: %s", vertex->GetTitle());
167 vertex->Print();
168 }
169
170 if(initVertex) delete initVertex; initVertex=NULL;
171 return vertex;
172}
173
bad4ba69 174//____________________________________________________________________
791aaf54 175Bool_t AlidNdPtHelper::TestRecVertex(const AliESDVertex* vertex, const AliESDVertex* vertexSPD, AnalysisMode analysisMode, Bool_t debug)
bad4ba69 176{
177 // Checks if a vertex meets the needed quality criteria
178 if(!vertex) return kFALSE;
774986ef 179 if(!vertex->GetStatus()) return kFALSE;
bad4ba69 180
181 Float_t requiredZResolution = -1;
774986ef 182 if (analysisMode == kSPD || analysisMode == kTPCITS ||
791aaf54 183 analysisMode == kTPCSPDvtx || analysisMode == kTPCSPDvtxUpdate || analysisMode == kTPCITSHybrid ||
184 analysisMode == kTPCTrackSPDvtx || analysisMode == kTPCTrackSPDvtxUpdate || analysisMode == kTPCITSHybridTrackSPDvtx || analysisMode == kTPCITSHybridTrackSPDvtxDCArPt
695facdf 185 || analysisMode == kITSStandAloneTrackSPDvtx || analysisMode == kITSStandAloneTPCTrackSPDvtx)
bad4ba69 186 {
774986ef 187 requiredZResolution = 1000;
bad4ba69 188 }
847e74b2 189 else if (analysisMode == kTPC)
bad4ba69 190 requiredZResolution = 10.;
191
774986ef 192 // check resolution
193 Double_t zRes = vertex->GetZRes();
194
195 if (zRes > requiredZResolution) {
196 if (debug)
197 Printf("AlidNdPtHelper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
198 return kFALSE;
199 }
791aaf54 200
201 // always check for SPD vertex
202 if(!vertexSPD) return kFALSE;
203 if(!vertexSPD->GetStatus()) return kFALSE;
204 if (vertexSPD->IsFromVertexerZ())
774986ef 205 {
791aaf54 206 if (vertexSPD->GetDispersion() > 0.02)
774986ef 207 {
208 if (debug)
209 Printf("AliPWG0Helper::TestVertex: Delta Phi too large in Vertexer Z: %f (required: %f", vertex->GetDispersion(), 0.02);
210 return kFALSE;
211 }
212 }
213
bad4ba69 214
215 return kTRUE;
216}
217
791aaf54 218//____________________________________________________________________
e9f7cb15 219Bool_t AlidNdPtHelper::IsGoodImpPar(const AliESDtrack *const track)
791aaf54 220{
221//
222// check whether particle has good DCAr(Pt) impact
223// parameter. Only for TPC+ITS tracks (7*sigma cut)
224// Origin: Andrea Dainese
225//
226
227Float_t d0z0[2],covd0z0[3];
228track->GetImpactParameters(d0z0,covd0z0);
229Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
230Float_t d0max = 7.*sigma;
231if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
232
233return kFALSE;
234}
235
0aaa8b91 236//____________________________________________________________________
f2dec884 237Bool_t AlidNdPtHelper::IsPrimaryParticle(AliStack* const stack, Int_t idx, ParticleMode particleMode)
0aaa8b91 238{
774986ef 239//
0aaa8b91 240// check primary particles
847e74b2 241// depending on the particle mode
0aaa8b91 242//
243 if(!stack) return kFALSE;
244
245 TParticle* particle = stack->Particle(idx);
246 if (!particle) return kFALSE;
247
248 // only charged particles
249 Double_t charge = particle->GetPDG()->Charge()/3.;
f2dec884 250 if (TMath::Abs(charge) < 0.001) return kFALSE;
0aaa8b91 251
252 Int_t pdg = TMath::Abs(particle->GetPdgCode());
253
254 // physical primary
255 Bool_t prim = stack->IsPhysicalPrimary(idx);
256
847e74b2 257 if(particleMode==kMCPion) {
0aaa8b91 258 if(prim && pdg==kPiPlus) return kTRUE;
259 else return kFALSE;
260 }
261
847e74b2 262 if (particleMode==kMCKaon) {
0aaa8b91 263 if(prim && pdg==kKPlus) return kTRUE;
264 else return kFALSE;
265 }
266
847e74b2 267 if (particleMode==kMCProton) {
0aaa8b91 268 if(prim && pdg==kProton) return kTRUE;
269 else return kFALSE;
270 }
271
17e8c701 272 if(particleMode==kMCRest) {
273 if(prim && pdg!=kPiPlus && pdg!=kKPlus && pdg!=kProton) return kTRUE;
274 else return kFALSE;
275 }
276
0aaa8b91 277return prim;
278}
279
d269b0e6 280//____________________________________________________________________
f2dec884 281Bool_t AlidNdPtHelper::IsCosmicTrack(AliESDtrack *const track1, AliESDtrack *const track2)
d269b0e6 282{
283//
284// check cosmic tracks
285//
286 if(!track1) return kFALSE;
287 if(!track2) return kFALSE;
d269b0e6 288
774986ef 289 //
290 // cosmic tracks in TPC
291 //
292 //if( TMath::Abs( track1->Theta() - track2->Theta() ) < 0.004 &&
293 // ((TMath::Abs(track1->Phi()-track2->Phi()) - TMath::Pi() )<0.004) &&
294 //if( track1->Pt() > 4.0 && (TMath::Abs(track1->Phi()-track2->Phi())-TMath::Pi())<0.1
295 // && (TMath::Abs(track1->Eta()+track2->Eta())-0.1) < 0.0 && (track1->Charge()+track2->Charge()) == 0)
296
297 //Float_t scaleF= 6.0;
298 if ( track1->Pt() > 4 && track2->Pt() > 4 &&
299 //(TMath::Abs(track1->GetSnp()-track2->GetSnp())-2.) < scaleF * TMath::Sqrt(track1->GetSigmaSnp2()+track2->GetSigmaSnp2()) &&
300 //TMath::Abs(track1->GetTgl()-track2->GetTgl()) < scaleF * TMath::Sqrt(track1->GetSigmaTgl2()+track2->GetSigmaTgl2()) &&
301 //TMath::Abs(track1->OneOverPt()-track2->OneOverPt()) < scaleF * TMath::Sqrt(track1->GetSigma1Pt2()+track2->GetSigma1Pt2()) &&
302 (track1->Charge()+track2->Charge()) == 0 &&
303 track1->Eta()*track2->Eta()<0.0 && TMath::Abs(track1->Eta()+track2->Eta())<0.03 &&
304 TMath::Abs(TMath::Abs(track1->Phi()-track2->Phi())-TMath::Pi())<0.1
305 )
d269b0e6 306 {
774986ef 307 printf("COSMIC candidate \n");
d269b0e6 308
774986ef 309 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());
310 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());
311 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());
312 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()));
313 return kTRUE;
314 }
847e74b2 315
316return kFALSE;
317}
318
0aaa8b91 319//____________________________________________________________________
70fdd197 320void AlidNdPtHelper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger)
0aaa8b91 321{
322 //
323 // Prints the given configuration
324 //
325
326 TString str(">>>> Running with ");
327
328 switch (analysisMode)
329 {
330 case kInvalid: str += "invalid setting"; break;
331 case kSPD : str += "SPD-only"; break;
332 case kTPC : str += "TPC-only"; break;
333 case kTPCITS : str += "Global tracking"; break;
334 case kTPCSPDvtx : str += "TPC tracking + SPD event vertex"; break;
847e74b2 335 case kTPCSPDvtxUpdate : str += "TPC tracks updated with SPD event vertex point"; break;
791aaf54 336 case kTPCTrackSPDvtx : str += "TPC tracking + Tracks event vertex or SPD event vertex"; break;
337 case kTPCTrackSPDvtxUpdate : str += "TPC tracks updated with Track or SPD event vertex point"; break;
774986ef 338 case kTPCITSHybrid : str += "TPC tracking + ITS refit + >1 SPD cluster"; break;
791aaf54 339 case kTPCITSHybridTrackSPDvtx : str += "TPC tracking + ITS refit + >1 SPD cluster + Tracks event vertex or SPD event vertex"; break;
340 case kTPCITSHybridTrackSPDvtxDCArPt : str += "TPC tracking + ITS refit + >1 SPD cluster + Tracks event vertex or SPD event vertex + DCAr(pt)"; break;
695facdf 341 case kITSStandAloneTrackSPDvtx : str += "ITS standalone + Tracks event vertex or SPD event vertex + DCAr(pt)"; break;
342 case kITSStandAloneTPCTrackSPDvtx : str += "kITSStandAloneTPCTrackSPDvtx + TPC stand alone track + Tracks event vertex or SPD event vertex + DCAr(pt)"; break;
0aaa8b91 343 case kMCRec : str += "TPC tracking + Replace rec. with MC values"; break;
0aaa8b91 344 }
345 str += " and trigger ";
346
70fdd197 347 str += AliTriggerAnalysis::GetTriggerName(trigger);
0aaa8b91 348
349 str += " <<<<";
350
351 Printf("%s", str.Data());
352}
353
354//____________________________________________________________________
e9f7cb15 355Int_t AlidNdPtHelper::ConvertPdgToPid(const TParticle *const particle) {
0aaa8b91 356//
357// Convert Abs(pdg) to pid
358// (0 - e, 1 - muons, 2 - pions, 3 - kaons, 4 - protons, 5 -all rest)
359//
360Int_t pid=-1;
361
362 if (TMath::Abs(particle->GetPdgCode()) == kElectron) { pid = 0; }
363 else if (TMath::Abs(particle->GetPdgCode()) == kMuonMinus) { pid = 1; }
364 else if (TMath::Abs(particle->GetPdgCode()) == kPiPlus) { pid = 2; }
365 else if (TMath::Abs(particle->GetPdgCode()) == kKPlus) { pid = 3; }
366 else if (TMath::Abs(particle->GetPdgCode()) == kProton) { pid = 4; }
367 else { pid = 5; }
368
369return pid;
370}
371
372//_____________________________________________________________________________
f2dec884 373TH1F* AlidNdPtHelper::CreateResHisto(TH2F* const hRes2, TH1F **phMean, Int_t integ, Bool_t drawBinFits, Int_t minHistEntries)
0aaa8b91 374{
bad4ba69 375//
376// Create mean and resolution
377// histograms
378//
0aaa8b91 379 TVirtualPad* currentPad = gPad;
380 TAxis* axis = hRes2->GetXaxis();
381 Int_t nBins = axis->GetNbins();
62e3b4b6 382 //Bool_t overflowBinFits = kFALSE;
0aaa8b91 383 TH1F* hRes, *hMean;
384 if (axis->GetXbins()->GetSize()){
385 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
386 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
387 }
388 else{
389 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
390 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
391
392 }
393 hRes->SetStats(false);
394 hRes->SetOption("E");
395 hRes->SetMinimum(0.);
396 //
397 hMean->SetStats(false);
398 hMean->SetOption("E");
399
400 // create the fit function
401 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
402
403 fitFunc->SetLineWidth(2);
404 fitFunc->SetFillStyle(0);
405 // create canvas for fits
406 TCanvas* canBinFits = NULL;
62e3b4b6 407 //Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
408 Int_t nPads = nBins;
0aaa8b91 409 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
410 Int_t ny = (nPads-1) / nx + 1;
411 if (drawBinFits) {
412 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
413 if (canBinFits) delete canBinFits;
414 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
415 canBinFits->Divide(nx, ny);
416 }
417
418 // loop over x bins and fit projection
62e3b4b6 419 //Int_t dBin = ((overflowBinFits) ? 1 : 0);
420 Int_t dBin = 0;
0aaa8b91 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 }
695facdf 533 else if(analysisMode == AlidNdPtHelper::kTPCITSHybrid )
774986ef 534 {
535 track = AlidNdPtHelper::GetTrackSPDvtx(esdEvent,iTrack,kFALSE);
536 }
68f10917 537 else if(analysisMode == AlidNdPtHelper::kTPCITSHybridTrackSPDvtx || analysisMode == AlidNdPtHelper::kTPCITSHybridTrackSPDvtxDCArPt || analysisMode == AlidNdPtHelper::kITSStandAloneTrackSPDvtx || analysisMode ==AlidNdPtHelper::kITSStandAloneTPCTrackSPDvtx)
791aaf54 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//_____________________________________________________________________________
e9f7cb15 573AliESDtrack *AlidNdPtHelper::GetTPCOnlyTrackSPDvtx(const AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
847e74b2 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//_____________________________________________________________________________
e9f7cb15 629AliESDtrack *AlidNdPtHelper::GetTPCOnlyTrackTrackSPDvtx(const AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
791aaf54 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//_____________________________________________________________________________
e9f7cb15 688AliESDtrack *AlidNdPtHelper::GetTrackSPDvtx(const AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
774986ef 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//_____________________________________________________________________________
e9f7cb15 729AliESDtrack *AlidNdPtHelper::GetTrackTrackSPDvtx(const AliESDEvent* esdEvent, Int_t iTrack, Bool_t bUpdate)
791aaf54 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
55468faf 774//_____________________________________________________________________________
775Bool_t AlidNdPtHelper::SelectEvent(const AliESDEvent* const esdEvent, AliESDtrackCuts* const esdTrackCuts) {
776// select events with at least
777// one reconstructed primary track in acceptance
778// pT>0.5 GeV/c, |eta|<0.8 for cross section studies
779
780if(!esdEvent) return kFALSE;
781if(!esdTrackCuts) return kFALSE;
782
783 AliESDtrack *track=0;
784 Int_t count = 0;
785 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
786 {
787 track = esdEvent->GetTrack(iTrack);
788 if(!track) continue;
789 if(track->Charge()==0) continue;
790 if(!esdTrackCuts->AcceptTrack(track)) continue;
791 if(track->Pt() < 0.5) continue;
792 if(TMath::Abs(track->Eta()) > 0.8) continue;
793
794 count++;
795 }
796
797 if(count > 0) return kTRUE;
798 else return kFALSE;
799
800return kFALSE;
801}
802
803//_____________________________________________________________________________
804Bool_t AlidNdPtHelper::SelectMCEvent(AliMCEvent* const mcEvent) {
805//
806// select events with at least
807// one prompt (MC primary) track in acceptance
808// pT>0.5 GeV/c, |eta|<0.8 for cross section studies
809//
810
811if(!mcEvent) return kFALSE;
812AliStack* stack = mcEvent->Stack();
813if(!stack) return kFALSE;
814
815 Int_t count = 0;
816 for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
817 {
818 TParticle* particle = stack->Particle(iMc);
819 if (!particle) continue;
820
821 // only charged particles
822 if(!particle->GetPDG()) continue;
823 Double_t charge = particle->GetPDG()->Charge()/3.;
824 if(charge == 0) continue;
825
826 // physical primary
827 Bool_t prim = stack->IsPhysicalPrimary(iMc);
828 if(!prim) continue;
829
830 if(particle->Pt() < 0.5) continue;
831 if(TMath::Abs(particle->Eta()) > 0.8) continue;
832
833 count++;
834 }
835
836 if(count > 0) return kTRUE;
837 else return kFALSE;
838
839return kFALSE;
840}
841
0aaa8b91 842//_____________________________________________________________________________
e9f7cb15 843Int_t AlidNdPtHelper::GetTPCMBTrackMult(const AliESDEvent *const esdEvent,const AlidNdPtEventCuts *const evtCuts, const AlidNdPtAcceptanceCuts *const accCuts,const AliESDtrackCuts *const trackCuts)
0aaa8b91 844{
845 //
846 // get MB event track multiplicity
847 //
848 if(!esdEvent)
849 {
850 ::Error("AlidNdPtHelper::GetTPCMBTrackMult()","esd event is NULL");
851 return 0;
852 }
853
854 if(!evtCuts || !accCuts || !trackCuts)
855 {
856 ::Error("AlidNdPtHelper::GetTPCMBTrackMult()","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 Int_t ntracks = esdEvent->GetNumberOfTracks();
871 Double_t dca[2],cov[3];
872 Int_t mult=0;
873 for (Int_t i=0;i <ntracks; i++){
874 AliESDtrack *t = esdEvent->GetTrack(i);
875 if (!t) continue;
876 if (t->Charge() == 0) continue;
877 if (!t->GetTPCInnerParam()) continue;
878 if (t->GetTPCNcls()<minTPCClust) continue;
879 //
774986ef 880 Double_t x[3]; t->GetXYZ(x);
881 Double_t b[3]; AliTracker::GetBxByBz(x,b);
882 const Double_t kMaxStep = 1; //max step over the material
883
0aaa8b91 884 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
774986ef 885 if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov))
0aaa8b91 886 {
887 if(tpcTrack) delete tpcTrack;
888 continue;
889 }
890 //
891 if (TMath::Abs(dca[0])>maxDCAr || TMath::Abs(dca[1])>maxDCAz) {
892 if(tpcTrack) delete tpcTrack;
893 continue;
894 }
895
896 mult++;
897
898 if(tpcTrack) delete tpcTrack;
899 }
900
901return mult;
902}
903
904//_____________________________________________________________________________
e9f7cb15 905Int_t AlidNdPtHelper::GetTPCMBPrimTrackMult(const AliESDEvent *const esdEvent, AliStack *const stack, const AlidNdPtEventCuts *const evtCuts, const AlidNdPtAcceptanceCuts *const accCuts, const AliESDtrackCuts *const trackCuts)
0aaa8b91 906{
907 //
908 // get MB primary event track multiplicity
909 //
910 if(!esdEvent)
911 {
912 ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","esd event is NULL");
913 return 0;
914 }
915
916 if(!stack)
917 {
918 ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","esd event is NULL");
919 return 0;
920 }
921
922 if(!evtCuts || !accCuts || !trackCuts)
923 {
924 ::Error("AlidNdPtHelper::GetTPCMBPrimTrackMult()","cuts not available");
925 return 0;
926 }
927
928 //
929 Double_t pos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
930 Double_t err[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
931 AliESDVertex vtx0(pos,err);
932
933 //
934 Float_t maxDCAr = accCuts->GetMaxDCAr();
935 Float_t maxDCAz = accCuts->GetMaxDCAz();
936 Float_t minTPCClust = trackCuts->GetMinNClusterTPC();
937
938 //
939 Int_t ntracks = esdEvent->GetNumberOfTracks();
940 Double_t dca[2],cov[3];
941 Int_t mult=0;
942 for (Int_t i=0;i <ntracks; i++){
943 AliESDtrack *t = esdEvent->GetTrack(i);
944 if (!t) continue;
945 if (t->Charge() == 0) continue;
946 if (!t->GetTPCInnerParam()) continue;
947 if (t->GetTPCNcls()<minTPCClust) continue;
948 //
774986ef 949 Double_t x[3]; t->GetXYZ(x);
950 Double_t b[3]; AliTracker::GetBxByBz(x,b);
951 const Double_t kMaxStep = 1; //max step over the material
952
0aaa8b91 953 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
774986ef 954 if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov))
0aaa8b91 955 {
956 if(tpcTrack) delete tpcTrack;
957 continue;
958 }
959 //
960 if (TMath::Abs(dca[0])>maxDCAr || TMath::Abs(dca[1])>maxDCAz) {
961 if(tpcTrack) delete tpcTrack;
962 continue;
963 }
964
965 Int_t label = TMath::Abs(t->GetLabel());
966 TParticle *part = stack->Particle(label);
967 if(!part) {
968 if(tpcTrack) delete tpcTrack;
969 continue;
970 }
971 if(!stack->IsPhysicalPrimary(label))
972 {
973 if(tpcTrack) delete tpcTrack;
974 continue;
975 }
976
977 mult++;
978
979 if(tpcTrack) delete tpcTrack;
980 }
981
982return mult;
983}
984
985
791aaf54 986//_____________________________________________________________________________
987Double_t AlidNdPtHelper::GetStrangenessCorrFactor(const Double_t pt)
988{
989// data driven correction factor for secondaries
990// underestimated secondaries with strangeness in Pythia (A. Dainese)
0aaa8b91 991
791aaf54 992//
993// pt=0.17; fact=1
994// pt=0.4; fact=1.07
995// pt=0.6; fact=1.25
996// pt>=1.2; fact=1.5
997//
0aaa8b91 998
791aaf54 999if (pt <= 0.17) return 1.0;
1000if (pt <= 0.4) return GetLinearInterpolationValue(0.17,1.0,0.4,1.07, pt);
1001if (pt <= 0.6) return GetLinearInterpolationValue(0.4,1.07,0.6,1.25, pt);
1002if (pt <= 1.2) return GetLinearInterpolationValue(0.6,1.25,1.2,1.5, pt);
1003return 1.5;
1004
1005}
1006
d4640855 1007//_____________________________________________________________________________
1008Double_t AlidNdPtHelper::GetStrangenessCorrFactorPbPb(const Double_t pt)
1009{
1010// data driven correction factor for secondaries (PbPb)
1011
1012if (pt <= 0.25) return 1.0;
1013if (pt <= 0.5) return GetLinearInterpolationValue(0.25,1.0,0.5,1.4, pt);
1014if (pt <= 1.0) return GetLinearInterpolationValue(0.5,1.4,1.0,1.47, pt);
1015if (pt <= 2.0) return GetLinearInterpolationValue(1.0,1.47,2.0,1.56, pt);
1016if (pt <= 5.0) return GetLinearInterpolationValue(2.0,1.56,5.0,1.67, pt);
1017return 1.67;
1018
1019}
1020
791aaf54 1021//___________________________________________________________________________
e9f7cb15 1022Double_t AlidNdPtHelper::GetLinearInterpolationValue(const Double_t x1,const Double_t y1,const Double_t x2,const Double_t y2, const Double_t pt)
791aaf54 1023{
1024//
1025// linear interpolation
1026//
1027 return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
1028}
0aaa8b91 1029
1030//_____________________________________________________________________________
f2dec884 1031Int_t AlidNdPtHelper::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AlidNdPtEventCuts *const evtCuts, AlidNdPtAcceptanceCuts *const accCuts)
0aaa8b91 1032{
1033 //
1034 // calculate mc event true track multiplicity
1035 //
1036 if(!mcEvent) return 0;
1037
1038 AliStack* stack = 0;
1039 Int_t mult = 0;
1040
1041 // MC particle stack
1042 stack = mcEvent->Stack();
1043 if (!stack) return 0;
1044
695facdf 1045 //
1046 //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
1047 //
1048
0aaa8b91 1049 Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
1050 if(!isEventOK) return 0;
1051
1052 Int_t nPart = stack->GetNtrack();
1053 for (Int_t iMc = 0; iMc < nPart; ++iMc)
1054 {
1055 TParticle* particle = stack->Particle(iMc);
1056 if (!particle)
1057 continue;
1058
1059 // only charged particles
09abaebb 1060 if(!particle->GetPDG()) continue;
0aaa8b91 1061 Double_t charge = particle->GetPDG()->Charge()/3.;
f2dec884 1062 if (TMath::Abs(charge) < 0.001)
0aaa8b91 1063 continue;
1064
1065 // physical primary
1066 Bool_t prim = stack->IsPhysicalPrimary(iMc);
1067 if(!prim) continue;
1068
774986ef 1069 // checked accepted without pt cut
1070 //if(accCuts->AcceptTrack(particle))
1071 if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
0aaa8b91 1072 {
1073 mult++;
1074 }
1075 }
1076
1077return mult;
1078}
1079
1080//_______________________________________________________________________
f2dec884 1081void AlidNdPtHelper::PrintMCInfo(AliStack *const pStack,Int_t label)
0aaa8b91 1082{
1083// print information about particles in the stack
1084
1085 if(!pStack)return;
1086 label = TMath::Abs(label);
1087 TParticle *part = pStack->Particle(label);
1088 Printf("########################");
1089 Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
1090 part->Print();
1091 TParticle* mother = part;
1092 Int_t imo = part->GetFirstMother();
1093 Int_t nprim = pStack->GetNprimary();
1094
1095 while((imo >= nprim)) {
1096 mother = pStack->Particle(imo);
1097 Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
1098 mother->Print();
1099 imo = mother->GetFirstMother();
1100 }
1101
1102 Printf("########################");
1103}
1104
1105
1106//_____________________________________________________________________________
f2dec884 1107TH1* AlidNdPtHelper::GetContCorrHisto(TH1 *const hist)
0aaa8b91 1108{
1109//
1110// get contamination histogram
1111//
1112 if(!hist) return 0;
1113
1114 Int_t nbins = hist->GetNbinsX();
f2dec884 1115 TH1 *hCont = (TH1D *)hist->Clone();
0aaa8b91 1116
1117 for(Int_t i=0; i<=nbins+1; i++) {
1118 Double_t binContent = hist->GetBinContent(i);
1119 Double_t binError = hist->GetBinError(i);
1120
f2dec884 1121 hCont->SetBinContent(i,1.-binContent);
1122 hCont->SetBinError(i,binError);
0aaa8b91 1123 }
1124
f2dec884 1125return hCont;
0aaa8b91 1126}
1127
1128
1129//_____________________________________________________________________________
f2dec884 1130TH1* AlidNdPtHelper::ScaleByBinWidth(TH1 *const hist)
0aaa8b91 1131{
1132//
1133// scale by bin width
1134//
1135 if(!hist) return 0;
1136
f2dec884 1137 TH1 *hScale = (TH1D *)hist->Clone();
1138 hScale->Scale(1.,"width");
0aaa8b91 1139
f2dec884 1140return hScale;
0aaa8b91 1141}
1142
1143//_____________________________________________________________________________
e9f7cb15 1144TH1* AlidNdPtHelper::CalcRelativeDifference(const TH1 *const hist1, const TH1 *const hist2)
0aaa8b91 1145{
1146//
1147// calculate rel. difference
1148//
1149
1150 if(!hist1) return 0;
1151 if(!hist2) return 0;
1152
f2dec884 1153 TH1 *h1Clone = (TH1D *)hist1->Clone();
1154 h1Clone->Sumw2();
0aaa8b91 1155
1156 // (rec-mc)/mc
f2dec884 1157 h1Clone->Add(hist2,-1);
1158 h1Clone->Divide(hist2);
0aaa8b91 1159
f2dec884 1160return h1Clone;
0aaa8b91 1161}
1162
1163//_____________________________________________________________________________
e9f7cb15 1164TH1* AlidNdPtHelper::CalcRelativeDifferenceFun(const TH1 *const hist1, TF1 *const fun)
0aaa8b91 1165{
1166//
08a80bfe 1167// calculate rel. difference
0aaa8b91 1168// between histogram and function
1169//
1170 if(!hist1) return 0;
1171 if(!fun) return 0;
1172
f2dec884 1173 TH1 *h1Clone = (TH1D *)hist1->Clone();
1174 h1Clone->Sumw2();
0aaa8b91 1175
1176 //
f2dec884 1177 h1Clone->Add(fun,-1);
1178 h1Clone->Divide(hist1);
0aaa8b91 1179
f2dec884 1180return h1Clone;
0aaa8b91 1181}
1182
1183//_____________________________________________________________________________
e9f7cb15 1184TH1* AlidNdPtHelper::NormalizeToEvent(const TH2 *const hist1, const TH1 *const hist2)
0aaa8b91 1185{
1186// normalise to event for a given multiplicity bin
1187// return pt histogram
1188
1189 if(!hist1) return 0;
1190 if(!hist2) return 0;
1191 char name[256];
1192
1193 Int_t nbinsX = hist1->GetNbinsX();
1194 //Int_t nbinsY = hist1->GetNbinsY();
1195
f2dec884 1196 TH1D *histNorm = 0;
0aaa8b91 1197 for(Int_t i=0; i<=nbinsX+1; i++) {
1198 sprintf(name,"mom_%d",i);
1199 TH1D *hist = (TH1D*)hist1->ProjectionY(name,i+1,i+1);
1200
1201 sprintf(name,"mom_norm");
1202 if(i==0) {
f2dec884 1203 histNorm = (TH1D *)hist->Clone(name);
1204 histNorm->Reset();
0aaa8b91 1205 }
1206
1207 Double_t nbEvents = hist2->GetBinContent(i);
1208 if(!nbEvents) { nbEvents = 1.; };
1209
1210 hist->Scale(1./nbEvents);
f2dec884 1211 histNorm->Add(hist);
0aaa8b91 1212 }
1213
f2dec884 1214return histNorm;
0aaa8b91 1215}
1216
1217//_____________________________________________________________________________
e9f7cb15 1218//THnSparse* AlidNdPtHelper::GenerateCorrMatrix(THnSparse *const hist1, const THnSparse *const hist2, char *const name) {
1219THnSparse* AlidNdPtHelper::GenerateCorrMatrix(THnSparse *const hist1, const THnSparse *const hist2, const char *name) {
0aaa8b91 1220// generate correction matrix
1221if(!hist1 || !hist2) return 0;
1222
1223THnSparse *h =(THnSparse*)hist1->Clone(name);;
1224h->Divide(hist1,hist2,1,1,"B");
1225
1226return h;
1227}
1228
1229//_____________________________________________________________________________
a9b12ce1 1230//TH2* AlidNdPtHelper::GenerateCorrMatrix(TH2 *const hist1, TH2 *const hist2, char *const name) {
1231TH2* AlidNdPtHelper::GenerateCorrMatrix(TH2 *const hist1, TH2 *const hist2, const char *name) {
0aaa8b91 1232// generate correction matrix
1233if(!hist1 || !hist2) return 0;
1234
1235TH2D *h =(TH2D*)hist1->Clone(name);;
1236h->Divide(hist1,hist2,1,1,"B");
1237
1238return h;
1239}
1240
1241//_____________________________________________________________________________
a9b12ce1 1242//TH1* AlidNdPtHelper::GenerateCorrMatrix(TH1 *const hist1, TH1 *const hist2, char *const name) {
1243TH1* AlidNdPtHelper::GenerateCorrMatrix(TH1 *const hist1, TH1 *const hist2, const char* name) {
0aaa8b91 1244// generate correction matrix
1245if(!hist1 || !hist2) return 0;
1246
1247TH1D *h =(TH1D*)hist1->Clone(name);;
1248h->Divide(hist1,hist2,1,1,"B");
1249
1250return h;
1251}
1252
1253//_____________________________________________________________________________
a9b12ce1 1254//THnSparse* AlidNdPtHelper::GenerateContCorrMatrix(THnSparse *const hist1, THnSparse *const hist2, char *const name) {
e9f7cb15 1255THnSparse* AlidNdPtHelper::GenerateContCorrMatrix(THnSparse *const hist1, const THnSparse *const hist2, const char* name) {
0aaa8b91 1256// generate contamination correction matrix
1257if(!hist1 || !hist2) return 0;
1258
1259THnSparse *hist = GenerateCorrMatrix(hist1, hist2, name);
1260if(!hist) return 0;
1261
1262// only for non ZERO bins!!!!
1263
1264Int_t* coord = new Int_t[hist->GetNdimensions()];
1265memset(coord, 0, sizeof(Int_t) * hist->GetNdimensions());
1266
1267 for (Long64_t i = 0; i < hist->GetNbins(); ++i) {
1268 Double_t v = hist->GetBinContent(i, coord);
1269 hist->SetBinContent(coord, 1.0-v);
1270 //printf("v %f, hist->GetBinContent(i, coord) %f \n",v,hist->GetBinContent(i, coord));
1271 Double_t err = hist->GetBinError(coord);
1272 hist->SetBinError(coord, err);
1273 }
1274
1275delete [] coord;
1276
1277return hist;
1278}
1279
1280//_____________________________________________________________________________
a9b12ce1 1281//TH2* AlidNdPtHelper::GenerateContCorrMatrix(TH2 *const hist1, TH2 *const hist2, char *const name) {
1282TH2* AlidNdPtHelper::GenerateContCorrMatrix(TH2 *const hist1, TH2 *const hist2, const char* name) {
0aaa8b91 1283// generate contamination correction matrix
1284if(!hist1 || !hist2) return 0;
1285
1286TH2 *hist = GenerateCorrMatrix(hist1, hist2, name);
1287if(!hist) return 0;
1288
1289Int_t nBinsX = hist->GetNbinsX();
1290Int_t nBinsY = hist->GetNbinsY();
1291
1292 for (Int_t i = 0; i < nBinsX+1; i++) {
1293 for (Int_t j = 0; j < nBinsY+1; j++) {
1294 Double_t cont = hist->GetBinContent(i,j);
1295 hist->SetBinContent(i,j,1.-cont);
1296 Double_t err = hist->GetBinError(i,j);
1297 hist->SetBinError(i,j,err);
1298 }
1299 }
1300
1301return hist;
1302}
1303
1304//_____________________________________________________________________________
a9b12ce1 1305//TH1* AlidNdPtHelper::GenerateContCorrMatrix(TH1 *const hist1, TH1 *const hist2, char *const name) {
1306TH1* AlidNdPtHelper::GenerateContCorrMatrix(TH1 *const hist1, TH1 *const hist2, const char* name) {
0aaa8b91 1307// generate contamination correction matrix
1308if(!hist1 || !hist2) return 0;
1309
1310TH1 *hist = GenerateCorrMatrix(hist1, hist2, name);
1311if(!hist) return 0;
1312
1313Int_t nBinsX = hist->GetNbinsX();
1314
1315 for (Int_t i = 0; i < nBinsX+1; i++) {
1316 Double_t cont = hist->GetBinContent(i);
1317 hist->SetBinContent(i,1.-cont);
1318 Double_t err = hist->GetBinError(i);
1319 hist->SetBinError(i,err);
1320 }
1321
1322return hist;
1323}
1324
1325//_____________________________________________________________________________
e9f7cb15 1326const AliESDVertex* AlidNdPtHelper::GetTPCVertexZ(const AliESDEvent* const esdEvent, const AlidNdPtEventCuts *const evtCuts, const AlidNdPtAcceptanceCuts *const accCuts, const AliESDtrackCuts *const trackCuts, Float_t fraction, Int_t ntracksMin){
0aaa8b91 1327 //
1328 // TPC Z vertexer
1329 //
bad4ba69 1330 if(!esdEvent)
1331 {
1332 ::Error("AlidNdPtHelper::GetTPCVertexZ()","cuts not available");
1333 return NULL;
1334 }
1335
1336 if(!evtCuts || !accCuts || !trackCuts)
1337 {
1338 ::Error("AlidNdPtHelper::GetTPCVertexZ()","cuts not available");
1339 return NULL;
1340 }
1341
1342 Double_t vtxpos[3]={evtCuts->GetMeanXv(),evtCuts->GetMeanYv(),evtCuts->GetMeanZv()};
1343 Double_t vtxsigma[3]={evtCuts->GetSigmaMeanXv(),evtCuts->GetSigmaMeanYv(),evtCuts->GetSigmaMeanZv()};
0aaa8b91 1344 AliESDVertex vtx0(vtxpos,vtxsigma);
bad4ba69 1345
1346 Double_t maxDCAr = accCuts->GetMaxDCAr();
1347 Double_t maxDCAz = accCuts->GetMaxDCAz();
1348 Int_t minTPCClust = trackCuts->GetMinNClusterTPC();
1349
0aaa8b91 1350 //
1351 Int_t ntracks = esdEvent->GetNumberOfTracks();
1352 TVectorD ztrack(ntracks);
0aaa8b91 1353 Double_t dca[2],cov[3];
1354 Int_t counter=0;
1355 for (Int_t i=0;i <ntracks; i++){
1356 AliESDtrack *t = esdEvent->GetTrack(i);
1357 if (!t) continue;
1358 if (!t->GetTPCInnerParam()) continue;
bad4ba69 1359 if (t->GetTPCNcls()<minTPCClust) continue;
0aaa8b91 1360 //
774986ef 1361
1362 Double_t x[3]; t->GetXYZ(x);
1363 Double_t b[3]; AliTracker::GetBxByBz(x,b);
1364 const Double_t kMaxStep = 1; //max step over the material
1365
0aaa8b91 1366 AliExternalTrackParam *tpcTrack = new AliExternalTrackParam(*(t->GetTPCInnerParam()));
774986ef 1367 if (!tpcTrack->PropagateToDCABxByBz(&vtx0,b,kMaxStep,dca,cov)) continue;
0aaa8b91 1368
1369 //
bad4ba69 1370 if (TMath::Abs(dca[0])>maxDCAr) continue;
1371 //if (TMath::Sqrt(cov[0])>sigmaXYcut) continue;
1372 if (TMath::Abs(tpcTrack->GetZ())>maxDCAz) continue;
0aaa8b91 1373
0aaa8b91 1374 ztrack[counter]=tpcTrack->GetZ();
1375 counter++;
1376
1377 if(tpcTrack) delete tpcTrack;
1378 }
bad4ba69 1379
0aaa8b91 1380 //
1381 // Find LTM z position
1382 //
1383 Double_t mean=0, sigma=0;
1384 if (counter<ntracksMin) return 0;
1385 //
1386 Int_t nused = TMath::Nint(counter*fraction);
1387 if (nused==counter) nused=counter-1;
1388 if (nused>1){
1389 AliMathBase::EvaluateUni(counter, ztrack.GetMatrixArray(), mean,sigma, TMath::Nint(counter*fraction));
1390 sigma/=TMath::Sqrt(nused);
1391 }else{
1392 mean = TMath::Mean(counter, ztrack.GetMatrixArray());
1393 sigma = TMath::RMS(counter, ztrack.GetMatrixArray());
1394 sigma/=TMath::Sqrt(counter-1);
1395 }
1396 vtxpos[2]=mean;
1397 vtxsigma[2]=sigma;
1398 const AliESDVertex* vertex = new AliESDVertex(vtxpos, vtxsigma);
1399 return vertex;
1400}
1401
1402//_____________________________________________________________________________
e9f7cb15 1403Int_t AlidNdPtHelper::GetSPDMBTrackMult(const AliESDEvent* const esdEvent, Float_t deltaThetaCut, Float_t deltaPhiCut)
0aaa8b91 1404{
1405 //
1406 // SPD track multiplicity
1407 //
1408
1409 // get tracklets
1410 const AliMultiplicity* mult = esdEvent->GetMultiplicity();
1411 if (!mult)
1412 return 0;
1413
1414 // get multiplicity from SPD tracklets
1415 Int_t inputCount = 0;
1416 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
1417 {
1418 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
1419
1420 Float_t phi = mult->GetPhi(i);
1421 if (phi < 0)
1422 phi += TMath::Pi() * 2;
1423 Float_t deltaPhi = mult->GetDeltaPhi(i);
1424 Float_t deltaTheta = mult->GetDeltaTheta(i);
1425
1426 if (TMath::Abs(deltaPhi) > 1)
1427 printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
1428
1429 if (deltaThetaCut > 0. && TMath::Abs(deltaTheta) > deltaThetaCut)
1430 continue;
1431
1432 if (deltaPhiCut > 0. && TMath::Abs(deltaPhi) > deltaPhiCut)
1433 continue;
1434
1435 ++inputCount;
1436 }
1437
1438return inputCount;
1439}
1440
1441//_____________________________________________________________________________
e9f7cb15 1442Int_t AlidNdPtHelper::GetSPDMBPrimTrackMult(const AliESDEvent* const esdEvent, AliStack* const stack, Float_t deltaThetaCut, Float_t deltaPhiCut)
0aaa8b91 1443{
1444 //
1445 // SPD track multiplicity
1446 //
1447
1448 // get tracklets
1449 const AliMultiplicity* mult = esdEvent->GetMultiplicity();
1450 if (!mult)
1451 return 0;
1452
1453 // get multiplicity from SPD tracklets
1454 Int_t inputCount = 0;
1455 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
1456 {
1457 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
1458
1459 Float_t phi = mult->GetPhi(i);
1460 if (phi < 0)
1461 phi += TMath::Pi() * 2;
1462 Float_t deltaPhi = mult->GetDeltaPhi(i);
1463 Float_t deltaTheta = mult->GetDeltaTheta(i);
1464
1465 if (TMath::Abs(deltaPhi) > 1)
1466 printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
1467
1468 if (deltaThetaCut > 0. && TMath::Abs(deltaTheta) > deltaThetaCut)
1469 continue;
1470
1471 if (deltaPhiCut > 0. && TMath::Abs(deltaPhi) > deltaPhiCut)
1472 continue;
1473
1474
1475 if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) ||
1476 !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
1477 continue;
1478
1479
1480 ++inputCount;
1481 }
1482
1483return inputCount;
1484}
fc98fbb5 1485//_____________________________________________________________________________
1486
1487THnSparse* AlidNdPtHelper::RebinTHnSparse(const THnSparse* hist1, THnSparse* hist2, const Char_t* newname, Option_t* option)
1488{
1489 THnSparse* htemp = 0;
1490 const THnSparse* hist = 0;
1491 TString opt = option;
1492 opt.ToLower();
1493 Bool_t calcErrors = kFALSE;
1494 Bool_t useRange = kFALSE;
1495 Bool_t overwrite = kFALSE;
1496 if (opt.Contains("e")) { calcErrors = kTRUE; } // calcluate correct errors (not implemented)
1497 if (opt.Contains("r")) { useRange = kTRUE; } // use the axis range given in hist1
1498 if (opt.Contains("o")) { overwrite = kTRUE; } // overwrite hist2 instead of creating a new one
1499 Int_t ndim = hist1->GetNdimensions();
1500 if (ndim != hist2->GetNdimensions()) {
1501 printf("AlidNdPtHelper::RebinTHnSparse: ERROR: Histograms have different dimensions \n");
1502 return 0;
1503 }
1504 Int_t* dims = new Int_t[ndim];
1505 for (Int_t i = 0; i < ndim; i++) { dims[i] = i; }
1506 if (useRange) {
1507 htemp = hist1->Projection(ndim,dims,"e");
1508 hist = htemp;
1509 } else { hist = hist1; }
1510 //THnSparse* hnew = hist2->Projection(ndim,dims,"o");
1511 //hnew->SetName(newname);
1512 THnSparse* hnew = 0;
1513 if (overwrite) {
1514 hnew = hist2;
1515 } else {
1516 hnew = (THnSparse*) hist2->Clone(newname);
1517 }
1518 for (Int_t i = 0; i < ndim; i++) { hnew->GetAxis(i)->SetRange(); }
1519 hnew->SetTitle(hist1->GetTitle());
1520 hnew->Reset();
1521 hnew->Sumw2();
1522 Double_t content;
1523 Double_t error;
1524 Int_t* c = new Int_t[ndim];
1525 Double_t* x = new Double_t[ndim];
1526 Long64_t n = hist->GetNbins();
1527 for (Long64_t j = 0; j < n; j++) {
1528 content = hist->GetBinContent(j,c);
1529 error = hist->GetBinError(j);
1530 for (Int_t i = 0; i < ndim; i++) {
1531 x[i] = hist->GetAxis(i)->GetBinCenter(c[i]);
1532 }
1533 /* function IsInRange is protected, shit!
1534 if (useRange) {
1535 if (! hist1->IsInRange(c)) continue;
1536 }
1537 */
1538 if (calcErrors) {
1539 // implementation to be done
1540 } else {
1541 hnew->Fill(x,content);
1542 }
1543 }
1544 delete[] c; c=0;
1545 delete[] x; x=0;
1546 delete[] dims; dims=0;
1547 if (htemp) { delete htemp; htemp = 0;}
1548 return hnew;
1549}
0aaa8b91 1550
1551