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