]>
Commit | Line | Data |
---|---|---|
600b313f | 1 | // |
2 | // Calculate the corrections in the central regions | |
3 | // | |
4 | // Inputs: | |
5 | // - AliESDEvent | |
6 | // | |
7 | // Outputs: | |
8 | // - AliAODCentralMult | |
9 | // | |
10 | // Histograms | |
11 | // | |
12 | // Corrections used | |
13 | // | |
14 | #include "AliCentralMCCorrectionsTask.h" | |
c8b1a7db | 15 | #include "AliCentralCorrectionManager.h" |
600b313f | 16 | #include "AliTriggerAnalysis.h" |
17 | #include "AliPhysicsSelection.h" | |
18 | #include "AliLog.h" | |
19 | #include "AliHeader.h" | |
20 | #include "AliGenEventHeader.h" | |
21 | #include "AliESDEvent.h" | |
22 | #include "AliAODHandler.h" | |
23 | #include "AliMultiplicity.h" | |
24 | #include "AliInputEventHandler.h" | |
25 | #include "AliStack.h" | |
26 | #include "AliMCEvent.h" | |
27 | #include "AliAODForwardMult.h" | |
28 | #include "AliCentralCorrSecondaryMap.h" | |
ea6f96fb | 29 | #include "AliCentralCorrAcceptance.h" |
f53fb4f6 | 30 | #include "AliForwardUtil.h" |
ea6f96fb | 31 | #include "AliMultiplicity.h" |
600b313f | 32 | #include <TH1.h> |
33 | #include <TH2D.h> | |
34 | #include <TDirectory.h> | |
35 | #include <TList.h> | |
36 | #include <TROOT.h> | |
37 | #include <iostream> | |
38 | ||
600b313f | 39 | //==================================================================== |
40 | AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask() | |
c8b1a7db | 41 | : AliBaseMCCorrectionsTask(), |
600b313f | 42 | fTrackDensity(), |
c8b1a7db | 43 | fSecCorr(0), |
44 | fAccCorr(0), | |
ea6f96fb | 45 | fNPhiBins(20), |
bfab35d9 | 46 | fEffectiveCorr(true), |
47 | fEtaCut(1.9), | |
48 | fCorrCut(0.6) | |
600b313f | 49 | { |
50 | // | |
51 | // Constructor | |
52 | // | |
53 | // Parameters: | |
54 | // name Name of task | |
55 | // | |
5ca83fee | 56 | DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask"); |
600b313f | 57 | } |
58 | ||
59 | //____________________________________________________________________ | |
60 | AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name) | |
c8b1a7db | 61 | : AliBaseMCCorrectionsTask(name, &(AliCentralCorrectionManager::Instance())), |
600b313f | 62 | fTrackDensity("trackDensity"), |
c8b1a7db | 63 | fSecCorr(0), |
64 | fAccCorr(0), | |
ea6f96fb | 65 | fNPhiBins(20), |
bfab35d9 | 66 | fEffectiveCorr(true), |
67 | fEtaCut(1.9), | |
68 | fCorrCut(0.6) | |
600b313f | 69 | { |
70 | // | |
71 | // Constructor | |
72 | // | |
73 | // Parameters: | |
74 | // name Name of task | |
75 | // | |
5ca83fee | 76 | DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name); |
600b313f | 77 | } |
78 | ||
79 | //____________________________________________________________________ | |
c8b1a7db | 80 | AliBaseMCCorrectionsTask::VtxBin* |
81 | AliCentralMCCorrectionsTask::CreateVtxBin(Double_t low, Double_t high) | |
600b313f | 82 | { |
c8b1a7db | 83 | return new AliCentralMCCorrectionsTask::VtxBin(low,high,fEtaAxis,fNPhiBins); |
600b313f | 84 | } |
85 | ||
86 | //____________________________________________________________________ | |
c8b1a7db | 87 | Bool_t |
88 | AliCentralMCCorrectionsTask::ProcessESD(const AliESDEvent& esd, | |
89 | const AliMCEvent& mc, | |
90 | AliBaseMCCorrectionsTask::VtxBin& bin, | |
91 | Double_t vz) | |
600b313f | 92 | { |
c8b1a7db | 93 | AliCentralMCCorrectionsTask::VtxBin& vb = |
94 | static_cast<AliCentralMCCorrectionsTask::VtxBin&>(bin); | |
600b313f | 95 | |
96 | // Now process our input data and store in own ESD object | |
c8b1a7db | 97 | fTrackDensity.Calculate(mc, vz, *vb.fHits, bin.fPrimary); |
98 | ||
ea6f96fb | 99 | // Get the ESD object |
c8b1a7db | 100 | const AliMultiplicity* spdmult = esd.GetMultiplicity(); |
101 | ||
ea6f96fb | 102 | // Count number of tracklets per bin |
103 | for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) | |
c8b1a7db | 104 | vb.fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j)); |
ea6f96fb | 105 | //...and then the unused clusters in layer 1 |
106 | for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) { | |
c8b1a7db | 107 | Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)); |
108 | vb.fClusters->Fill(eta, spdmult->GetPhiSingle(j)); | |
ea6f96fb | 109 | } |
c8b1a7db | 110 | |
ea6f96fb | 111 | // Count events |
c8b1a7db | 112 | bin.fCounts->Fill(0.5); |
bcd522ff | 113 | |
c8b1a7db | 114 | return true; |
600b313f | 115 | } |
600b313f | 116 | //____________________________________________________________________ |
117 | void | |
c8b1a7db | 118 | AliCentralMCCorrectionsTask::CreateCorrections(TList* results) |
600b313f | 119 | { |
c8b1a7db | 120 | fSecCorr = new AliCentralCorrSecondaryMap; |
121 | fSecCorr->SetVertexAxis(fVtxAxis); | |
600b313f | 122 | |
c8b1a7db | 123 | fAccCorr = new AliCentralCorrAcceptance; |
124 | fAccCorr->SetVertexAxis(fVtxAxis); | |
ea6f96fb | 125 | |
c8b1a7db | 126 | results->Add(fSecCorr); |
127 | results->Add(fAccCorr); | |
128 | } | |
600b313f | 129 | |
c8b1a7db | 130 | //____________________________________________________________________ |
131 | Bool_t | |
132 | AliCentralMCCorrectionsTask::FinalizeVtxBin(AliBaseMCCorrectionsTask::VtxBin* | |
133 | bin, UShort_t iVz) | |
134 | { | |
135 | ||
136 | AliCentralMCCorrectionsTask::VtxBin* vb = | |
137 | static_cast<AliCentralMCCorrectionsTask::VtxBin*>(bin); | |
138 | vb->Terminate(fList, fResults, iVz, fEffectiveCorr, | |
139 | fEtaCut, fCorrCut, fSecCorr,fAccCorr); | |
140 | return true; | |
600b313f | 141 | } |
c8b1a7db | 142 | |
600b313f | 143 | |
144 | //____________________________________________________________________ | |
145 | void | |
146 | AliCentralMCCorrectionsTask::Print(Option_t* option) const | |
147 | { | |
c8b1a7db | 148 | AliBaseMCCorrectionsTask::Print(option); |
149 | std::cout << " # of phi bins: " << fNPhiBins << "\n" | |
bfab35d9 | 150 | << " Effective corr.: " << fEffectiveCorr << "\n" |
151 | << " Eta cut-off: " << fEtaCut << "\n" | |
152 | << " Acceptance cut: " << fCorrCut | |
600b313f | 153 | << std::endl; |
154 | gROOT->IncreaseDirLevel(); | |
600b313f | 155 | fTrackDensity.Print(option); |
156 | gROOT->DecreaseDirLevel(); | |
157 | } | |
158 | ||
159 | //==================================================================== | |
600b313f | 160 | AliCentralMCCorrectionsTask::VtxBin::VtxBin() |
c8b1a7db | 161 | : AliBaseMCCorrectionsTask::VtxBin(), |
162 | fHits(0), | |
163 | fClusters(0) | |
600b313f | 164 | { |
165 | } | |
166 | //____________________________________________________________________ | |
167 | AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t low, | |
168 | Double_t high, | |
169 | const TAxis& axis, | |
170 | UShort_t nPhi) | |
c8b1a7db | 171 | : AliBaseMCCorrectionsTask::VtxBin(low, high, axis, nPhi), |
600b313f | 172 | fHits(0), |
c8b1a7db | 173 | fClusters(0) |
600b313f | 174 | { |
600b313f | 175 | fHits = static_cast<TH2D*>(fPrimary->Clone("hits")); |
176 | fHits->SetTitle("Hits"); | |
177 | fHits->SetDirectory(0); | |
178 | ||
ea6f96fb | 179 | fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters")); |
180 | fClusters->SetTitle("Clusters"); | |
181 | fClusters->SetDirectory(0); | |
600b313f | 182 | } |
183 | ||
600b313f | 184 | |
185 | //____________________________________________________________________ | |
c8b1a7db | 186 | TList* |
5934a3e3 | 187 | AliCentralMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l) |
600b313f | 188 | { |
c8b1a7db | 189 | TList* d = AliBaseMCCorrectionsTask::VtxBin::CreateOutputObjects(l); |
600b313f | 190 | |
191 | d->Add(fHits); | |
ea6f96fb | 192 | d->Add(fClusters); |
c8b1a7db | 193 | |
194 | return d; | |
600b313f | 195 | } |
196 | //____________________________________________________________________ | |
197 | void | |
5934a3e3 | 198 | AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input, |
bfab35d9 | 199 | TList* output, |
200 | UShort_t iVz, | |
201 | Bool_t effectiveCorr, | |
202 | Double_t etaCut, | |
203 | Double_t accCut, | |
204 | AliCentralCorrSecondaryMap* map, | |
205 | AliCentralCorrAcceptance* acorr) | |
600b313f | 206 | { |
207 | TList* out = new TList; | |
208 | out->SetName(GetName()); | |
209 | out->SetOwner(); | |
210 | output->Add(out); | |
bfab35d9 | 211 | |
600b313f | 212 | TList* l = static_cast<TList*>(input->FindObject(GetName())); |
213 | if (!l) { | |
214 | AliError(Form("List %s not found in %s", GetName(), input->GetName())); | |
215 | return; | |
216 | } | |
217 | ||
bfab35d9 | 218 | // Get the sums |
600b313f | 219 | TH2D* hits = static_cast<TH2D*>(l->FindObject("hits")); |
ea6f96fb | 220 | TH2D* clus = static_cast<TH2D*>(l->FindObject("clusters")); |
600b313f | 221 | TH2D* prim = static_cast<TH2D*>(l->FindObject("primary")); |
222 | if (!hits || !prim) { | |
223 | AliError(Form("Missing histograms: %p, %p", hits, prim)); | |
224 | return; | |
225 | } | |
226 | ||
bfab35d9 | 227 | // Clone cluster and hit map |
228 | TH2D* secMapEff = static_cast<TH2D*>(clus->Clone("secMapEff")); | |
229 | TH2D* secMapHit = static_cast<TH2D*>(hits->Clone("secMapHit")); | |
230 | secMapEff->SetTitle("2^{nd} map from clusters"); | |
231 | secMapEff->SetDirectory(0); | |
232 | secMapHit->SetTitle("2^{nd} map from MC hits"); | |
233 | secMapHit->SetDirectory(0); | |
234 | ||
235 | // Divide cluster and hit map with number of primaries | |
236 | secMapEff->Divide(prim); | |
237 | secMapHit->Divide(prim); | |
238 | ||
239 | // Create acceptance histograms | |
240 | TH1D* accEff = new TH1D("accEff", | |
241 | "Acceptance correction for SPD (from 2^{nd} map)" , | |
242 | fPrimary->GetXaxis()->GetNbins(), | |
243 | fPrimary->GetXaxis()->GetXmin(), | |
244 | fPrimary->GetXaxis()->GetXmax()); | |
245 | TH1D* accHit = static_cast<TH1D*>(accEff->Clone("accHit")); | |
246 | accHit->SetTitle("Acceptance correction for SPD (from clusters)"); | |
247 | ||
248 | // Diagnostics histogra, | |
249 | TH2* dia = static_cast<TH2D*>(clus->Clone("diagnostics")); | |
250 | dia->SetTitle("Scaled cluster density"); | |
251 | ||
252 | // Total number of channels along phi and # of eta bins | |
253 | Int_t nTotal = secMapHit->GetNbinsY(); | |
254 | Int_t nEta = secMapHit->GetNbinsX(); | |
255 | ||
256 | for(Int_t xx = 1; xx <= nEta; xx++) { | |
257 | Double_t eta = secMapHit->GetXaxis()->GetBinCenter(xx); | |
258 | Bool_t ins = TMath::Abs(eta) <= etaCut; | |
259 | Double_t mm = 0; | |
260 | if (ins) { | |
261 | // Find the maximum cluster signal in this phi range | |
262 | for (Int_t yy = 1; yy <= nTotal; yy++) { | |
263 | Double_t c = clus->GetBinContent(xx,yy); | |
264 | mm = TMath::Max(mm, c); | |
265 | } | |
266 | } | |
267 | // Count number of phi bins with enough clusters or high enough | |
268 | // correction. | |
269 | Int_t nOKEff = 0; | |
270 | Int_t nOKHit = 0; | |
271 | for(Int_t yy = 1; yy <=nTotal; yy++) { | |
272 | if (!ins) { // Not inside Eta cut | |
273 | secMapEff->SetBinContent(xx,yy,0.); | |
274 | secMapEff->SetBinError(xx,yy,0.); | |
275 | secMapHit->SetBinContent(xx,yy,0.); | |
276 | secMapHit->SetBinError(xx,yy,0.); | |
277 | dia->SetBinContent(xx,yy,0); | |
278 | continue; | |
279 | } | |
280 | ||
281 | // Check if the background correction is large enough, or zero map | |
282 | if(secMapEff->GetBinContent(xx,yy) > 0.9) { | |
283 | // acc->Fill(h->GetXaxis()->GetBinCenter(xx)); | |
284 | nOKEff++; | |
285 | } | |
286 | else { | |
287 | secMapEff->SetBinContent(xx,yy,0.); | |
288 | secMapEff->SetBinError(xx,yy,0.); | |
289 | } | |
290 | ||
291 | // Check if the number of cluster is large enough, or zero map | |
292 | Double_t c = clus->GetBinContent(xx,yy); | |
293 | Double_t s = (mm < 1e-6) ? 0 : c / mm; | |
294 | dia->SetBinContent(xx,yy,s); | |
295 | if (s >= accCut) { | |
296 | nOKHit++; | |
ea6f96fb | 297 | } |
ea6f96fb | 298 | else { |
bfab35d9 | 299 | secMapHit->SetBinContent(xx,yy,0); |
300 | secMapHit->SetBinError(xx,yy,0); | |
ea6f96fb | 301 | } |
ea6f96fb | 302 | } |
bfab35d9 | 303 | |
304 | // Calculate acceptance as ratio of bins with enough clusters and | |
305 | // total number of phi bins. | |
306 | Double_t accXX = float(nOKHit) / nTotal; | |
307 | if (accXX < 0.2) accXX = 0; | |
308 | accHit->SetBinContent(xx, accXX); | |
309 | ||
310 | // Calculate acceptance as ratio of bins with large enough | |
311 | // correction and total number of phi bins. | |
312 | accXX = float(nOKEff) / nTotal; | |
313 | if (accXX < 0.2) accXX = 0; | |
314 | accEff->SetBinContent(xx, accXX); | |
ea6f96fb | 315 | } |
600b313f | 316 | |
bfab35d9 | 317 | TH2D* secMap = (effectiveCorr ? secMapEff : secMapHit); |
318 | TH2D* secMapAlt = (effectiveCorr ? secMapHit : secMapEff); | |
319 | TH1D* acc = (effectiveCorr ? accEff : accHit); | |
320 | TH1D* accAlt = (effectiveCorr ? accHit : accEff); | |
321 | out->Add(secMap->Clone("secMap")); | |
322 | out->Add(secMapAlt->Clone()); | |
323 | out->Add(acc->Clone("acc")); | |
324 | out->Add(accAlt->Clone()); | |
325 | out->Add(dia->Clone()); | |
326 | ||
327 | map->SetCorrection(iVz, secMap); | |
ea6f96fb | 328 | acorr->SetCorrection(iVz, acc); |
600b313f | 329 | } |
330 | ||
331 | // | |
332 | // EOF | |
333 | // |