]>
Commit | Line | Data |
---|---|---|
8309c1ab | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
18 | /////////////////////////////////////////////////////////////////////////////// | |
19 | // // | |
20 | // class for ZDC reconstruction // | |
21 | // // | |
22 | /////////////////////////////////////////////////////////////////////////////// | |
23 | ||
24 | ||
25 | #include <TF1.h> | |
26 | ||
27 | #include "AliRunLoader.h" | |
28 | #include "AliRawReader.h" | |
29 | #include "AliESD.h" | |
30 | #include "AliZDCDigit.h" | |
31 | #include "AliZDCRawStream.h" | |
32 | #include "AliZDCReco.h" | |
33 | #include "AliZDCReconstructor.h" | |
34 | ||
35 | ||
36 | ClassImp(AliZDCReconstructor) | |
37 | ||
38 | ||
39 | //_____________________________________________________________________________ | |
40 | AliZDCReconstructor:: AliZDCReconstructor() | |
41 | { | |
42 | // default constructor | |
43 | ||
44 | // --- Number of generated spectator nucleons and impact parameter | |
45 | // -------------------------------------------------------------------------------------------------- | |
46 | // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001) | |
47 | /*// Fit results for neutrons (Nspectator n true vs. EZN) | |
48 | fZNCen = new TF1("fZNCen", | |
49 | "(-2.116909+sqrt(2.116909*2.116909-4*(-0.00651)*(14.556798-x)))/(2*(-0.00651))",0.,158.5); | |
50 | fZNPer = new TF1("fZNPer", | |
51 | "(-34.695134-sqrt(34.695134*34.695134-4*(-0.174780)*(-1562.283443-x)))/(2*(-0.174780))",0.,158.5); | |
52 | // Fit results for protons (Nspectator p true vs. EZP) | |
53 | fZPCen = new TF1("fZPCen", | |
54 | "(-1.3217+sqrt(1.3217*1.3217-4*(-0.007934)*(4.742873-x)))/(2*(-0.007934))",0.,58.91); | |
55 | fZPPer = new TF1("fZPPer", | |
56 | "(-15.788267-sqrt(15.788267*15.788267-4*(-0.133359)*(-383.800673-x)))/(2*(-0.133359))",0.,58.91); | |
57 | // Fit results for total number of spectators (Nspectators true vs. EZDC) | |
58 | fZDCCen = new TF1("fZDCCen", | |
59 | "(-1.867335+sqrt(1.867335*1.867335-4*(-0.004119)*(19.100289-x)))/(2*(-0.004119))",0.,220.4); | |
60 | fZDCPer = new TF1("fZDCPer", | |
61 | "(-22.429097-sqrt(22.429097*22.429097-4*(-0.072435)*(-1482.034526-x)))/(2*(-0.072435))",0.,220.4);*/ | |
62 | // -------------------------------------------------------------------------------------------------- | |
63 | // [1] ### Results from a new production -> 0<b<18 fm (Apr 2002) | |
64 | // Fit results for neutrons (Nspectator n true vs. EZN) | |
65 | fZNCen = new TF1("fZNCen", | |
66 | "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.); | |
67 | fZNPer = new TF1("fZNPer", | |
68 | "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.); | |
69 | // Fit results for protons (Nspectator p true vs. EZP) | |
70 | fZPCen = new TF1("fZPCen", | |
71 | "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.); | |
72 | fZPPer = new TF1("fZPPer", | |
73 | "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.); | |
74 | // Fit results for total number of spectators (Nspectators true vs. EZDC) | |
75 | fZDCCen = new TF1("fZDCCen", | |
76 | "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.); | |
77 | fZDCPer = new TF1("fZDCPer", | |
78 | "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.); | |
79 | // -------------------------------------------------------------------------------------------------- | |
80 | // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001) | |
81 | /*// Fit results for b (b vs. EZDC) | |
82 | //fbCen = new TF1("fbCen","0.611543+0.052231*x-0.000112*x*x+0.000000374*x*x*x",0.,222.); | |
83 | //fbPer = new TF1("fbPer","16.552010-0.023866*x-0.00001*x*x",0.,222.); | |
84 | fbCen = new TF1("fbCen","0.612769+0.051929*x-0.0001074*x*x+0.0000003724*x*x*x",0.,225.); | |
85 | fbPer = new TF1("fbPer","16.6131016-0.026053*x+0.000006893*x*x",0.,225.);*/ | |
86 | // -------------------------------------------------------------------------------------------------- | |
87 | // [2] ### Results from a new production -> 0<b<18 fm (Apr 2002) | |
88 | fbCen = new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.); | |
89 | fbPer = new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.); | |
90 | // -------------------------------------------------------------------------------------------------- | |
91 | // Evaluating Nspectators and b from ZEM energy | |
92 | // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001) | |
93 | /*fZEMn = new TF1("fZEMn","124.2-0.0566*x+0.000006014*x*x",0.,3500.); | |
94 | fZEMp = new TF1("fZEMp","81.3-0.03834*x+0.000004359*x*x",0.,3500.); | |
95 | fZEMsp = new TF1("fZEMsp","205.6-0.09567*x+0.00001056*x*x",0.,3500.); | |
96 | fZEMb = new TF1("fZEMb","15.8-0.02084*x+2.802e-5*x*x-2.007e-8*x*x*x+6.586e-12*x*x*x*x-8.042e-16*x*x*x*x*x",0.,3500.);*/ | |
97 | // -------------------------------------------------------------------------------------------------- | |
98 | // [2] ### Results from a new production -> 0<b<18 fm (Apr 2002) | |
99 | fZEMn = new TF1("fZEMn","126.2-0.05399*x+0.000005679*x*x",0.,4000.); | |
100 | fZEMp = new TF1("fZEMp","82.49-0.03611*x+0.00000385*x*x",0.,4000.); | |
101 | fZEMsp = new TF1("fZEMsp","208.7-0.09006*x+0.000009526*x*x",0.,4000.); | |
102 | fZEMb = new TF1("fZEMb","16.06-0.01633*x+1.44e-5*x*x-6.778e-9*x*x*x+1.438e-12*x*x*x*x-1.112e-16*x*x*x*x*x",0.,4000.); | |
103 | } | |
104 | ||
105 | //_____________________________________________________________________________ | |
106 | AliZDCReconstructor::AliZDCReconstructor(const AliZDCReconstructor& | |
107 | reconstructor): | |
108 | AliReconstructor(reconstructor) | |
109 | { | |
110 | // copy constructor | |
111 | ||
112 | Fatal("AliZDCReconstructor", "copy constructor not implemented"); | |
113 | } | |
114 | ||
115 | //_____________________________________________________________________________ | |
116 | AliZDCReconstructor& AliZDCReconstructor::operator = | |
117 | (const AliZDCReconstructor& /*reconstructor*/) | |
118 | { | |
119 | // assignment operator | |
120 | ||
121 | Fatal("operator =", "assignment operator not implemented"); | |
122 | return *this; | |
123 | } | |
124 | ||
125 | //_____________________________________________________________________________ | |
126 | AliZDCReconstructor::~AliZDCReconstructor() | |
127 | { | |
128 | // destructor | |
129 | ||
130 | delete fZNCen; | |
131 | delete fZNPer; | |
132 | delete fZPCen; | |
133 | delete fZPPer; | |
134 | delete fZDCCen; | |
135 | delete fZDCPer; | |
136 | delete fbCen; | |
137 | delete fbPer; | |
138 | delete fZEMn; | |
139 | delete fZEMp; | |
140 | delete fZEMsp; | |
141 | delete fZEMb; | |
142 | } | |
143 | ||
144 | ||
145 | //_____________________________________________________________________________ | |
146 | void AliZDCReconstructor::Reconstruct(AliRunLoader* runLoader) const | |
147 | { | |
148 | // local ZDC reconstruction | |
149 | ||
150 | AliLoader* loader = runLoader->GetLoader("ZDCLoader"); | |
151 | if (!loader) return; | |
152 | loader->LoadDigits("read"); | |
153 | loader->LoadRecPoints("recreate"); | |
154 | AliZDCDigit digit; | |
155 | AliZDCDigit* pdigit = &digit; | |
156 | ||
157 | // Event loop | |
158 | for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { | |
159 | runLoader->GetEvent(iEvent); | |
160 | ||
161 | // load digits | |
162 | loader->LoadDigits(); | |
163 | TTree* treeD = loader->TreeD(); | |
164 | if (!treeD) continue; | |
165 | treeD->SetBranchAddress("ZDC", &pdigit); | |
166 | ||
167 | // loop over digits | |
168 | Int_t znraw=0, zpraw=0, zemraw=0; | |
169 | for (Int_t iDigit = 0; iDigit < treeD->GetEntries(); iDigit++) { | |
170 | treeD->GetEntry(iDigit); | |
171 | if (!pdigit) continue; | |
172 | ||
173 | if(digit.GetSector(0) == 1) znraw += digit.GetADCValue(0); | |
174 | else if(digit.GetSector(0) == 2) zpraw += digit.GetADCValue(0); | |
175 | else if(digit.GetSector(0) == 3) zemraw += digit.GetADCValue(0); | |
176 | } | |
177 | ||
178 | // reconstruct the event | |
179 | ReconstructEvent(loader, znraw, zpraw, zemraw); | |
180 | } | |
181 | ||
182 | loader->UnloadDigits(); | |
183 | loader->UnloadRecPoints(); | |
184 | } | |
185 | ||
186 | //_____________________________________________________________________________ | |
187 | void AliZDCReconstructor::Reconstruct(AliRunLoader* runLoader, | |
188 | AliRawReader* rawReader) const | |
189 | { | |
190 | // local ZDC reconstruction for raw data | |
191 | ||
192 | AliLoader* loader = runLoader->GetLoader("ZDCLoader"); | |
193 | if (!loader) return; | |
194 | loader->LoadRecPoints("recreate"); | |
195 | // Event loop | |
196 | Int_t iEvent = 0; | |
197 | while (rawReader->NextEvent()) { | |
198 | runLoader->GetEvent(iEvent++); | |
199 | ||
200 | // loop over raw data digits | |
201 | Int_t znraw=0, zpraw=0, zemraw=0; | |
202 | AliZDCRawStream digit(rawReader); | |
203 | while (digit.Next()) { | |
204 | if(digit.IsADCDataWord()){ | |
205 | if(digit.GetADCGain() == 0){ | |
206 | printf("ADC value = %d\n",digit.GetADCValue()); | |
207 | if(digit.GetSector(0) == 1) znraw += digit.GetADCValue(); | |
208 | else if(digit.GetSector(0) == 2) zpraw += digit.GetADCValue(); | |
209 | else if(digit.GetSector(0) == 3) zemraw += digit.GetADCValue(); | |
210 | } | |
211 | } | |
212 | } | |
213 | // reconstruct the event | |
214 | ReconstructEvent(loader, znraw, zpraw, zemraw); | |
215 | } | |
216 | ||
217 | loader->UnloadRecPoints(); | |
218 | } | |
219 | ||
220 | //_____________________________________________________________________________ | |
221 | void AliZDCReconstructor::ReconstructEvent(AliLoader* loader, Int_t znraw, | |
222 | Int_t zpraw, Int_t zemraw) const | |
223 | { | |
224 | // reconstruct one event | |
225 | ||
226 | // if (GetDebug()) printf("\n --- znraw = %d, zpraw = %d, zemraw = %d\n",znraw, zpraw, zemraw); | |
227 | ||
228 | // --- Pedestal subtraction | |
229 | Int_t zncorr, zpcorr, zemcorr, meanPed=50; | |
230 | zncorr = znraw - 5*meanPed; | |
231 | zpcorr = zpraw - 5*meanPed; | |
232 | zemcorr = zemraw - 2*meanPed; | |
233 | if(zncorr<0) zncorr=0; | |
234 | if(zpcorr<0) zpcorr=0; | |
235 | if(zemcorr<0) zemcorr=0; | |
236 | // if (GetDebug()) printf("\n zncorr = %d, zpcorr = %d, zemcorr = %d\n",zncorr,zpcorr,zemcorr); | |
237 | ||
238 | // --- ADCchannel -> photoelectrons | |
239 | // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7) | |
240 | // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7) | |
241 | Float_t znphe, zpphe, zemphe, convFactor = 0.08; | |
242 | znphe = zncorr/convFactor; | |
243 | zpphe = zpcorr/convFactor; | |
244 | zemphe = zemcorr/convFactor; | |
245 | // if (GetDebug()) printf("\n znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe); | |
246 | ||
247 | // --- Energy calibration | |
248 | // Conversion factors for hadronic ZDCs goes from phe yield to TRUE | |
249 | // incident energy (conversion from GeV to TeV is included); while for EM | |
250 | // calos conversion is from light yield to detected energy calculated by | |
251 | // GEANT NB -> ZN and ZP conversion factors are constant since incident | |
252 | // spectators have all the same energy, ZEM energy is obtained through a | |
253 | // fit over the whole range of incident particle energies | |
254 | // (obtained with full HIJING simulations) | |
255 | Float_t znenergy, zpenergy, zemenergy, zdcenergy; | |
256 | Float_t znphexTeV=329., zpphexTeV=369.; | |
257 | znenergy = znphe/znphexTeV; | |
258 | zpenergy = zpphe/zpphexTeV; | |
259 | zdcenergy = znenergy+zpenergy; | |
260 | zemenergy = -4.81+0.3238*zemphe; | |
261 | if(zemenergy<0) zemenergy=0; | |
262 | // if (GetDebug()) printf(" znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, " | |
263 | // "\n zemenergy = %f TeV\n", znenergy, zpenergy, | |
264 | // zdcenergy, zemenergy); | |
265 | ||
266 | // if(zdcenergy==0) | |
267 | // if (GetDebug()) printf("\n\n ### ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, " | |
268 | // " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy); | |
269 | ||
270 | // --- Number of incident spectator nucleons | |
271 | Int_t nDetSpecN, nDetSpecP; | |
272 | nDetSpecN = (Int_t) (znenergy/2.760); | |
273 | nDetSpecP = (Int_t) (zpenergy/2.760); | |
274 | // if (GetDebug()) printf("\n nDetSpecN = %d, nDetSpecP = %d\n",nDetSpecN, nDetSpecP); | |
275 | ||
276 | Int_t nGenSpecN=0, nGenSpecP=0, nGenSpec=0; | |
277 | Double_t impPar=0; | |
278 | // Cut value for Ezem (GeV) | |
279 | // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001) | |
280 | //Float_t eZEMCut = 360.; | |
281 | // [2] ### Results from a new production -> 0<b<18 fm (Apr 2002) | |
282 | Float_t eZEMCut = 420.; | |
283 | Float_t deltaEZEMSup = 690.; | |
284 | Float_t deltaEZEMInf = 270.; | |
285 | if(zemenergy > (eZEMCut+deltaEZEMSup)){ | |
286 | nGenSpecN = (Int_t) (fZNCen->Eval(znenergy)); | |
287 | nGenSpecP = (Int_t) (fZPCen->Eval(zpenergy)); | |
288 | nGenSpec = (Int_t) (fZDCCen->Eval(zdcenergy)); | |
289 | impPar = fbCen->Eval(zdcenergy); | |
290 | //printf(" fZNCen = %f, fZPCen = %f, fZDCCen = %f\n",fZNCen->Eval(znenergy), | |
291 | // fZPCen->Eval(zpenergy),fZDCCen->Eval(zdcenergy)); | |
292 | } | |
293 | else if(zemenergy < (eZEMCut-deltaEZEMInf)){ | |
294 | nGenSpecN = (Int_t) (fZNPer->Eval(znenergy)); | |
295 | nGenSpecP = (Int_t) (fZPPer->Eval(zpenergy)); | |
296 | nGenSpec = (Int_t) (fZDCPer->Eval(zdcenergy)); | |
297 | impPar = fbPer->Eval(zdcenergy); | |
298 | //printf(" fZNPer = %f, fZPPer = %f, fZDCPer = %f\n",fZNPer->Eval(znenergy), | |
299 | // fZPPer->Eval(zpenergy),fZDCPer->Eval(zdcenergy)); | |
300 | } | |
301 | else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){ | |
302 | nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy)); | |
303 | nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy)); | |
304 | nGenSpec = (Int_t)(fZEMsp->Eval(zemenergy)); | |
305 | impPar = fZEMb->Eval(zemenergy); | |
306 | //printf(" Nspec ZEM = %f, Nspec ZDC = %f\n",fZEMsp->Eval(znenergy),fZDCPer->Eval(zdcenergy)); | |
307 | } | |
308 | // [1] ### Results in Chiara's PhD thesis -> 0<b<15 fm (Dec 2001) | |
309 | /*if(znenergy>158.5) nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy)); | |
310 | if(zpenergy>58.91) nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy)); | |
311 | if(zdcenergy>220.4) nGenSpec = (Int_t)(fZEMsp->Eval(zemenergy)); | |
312 | if(zdcenergy>225.) impPar = fZEMb->Eval(zemenergy);*/ | |
313 | // [2] ### Results from a new production -> 0<b<18 fm (Apr 2002) | |
314 | if(znenergy>162.) nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy)); | |
315 | if(zpenergy>59.75) nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy)); | |
316 | if(zdcenergy>221.5) nGenSpec = (Int_t)(fZEMsp->Eval(zemenergy)); | |
317 | if(zdcenergy>220.) impPar = fZEMb->Eval(zemenergy); | |
318 | ||
319 | if(nGenSpecN>125) nGenSpecN=125; | |
320 | else if(nGenSpecN<0) nGenSpecN=0; | |
321 | if(nGenSpecP>82) nGenSpecP=82; | |
322 | else if(nGenSpecP<0) nGenSpecP=0; | |
323 | if(nGenSpec>207) nGenSpec=207; | |
324 | else if(nGenSpec<0) nGenSpec=0; | |
325 | //printf(" NRecSpecN = %d, NRecSpecP = %d, NRecSpec = %d\n",nGenSpecN,nGenSpecP,nGenSpec); | |
326 | ||
327 | // --- Number of participants | |
328 | Int_t nPart, nPartTot; | |
329 | nPart = 207-nGenSpecN-nGenSpecP; | |
330 | nPartTot = 207-nGenSpec; | |
331 | //printf(" ### nPart(ZP+ZN) = %d, nPart(ZDC) = %d, b = %f fm\n",nPart,nPartTot,impPar); | |
332 | // if (GetDebug()) printf(" ### nPart = %d, b = %f fm\n",nPartTot,impPar); | |
333 | ||
334 | // create the output tree | |
335 | loader->MakeTree("R"); | |
336 | TTree* treeR = loader->TreeR(); | |
337 | AliZDCReco reco(znenergy, zpenergy, zdcenergy, zemenergy, | |
338 | nDetSpecN, nDetSpecP, nGenSpecN, nGenSpecP, nGenSpec, | |
339 | nPartTot, impPar); | |
340 | AliZDCReco* preco = &reco; | |
341 | const Int_t kBufferSize = 4000; | |
342 | treeR->Branch("ZDC", "AliZDCReco", &preco, kBufferSize); | |
343 | ||
344 | // write the output tree | |
345 | treeR->Fill(); | |
346 | loader->WriteRecPoints("OVERWRITE"); | |
347 | } | |
348 | ||
349 | //_____________________________________________________________________________ | |
350 | void AliZDCReconstructor::FillESD(AliRunLoader* runLoader, | |
351 | AliESD* esd) const | |
352 | { | |
353 | // fill energies and number of participants to the ESD | |
354 | ||
355 | AliLoader* loader = runLoader->GetLoader("ZDCLoader"); | |
356 | if (!loader) return; | |
357 | loader->LoadRecPoints(); | |
358 | ||
359 | TTree* treeR = loader->TreeR(); | |
360 | if (!treeR) return; | |
361 | AliZDCReco reco; | |
362 | AliZDCReco* preco = &reco; | |
363 | treeR->SetBranchAddress("ZDC", &preco); | |
364 | ||
365 | treeR->GetEntry(0); | |
366 | esd->SetZDC(reco.GetZNenergy(), reco.GetZPenergy(), reco.GetZEMenergy(), | |
367 | reco.GetNPart()); | |
368 | ||
369 | loader->UnloadRecPoints(); | |
370 | } |