1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // class for ZDC reconstruction //
22 ///////////////////////////////////////////////////////////////////////////////
27 #include "AliRunLoader.h"
28 #include "AliRawReader.h"
29 #include "AliESDEvent.h"
30 #include "AliZDCDigit.h"
31 #include "AliZDCRawStream.h"
32 #include "AliZDCReco.h"
33 #include "AliZDCReconstructor.h"
34 #include "AliZDCCalibData.h"
37 ClassImp(AliZDCReconstructor)
40 //_____________________________________________________________________________
41 AliZDCReconstructor:: AliZDCReconstructor() :
43 fZNCen(new TF1("fZNCen",
44 "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.)),
45 fZNPer(new TF1("fZNPer",
46 "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.)),
47 fZPCen(new TF1("fZPCen",
48 "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.)),
49 fZPPer(new TF1("fZPPer",
50 "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.)),
51 fZDCCen(new TF1("fZDCCen",
52 "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.)),
53 fZDCPer(new TF1("fZDCPer",
54 "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.)),
55 fbCen(new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.)),
56 fbPer(new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.)),
58 fZEMn(new TF1("fZEMn","121.7-0.1934*x+0.00007565*x*x",0.,1200.)),
59 fZEMp(new TF1("fZEMp","80.05-0.1315*x+0.00005327*x*x",0.,1200.)),
60 fZEMsp(new TF1("fZEMsp","201.7-0.325*x+0.0001292*x*x",0.,1200.)),
61 fZEMb(new TF1("fZEMb",
62 "13.83-0.02851*x+5.101e-5*x*x-7.305e-8*x*x*x+5.101e-11*x*x*x*x-1.25e-14*x*x*x*x*x",0.,1200.)),
64 fCalibData(GetCalibData())
67 // **** Default constructor
72 //_____________________________________________________________________________
73 AliZDCReconstructor::~AliZDCReconstructor()
93 //_____________________________________________________________________________
94 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
96 // *** Local ZDC reconstruction for digits
97 // Works on the current event
99 // Retrieving calibration data
101 for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
105 AliZDCDigit* pdigit = &digit;
106 digitsTree->SetBranchAddress("ZDC", &pdigit);
109 Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
110 Float_t dZEMCorrHG=0.;
111 Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
112 Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
113 Float_t dZEMCorrLG=0.;
114 Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
116 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
117 for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
118 digitsTree->GetEntry(iDigit);
119 if (!pdigit) continue;
122 Int_t det = digit.GetSector(0);
123 Int_t quad = digit.GetSector(1);
125 //printf("\n\t #%d det %d quad %d", iDigit, det, quad);
127 if(det == 1){ // *** ZN1
129 tZN1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
130 if(tZN1CorrHG[quad]<0.) tZN1CorrHG[quad] = 0.;
131 tZN1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
132 if(tZN1CorrLG[quad]<0.) tZN1CorrLG[quad] = 0.;
133 //printf("\t pedindex %d tZN1CorrHG[%d] = %1.0f tZN1CorrLG[%d] = %1.0f",
134 // pedindex, quad, tZN1CorrHG[quad], quad, tZN1CorrLG[quad]);
136 else if(det == 2){ // *** ZP1
138 tZP1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
139 if(tZP1CorrLG[quad]<0.) tZP1CorrLG[quad] = 0.;
140 tZP1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
141 if(tZP1CorrHG[quad]<0.) tZP1CorrHG[quad] = 0.;
142 //printf("\t pedindex %d tZP1CorrHG[%d] = %1.0f tZP1CorrLG[%d] = %1.0f",
143 // pedindex, quad, tZP1CorrHG[quad], quad, tZP1CorrLG[quad]);
146 if(quad == 1){ // *** ZEM1
148 dZEMCorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
149 if(dZEMCorrHG<0.) dZEMCorrHG = 0.;
150 dZEMCorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+2]);
151 if(dZEMCorrLG<0.) dZEMCorrLG = 0.;
152 //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n",
153 // pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEMCorrHG);
154 //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n",
155 // pedindex+2, digit.GetADCValue(1), meanPed[pedindex+2], dZEMCorrLG);
156 ////printf("\t pedindex %d dZEMCorrHG = %1.0f dZEMCorrLG = %1.0f\n", pedindex, dZEMCorrHG, dZEMCorrLG);
158 else if(quad == 2){ // *** ZEM1
160 dZEMCorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
161 if(dZEMCorrHG<0.) dZEMCorrHG = 0.;
162 dZEMCorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+2]);
163 if(dZEMCorrLG<0.) dZEMCorrLG = 0.;
164 //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n",
165 // pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEMCorrHG);
166 //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n",
167 // pedindex+2, digit.GetADCValue(1),meanPed[pedindex+2], dZEMCorrLG);
168 ////printf("\t pedindex %d dZEMCorrHG = %1.0f dZEMCorrLG = %1.0f\n", pedindex, dZEMCorrHG, dZEMCorrLG);
171 else if(det == 4){ // *** ZN2
173 tZN2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
174 if(tZN2CorrHG[quad]<0.) tZN2CorrHG[quad] = 0.;
175 tZN2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
176 if(tZN2CorrLG[quad]<0.) tZN2CorrLG[quad] = 0.;
177 //printf("\t pedindex %d tZN2CorrHG[%d] = %1.0f tZN2CorrLG[%d] = %1.0f\n",
178 // pedindex, quad, tZN2CorrHG[quad], quad, tZN2CorrLG[quad]);
180 else if(det == 5){ // *** ZP2
182 tZP2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
183 if(tZP2CorrHG[quad]<0.) tZP2CorrHG[quad] = 0.;
184 tZP2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
185 if(tZP2CorrLG[quad]<0.) tZP2CorrLG[quad] = 0.;
186 //printf("\t pedindex %d tZP2CorrHG[%d] = %1.0f tZP2CorrLG[%d] = %1.0f\n",
187 // pedindex, quad, tZP2CorrHG[quad], quad, tZP2CorrLG[quad]);
191 // reconstruct the event
192 ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG,
193 tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG,
194 tZP2CorrLG, dZEMCorrHG);
198 //_____________________________________________________________________________
199 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
201 // *** ZDC raw data reconstruction
202 // Works on the current event
204 // Retrieving calibration data
206 for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
210 // loop over raw data rawDatas
211 Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
212 Float_t dZEMCorrHG=0.;
213 Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
214 Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
215 Float_t dZEMCorrLG=0.;
216 Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
218 AliZDCRawStream rawData(rawReader);
219 while (rawData.Next()) {
220 if(rawData.IsADCDataWord()){
221 Int_t det = rawData.GetSector(0);
222 Int_t quad = rawData.GetSector(1);
223 Int_t gain = rawData.GetADCGain();
228 if(gain == 0) tZN1CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
229 else tZN1CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]);
233 if(gain == 0) tZP1CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
234 else tZP1CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]);
239 if(gain == 0) dZEMCorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
240 else dZEMCorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]);
243 pedindex = rawData.GetSector(1)+21;
244 if(gain == 0) dZEMCorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
245 else dZEMCorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]);
249 pedindex = rawData.GetSector(1)+24;
250 if(gain == 0) tZN2CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
251 else tZN2CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]);
254 pedindex = rawData.GetSector(1)+34;
255 if(gain == 0) tZP2CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
256 else tZP2CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]);
261 // reconstruct the event
262 ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG,
263 tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG,
264 tZP2CorrLG, dZEMCorrHG);
268 //_____________________________________________________________________________
269 void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree,
270 Float_t* ZN1ADCCorrHG, Float_t* ZP1ADCCorrHG,
271 Float_t* ZN2ADCCorrHG, Float_t* ZP2ADCCorrHG,
272 Float_t* ZN1ADCCorrLG, Float_t* ZP1ADCCorrLG,
273 Float_t* ZN2ADCCorrLG, Float_t* ZP2ADCCorrLG,
274 Float_t corrADCZEMHG) const
276 // ***** Reconstruct one event
278 // *** RECONSTRUCTION FROM SIMULATED DATA
279 // It passes trhough the no. of phe which is known from simulations
280 // --- ADCchannel -> photoelectrons
281 // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
282 // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
283 //Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
284 //zn1phe = ZN1Corr/convFactor;
285 //zp1phe = ZP1Corr/convFactor;
286 //zemphe = ZEMCorr/convFactor;
287 //zn2phe = ZN2Corr/convFactor;
288 //zp2phe = ZP2Corr/convFactor;
289 ////if AliDebug(1,Form("\n znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
291 //// --- Energy calibration
292 //// Conversion factors for hadronic ZDCs goes from phe yield to TRUE
293 //// incident energy (conversion from GeV to TeV is included); while for EM
294 //// calos conversion is from light yield to detected energy calculated by
295 //// GEANT NB -> ZN and ZP conversion factors are constant since incident
296 //// spectators have all the same energy, ZEM energy is obtained through a
297 //// fit over the whole range of incident particle energies
298 //// (obtained with full HIJING simulations)
299 //Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
300 //Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
301 //zn1energy = zn1phe/zn1phexTeV;
302 //zp1energy = zp1phe/zp1phexTeV;
303 //zdc1energy = zn1energy+zp1energy;
304 //zn2energy = zn2phe/zn2phexTeV;
305 //zp2energy = zp2phe/zp2phexTeV;
306 //zdc2energy = zn2energy+zp2energy;
307 //zemenergy = -4.81+0.3238*zemphe;
308 //if(zemenergy<0) zemenergy=0;
309 //// if AliDebug(1,Form(" znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
310 //// "\n zemenergy = %f TeV\n", znenergy, zpenergy,
311 //// zdcenergy, zemenergy);
312 //// if(zdcenergy==0)
313 //// if AliDebug(1,Form("\n\n ### ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
314 //// " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy);
317 // *** RECONSTRUCTION FROM "REAL" DATA
319 // Retrieving calibration data
320 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
321 for(Int_t ji=0; ji<5; ji++){
322 equalCoeffZN1[ji] = fCalibData->GetZN1EqualCoeff(ji);
323 equalCoeffZP1[ji] = fCalibData->GetZP1EqualCoeff(ji);
324 equalCoeffZN2[ji] = fCalibData->GetZN2EqualCoeff(ji);
325 equalCoeffZP2[ji] = fCalibData->GetZP2EqualCoeff(ji);
329 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fCalibData->GetEnCalib(ij);
331 Float_t endPointZEM = fCalibData->GetZEMEndValue();
332 Float_t cutFractionZEM = fCalibData->GetZEMCutFraction();
333 Float_t dZEMSup = fCalibData->GetDZEMSup();
334 Float_t dZEMInf = fCalibData->GetDZEMInf();
336 Float_t cutValueZEM = endPointZEM*cutFractionZEM;
337 Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
338 Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
340 Float_t maxValEZN1 = fCalibData->GetEZN1MaxValue();
341 Float_t maxValEZP1 = fCalibData->GetEZP1MaxValue();
342 Float_t maxValEZDC1 = fCalibData->GetEZDC1MaxValue();
343 Float_t maxValEZN2 = fCalibData->GetEZN2MaxValue();
344 Float_t maxValEZP2 = fCalibData->GetEZP2MaxValue();
345 Float_t maxValEZDC2 = fCalibData->GetEZDC2MaxValue();
347 //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
348 // " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
350 // Equalization of detector responses
351 Float_t equalTowZN1HG[5], equalTowZN2HG[5], equalTowZP1HG[5], equalTowZP2HG[5];
352 Float_t equalTowZN1LG[5], equalTowZN2LG[5], equalTowZP1LG[5], equalTowZP2LG[5];
353 for(Int_t gi=0; gi<5; gi++){
354 equalTowZN1HG[gi] = ZN1ADCCorrHG[gi]*equalCoeffZN1[gi];
355 equalTowZP1HG[gi] = ZP1ADCCorrHG[gi]*equalCoeffZP1[gi];
356 equalTowZN2HG[gi] = ZN2ADCCorrHG[gi]*equalCoeffZN2[gi];
357 equalTowZP2HG[gi] = ZP2ADCCorrHG[gi]*equalCoeffZP2[gi];
359 equalTowZN1LG[gi] = ZN1ADCCorrLG[gi]*equalCoeffZN1[gi];
360 equalTowZP1LG[gi] = ZP1ADCCorrLG[gi]*equalCoeffZP1[gi];
361 equalTowZN2LG[gi] = ZN2ADCCorrLG[gi]*equalCoeffZN2[gi];
362 equalTowZP2LG[gi] = ZP2ADCCorrLG[gi]*equalCoeffZP2[gi];
365 // Energy calibration of detector responses
366 Float_t calibTowZN1HG[5], calibTowZN2HG[5], calibTowZP1HG[5], calibTowZP2HG[5];
367 Float_t calibSumZN1HG=0., calibSumZN2HG=0., calibSumZP1HG=0., calibSumZP2HG=0.;
368 Float_t calibTowZN1LG[5], calibTowZN2LG[5], calibTowZP1LG[5], calibTowZP2LG[5];
369 Float_t calibSumZN1LG=0., calibSumZN2LG=0., calibSumZ12LG=0., calibSumZP2LG=0.;
370 for(Int_t gi=0; gi<5; gi++){
371 calibTowZN1HG[gi] = equalTowZN1HG[gi]*calibEne[0];
372 calibTowZP1HG[gi] = equalTowZP1HG[gi]*calibEne[1];
373 calibTowZN2HG[gi] = equalTowZN2HG[gi]*calibEne[2];
374 calibTowZP2HG[gi] = equalTowZP2HG[gi]*calibEne[3];
375 calibSumZN1HG += calibTowZN1HG[gi];
376 calibSumZP1HG += calibTowZP1HG[gi];
377 calibSumZN2HG += calibTowZN2HG[gi];
378 calibSumZP2HG += calibTowZP2HG[gi];
380 calibTowZN1LG[gi] = equalTowZN1LG[gi]*calibEne[0];
381 calibTowZP1LG[gi] = equalTowZP1LG[gi]*calibEne[1];
382 calibTowZN2LG[gi] = equalTowZN2LG[gi]*calibEne[2];
383 calibTowZP2LG[gi] = equalTowZP2LG[gi]*calibEne[3];
384 calibSumZN1LG += calibTowZN1LG[gi];
385 calibSumZ12LG += calibTowZP1LG[gi];
386 calibSumZN2LG += calibTowZN2LG[gi];
387 calibSumZP2LG += calibTowZP2LG[gi];
390 // --- Number of detected spectator nucleons
391 // *** N.B. -> It works only in Pb-Pb
392 Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
393 nDetSpecNLeft = (Int_t) (calibSumZN1HG/2.760);
394 nDetSpecPLeft = (Int_t) (calibSumZP1HG/2.760);
395 nDetSpecNRight = (Int_t) (calibSumZN2HG/2.760);
396 nDetSpecPRight = (Int_t) (calibSumZP2HG/2.760);
397 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
398 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
399 nDetSpecNRight, nDetSpecPRight);*/
401 // --- Number of generated spectator nucleons (from HIJING parameterization)
402 Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
403 Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
406 // *** RECONSTRUCTION FROM SIMULATED DATA
407 // Cut value for Ezem (GeV)
408 // ### Results from production -> 0<b<18 fm (Apr 2002)
409 /*Float_t eZEMCut = 420.;
410 Float_t deltaEZEMSup = 690.;
411 Float_t deltaEZEMInf = 270.;
412 if(zemenergy > (eZEMCut+deltaEZEMSup)){
413 nGenSpecNLeft = (Int_t) (fZNCen->Eval(ZN1CalibSum));
414 nGenSpecPLeft = (Int_t) (fZPCen->Eval(ZP1CalibSum));
415 nGenSpecLeft = (Int_t) (fZDCCen->Eval(ZN1CalibSum+ZP1CalibSum));
416 nGenSpecNRight = (Int_t) (fZNCen->Eval(ZN2CalibSum));
417 nGenSpecPRight = (Int_t) (fZNCen->Eval(ZP2CalibSum));
418 nGenSpecRight = (Int_t) (fZNCen->Eval(ZN2CalibSum+ZP2CalibSum));
419 impPar = fbCen->Eval(ZN1CalibSum+ZP1CalibSum);
421 else if(zemenergy < (eZEMCut-deltaEZEMInf)){
422 nGenSpecNLeft = (Int_t) (fZNPer->Eval(ZN1CalibSum));
423 nGenSpecPLeft = (Int_t) (fZPPer->Eval(ZP1CalibSum));
424 nGenSpecLeft = (Int_t) (fZDCPer->Eval(ZN1CalibSum+ZP1CalibSum));
425 impPar = fbPer->Eval(ZN1CalibSum+ZP1CalibSum);
427 else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
428 nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
429 nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
430 nGenSpecLeft = (Int_t)(fZEMsp->Eval(zemenergy));
431 impPar = fZEMb->Eval(zemenergy);
433 // ### Results from production -> 0<b<18 fm (Apr 2002)
434 if(ZN1CalibSum>162.) nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
435 if(ZP1CalibSum>59.75) nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
436 if(ZN1CalibSum+ZP1CalibSum>221.5) nGenSpecLeft = (Int_t)(fZEMsp->Eval(zemenergy));
437 if(ZN1CalibSum+ZP1CalibSum>220.) impPar = fZEMb->Eval(zemenergy);
441 // *** RECONSTRUCTION FROM REAL DATA
443 if(corrADCZEMHG > supValueZEM){
444 nGenSpecNLeft = (Int_t) (fZNCen->Eval(calibSumZN1HG));
445 nGenSpecPLeft = (Int_t) (fZPCen->Eval(calibSumZP1HG));
446 nGenSpecLeft = (Int_t) (fZDCCen->Eval(calibSumZN1HG+calibSumZP1HG));
447 nGenSpecNRight = (Int_t) (fZNCen->Eval(calibSumZN2HG));
448 nGenSpecPRight = (Int_t) (fZNCen->Eval(calibSumZP2HG));
449 nGenSpecRight = (Int_t) (fZNCen->Eval(calibSumZN2HG+calibSumZP2HG));
450 impPar = fbCen->Eval(calibSumZN1HG+calibSumZP1HG);
452 else if(corrADCZEMHG < infValueZEM){
453 nGenSpecNLeft = (Int_t) (fZNPer->Eval(calibSumZN1HG));
454 nGenSpecPLeft = (Int_t) (fZPPer->Eval(calibSumZP1HG));
455 nGenSpecLeft = (Int_t) (fZDCPer->Eval(calibSumZN1HG+calibSumZP1HG));
456 impPar = fbPer->Eval(calibSumZN1HG+calibSumZP1HG);
458 else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
459 nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
460 nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
461 nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
462 impPar = fZEMb->Eval(corrADCZEMHG);
465 if(calibSumZN1HG/maxValEZN1>1.) nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
466 if(calibSumZP1HG/maxValEZP1>1.) nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
467 if((calibSumZN1HG+calibSumZP1HG/maxValEZDC1)>1.){
468 nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
469 impPar = fZEMb->Eval(corrADCZEMHG);
471 if(calibSumZN2HG/maxValEZN2>1.) nGenSpecNRight = (Int_t) (fZEMn->Eval(corrADCZEMHG));
472 if(calibSumZP2HG/maxValEZP2>1.) nGenSpecPRight = (Int_t) (fZEMp->Eval(corrADCZEMHG));
473 if((calibSumZN2HG+calibSumZP2HG/maxValEZDC2)>1.) nGenSpecRight = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
475 if(nGenSpecNLeft>125) nGenSpecNLeft=125;
476 else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
477 if(nGenSpecPLeft>82) nGenSpecPLeft=82;
478 else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
479 if(nGenSpecLeft>207) nGenSpecLeft=207;
480 else if(nGenSpecLeft<0) nGenSpecLeft=0;
482 // --- Number of generated participants (from HIJING parameterization)
483 Int_t nPart, nPartTotLeft, nPartTotRight;
484 nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
485 nPartTotLeft = 207-nGenSpecLeft;
486 nPartTotRight = 207-nGenSpecRight;
488 if(nPartTotLeft<0) nPartTotLeft=0;
489 if(nPartTotRight<0) nPartTotRight=0;
492 printf("\n\t AliZDCReconstructor -> calibSumZN1HG %1.0f, calibSumZP1HG %1.0f,"
493 " calibSumZN2HG %1.0f, calibSumZP2HG %1.0f, corrADCZEMHG %1.0f\n",
494 calibSumZN1HG,calibSumZP1HG,calibSumZN2HG,calibSumZP2HG,corrADCZEMHG);
495 printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
496 "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n",
497 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft,
498 nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
499 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
501 // create the output tree
502 AliZDCReco reco(calibSumZN1HG, calibSumZP1HG, calibSumZN2HG, calibSumZP2HG,
503 calibTowZN1LG, calibTowZN2LG, calibTowZP1LG, calibTowZP2LG,
505 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
506 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
507 nGenSpecPRight, nGenSpecRight,
508 nPartTotLeft, nPartTotRight, impPar);
510 AliZDCReco* preco = &reco;
511 const Int_t kBufferSize = 4000;
512 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
514 // write the output tree
515 clustersTree->Fill();
518 //_____________________________________________________________________________
519 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
521 // fill energies and number of participants to the ESD
524 AliZDCReco* preco = &reco;
525 clustersTree->SetBranchAddress("ZDC", &preco);
527 clustersTree->GetEntry(0);
528 /*Double_t tZN1Ene[4], tZN2Ene[4];
529 for(Int_t i=0; i<4; i++){
530 tZN1Ene[i] = reco.GetZN1EnTow(i);
531 tZN2Ene[i] = reco.GetZN2EnTow(i);
533 esd->SetZDC(tZN1Ene, tZN2Ene, reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEMsignal(),
534 reco.GetZN2Energy(), reco.GetZP2Energy(),
535 reco.GetNPartLeft());
537 esd->SetZDC(reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEMsignal(),
538 reco.GetZN2Energy(), reco.GetZP2Energy(),
539 reco.GetNPartLeft());
541 /*Double_t tZN1Ene[4], tZN2Ene[4];
542 for(Int_t i=0; i<4; i++){
543 tZN1Ene[i] = reco.GetZN1EnTow(i);
544 tZN2Ene[i] = reco.GetZN2EnTow(i);
546 //esd->SetZN1TowerEnergy(tZN1Ene);
547 //esd->SetZN2TowerEnergy(tZN2Ene);
551 //_____________________________________________________________________________
552 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
554 // Setting the storage
556 Bool_t deleteManager = kFALSE;
558 AliCDBManager *manager = AliCDBManager::Instance();
559 AliCDBStorage *defstorage = manager->GetDefaultStorage();
561 if(!defstorage || !(defstorage->Contains("ZDC"))){
562 AliWarning("No default storage set or default storage doesn't contain ZDC!");
563 manager->SetDefaultStorage(uri);
564 deleteManager = kTRUE;
567 AliCDBStorage *storage = manager->GetDefaultStorage();
570 AliCDBManager::Instance()->UnsetDefaultStorage();
571 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
577 //_____________________________________________________________________________
578 AliZDCCalibData* AliZDCReconstructor::GetCalibData() const
581 // Getting calibration object for ZDC set
583 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Data");
584 if(!entry) AliFatal("No calibration data loaded!");
586 AliZDCCalibData *calibdata = dynamic_cast<AliZDCCalibData*> (entry->GetObject());
587 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
591 /**************************************************************************
592 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
594 * Author: The ALICE Off-line Project. *
595 * Contributors are mentioned in the code where appropriate. *
597 * Permission to use, copy, modify and distribute this software and its *
598 * documentation strictly for non-commercial purposes is hereby granted *
599 * without fee, provided that the above copyright notice appears in all *
600 * copies and that both the copyright notice and this permission notice *
601 * appear in the supporting documentation. The authors make no claims *
602 * about the suitability of this software for any purpose. It is *
603 * provided "as is" without express or implied warranty. *
604 **************************************************************************/
608 ///////////////////////////////////////////////////////////////////////////////
610 // class for ZDC reconstruction //
612 ///////////////////////////////////////////////////////////////////////////////
617 #include "AliRunLoader.h"
618 #include "AliRawReader.h"
619 #include "AliESDEvent.h"
620 #include "AliZDCDigit.h"
621 #include "AliZDCRawStream.h"
622 #include "AliZDCReco.h"
623 #include "AliZDCReconstructor.h"
624 #include "AliZDCCalibData.h"
627 ClassImp(AliZDCReconstructor)
630 //_____________________________________________________________________________
631 AliZDCReconstructor:: AliZDCReconstructor() :
633 fZNCen(new TF1("fZNCen",
634 "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.)),
635 fZNPer(new TF1("fZNPer",
636 "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.)),
637 fZPCen(new TF1("fZPCen",
638 "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.)),
639 fZPPer(new TF1("fZPPer",
640 "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.)),
641 fZDCCen(new TF1("fZDCCen",
642 "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.)),
643 fZDCPer(new TF1("fZDCPer",
644 "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.)),
645 fbCen(new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.)),
646 fbPer(new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.)),
648 fZEMn(new TF1("fZEMn","121.7-0.1934*x+0.00007565*x*x",0.,1200.)),
649 fZEMp(new TF1("fZEMp","80.05-0.1315*x+0.00005327*x*x",0.,1200.)),
650 fZEMsp(new TF1("fZEMsp","201.7-0.325*x+0.0001292*x*x",0.,1200.)),
651 fZEMb(new TF1("fZEMb",
652 "13.83-0.02851*x+5.101e-5*x*x-7.305e-8*x*x*x+5.101e-11*x*x*x*x-1.25e-14*x*x*x*x*x",0.,1200.)),
654 fCalibData(GetCalibData())
657 // **** Default constructor
662 //_____________________________________________________________________________
663 AliZDCReconstructor::~AliZDCReconstructor()
683 //_____________________________________________________________________________
684 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
686 // *** Local ZDC reconstruction for digits
687 // Works on the current event
689 // Retrieving calibration data
691 for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
695 AliZDCDigit* pdigit = &digit;
696 digitsTree->SetBranchAddress("ZDC", &pdigit);
699 Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
700 Float_t dZEMCorrHG=0.;
701 Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
702 Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
703 Float_t dZEMCorrLG=0.;
704 Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
706 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
707 for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
708 digitsTree->GetEntry(iDigit);
709 if (!pdigit) continue;
712 Int_t det = digit.GetSector(0);
713 Int_t quad = digit.GetSector(1);
715 //printf("\n\t #%d det %d quad %d", iDigit, det, quad);
717 if(det == 1){ // *** ZN1
719 tZN1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
720 if(tZN1CorrHG[quad]<0.) tZN1CorrHG[quad] = 0.;
721 tZN1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
722 if(tZN1CorrLG[quad]<0.) tZN1CorrLG[quad] = 0.;
723 //printf("\t pedindex %d tZN1CorrHG[%d] = %1.0f tZN1CorrLG[%d] = %1.0f",
724 // pedindex, quad, tZN1CorrHG[quad], quad, tZN1CorrLG[quad]);
726 else if(det == 2){ // *** ZP1
728 tZP1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
729 if(tZP1CorrLG[quad]<0.) tZP1CorrLG[quad] = 0.;
730 tZP1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
731 if(tZP1CorrHG[quad]<0.) tZP1CorrHG[quad] = 0.;
732 //printf("\t pedindex %d tZP1CorrHG[%d] = %1.0f tZP1CorrLG[%d] = %1.0f",
733 // pedindex, quad, tZP1CorrHG[quad], quad, tZP1CorrLG[quad]);
736 if(quad == 1){ // *** ZEM1
738 dZEMCorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
739 if(dZEMCorrHG<0.) dZEMCorrHG = 0.;
740 dZEMCorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+2]);
741 if(dZEMCorrLG<0.) dZEMCorrLG = 0.;
742 //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n",
743 // pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEMCorrHG);
744 //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n",
745 // pedindex+2, digit.GetADCValue(1), meanPed[pedindex+2], dZEMCorrLG);
746 ////printf("\t pedindex %d dZEMCorrHG = %1.0f dZEMCorrLG = %1.0f\n", pedindex, dZEMCorrHG, dZEMCorrLG);
748 else if(quad == 2){ // *** ZEM1
750 dZEMCorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
751 if(dZEMCorrHG<0.) dZEMCorrHG = 0.;
752 dZEMCorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+2]);
753 if(dZEMCorrLG<0.) dZEMCorrLG = 0.;
754 //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n",
755 // pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEMCorrHG);
756 //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n",
757 // pedindex+2, digit.GetADCValue(1),meanPed[pedindex+2], dZEMCorrLG);
758 ////printf("\t pedindex %d dZEMCorrHG = %1.0f dZEMCorrLG = %1.0f\n", pedindex, dZEMCorrHG, dZEMCorrLG);
761 else if(det == 4){ // *** ZN2
763 tZN2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
764 if(tZN2CorrHG[quad]<0.) tZN2CorrHG[quad] = 0.;
765 tZN2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
766 if(tZN2CorrLG[quad]<0.) tZN2CorrLG[quad] = 0.;
767 //printf("\t pedindex %d tZN2CorrHG[%d] = %1.0f tZN2CorrLG[%d] = %1.0f\n",
768 // pedindex, quad, tZN2CorrHG[quad], quad, tZN2CorrLG[quad]);
770 else if(det == 5){ // *** ZP2
772 tZP2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
773 if(tZP2CorrHG[quad]<0.) tZP2CorrHG[quad] = 0.;
774 tZP2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
775 if(tZP2CorrLG[quad]<0.) tZP2CorrLG[quad] = 0.;
776 //printf("\t pedindex %d tZP2CorrHG[%d] = %1.0f tZP2CorrLG[%d] = %1.0f\n",
777 // pedindex, quad, tZP2CorrHG[quad], quad, tZP2CorrLG[quad]);
781 // reconstruct the event
782 ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG,
783 tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG,
784 tZP2CorrLG, dZEMCorrHG);
788 //_____________________________________________________________________________
789 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
791 // *** ZDC raw data reconstruction
792 // Works on the current event
794 // Retrieving calibration data
796 for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
800 // loop over raw data rawDatas
801 Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
802 Float_t dZEMCorrHG=0.;
803 Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
804 Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
805 Float_t dZEMCorrLG=0.;
806 Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
808 AliZDCRawStream rawData(rawReader);
809 while (rawData.Next()) {
810 if(rawData.IsADCDataWord()){
811 Int_t det = rawData.GetSector(0);
812 Int_t quad = rawData.GetSector(1);
813 Int_t gain = rawData.GetADCGain();
818 if(gain == 0) tZN1CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
819 else tZN1CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]);
823 if(gain == 0) tZP1CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
824 else tZP1CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]);
829 if(gain == 0) dZEMCorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
830 else dZEMCorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]);
833 pedindex = rawData.GetSector(1)+21;
834 if(gain == 0) dZEMCorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
835 else dZEMCorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]);
839 pedindex = rawData.GetSector(1)+24;
840 if(gain == 0) tZN2CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
841 else tZN2CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]);
844 pedindex = rawData.GetSector(1)+34;
845 if(gain == 0) tZP2CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
846 else tZP2CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]);
851 // reconstruct the event
852 ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG,
853 tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG,
854 tZP2CorrLG, dZEMCorrHG);
858 //_____________________________________________________________________________
859 void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree,
860 Float_t* ZN1ADCCorrHG, Float_t* ZP1ADCCorrHG,
861 Float_t* ZN2ADCCorrHG, Float_t* ZP2ADCCorrHG,
862 Float_t* ZN1ADCCorrLG, Float_t* ZP1ADCCorrLG,
863 Float_t* ZN2ADCCorrLG, Float_t* ZP2ADCCorrLG,
864 Float_t corrADCZEMHG) const
866 // ***** Reconstruct one event
868 // *** RECONSTRUCTION FROM SIMULATED DATA
869 // It passes trhough the no. of phe which is known from simulations
870 // --- ADCchannel -> photoelectrons
871 // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
872 // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
873 //Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
874 //zn1phe = ZN1Corr/convFactor;
875 //zp1phe = ZP1Corr/convFactor;
876 //zemphe = ZEMCorr/convFactor;
877 //zn2phe = ZN2Corr/convFactor;
878 //zp2phe = ZP2Corr/convFactor;
879 ////if AliDebug(1,Form("\n znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
881 //// --- Energy calibration
882 //// Conversion factors for hadronic ZDCs goes from phe yield to TRUE
883 //// incident energy (conversion from GeV to TeV is included); while for EM
884 //// calos conversion is from light yield to detected energy calculated by
885 //// GEANT NB -> ZN and ZP conversion factors are constant since incident
886 //// spectators have all the same energy, ZEM energy is obtained through a
887 //// fit over the whole range of incident particle energies
888 //// (obtained with full HIJING simulations)
889 //Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
890 //Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
891 //zn1energy = zn1phe/zn1phexTeV;
892 //zp1energy = zp1phe/zp1phexTeV;
893 //zdc1energy = zn1energy+zp1energy;
894 //zn2energy = zn2phe/zn2phexTeV;
895 //zp2energy = zp2phe/zp2phexTeV;
896 //zdc2energy = zn2energy+zp2energy;
897 //zemenergy = -4.81+0.3238*zemphe;
898 //if(zemenergy<0) zemenergy=0;
899 //// if AliDebug(1,Form(" znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
900 //// "\n zemenergy = %f TeV\n", znenergy, zpenergy,
901 //// zdcenergy, zemenergy);
902 //// if(zdcenergy==0)
903 //// if AliDebug(1,Form("\n\n ### ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
904 //// " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy);
907 // *** RECONSTRUCTION FROM "REAL" DATA
909 // Retrieving calibration data
910 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
911 for(Int_t ji=0; ji<5; ji++){
912 equalCoeffZN1[ji] = fCalibData->GetZN1EqualCoeff(ji);
913 equalCoeffZP1[ji] = fCalibData->GetZP1EqualCoeff(ji);
914 equalCoeffZN2[ji] = fCalibData->GetZN2EqualCoeff(ji);
915 equalCoeffZP2[ji] = fCalibData->GetZP2EqualCoeff(ji);
919 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fCalibData->GetEnCalib(ij);
921 Float_t endPointZEM = fCalibData->GetZEMEndValue();
922 Float_t cutFractionZEM = fCalibData->GetZEMCutFraction();
923 Float_t dZEMSup = fCalibData->GetDZEMSup();
924 Float_t dZEMInf = fCalibData->GetDZEMInf();
926 Float_t cutValueZEM = endPointZEM*cutFractionZEM;
927 Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
928 Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
930 Float_t maxValEZN1 = fCalibData->GetEZN1MaxValue();
931 Float_t maxValEZP1 = fCalibData->GetEZP1MaxValue();
932 Float_t maxValEZDC1 = fCalibData->GetEZDC1MaxValue();
933 Float_t maxValEZN2 = fCalibData->GetEZN2MaxValue();
934 Float_t maxValEZP2 = fCalibData->GetEZP2MaxValue();
935 Float_t maxValEZDC2 = fCalibData->GetEZDC2MaxValue();
937 //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
938 // " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
940 // Equalization of detector responses
941 Float_t equalTowZN1HG[5], equalTowZN2HG[5], equalTowZP1HG[5], equalTowZP2HG[5];
942 Float_t equalTowZN1LG[5], equalTowZN2LG[5], equalTowZP1LG[5], equalTowZP2LG[5];
943 for(Int_t gi=0; gi<5; gi++){
944 equalTowZN1HG[gi] = ZN1ADCCorrHG[gi]*equalCoeffZN1[gi];
945 equalTowZP1HG[gi] = ZP1ADCCorrHG[gi]*equalCoeffZP1[gi];
946 equalTowZN2HG[gi] = ZN2ADCCorrHG[gi]*equalCoeffZN2[gi];
947 equalTowZP2HG[gi] = ZP2ADCCorrHG[gi]*equalCoeffZP2[gi];
949 equalTowZN1LG[gi] = ZN1ADCCorrLG[gi]*equalCoeffZN1[gi];
950 equalTowZP1LG[gi] = ZP1ADCCorrLG[gi]*equalCoeffZP1[gi];
951 equalTowZN2LG[gi] = ZN2ADCCorrLG[gi]*equalCoeffZN2[gi];
952 equalTowZP2LG[gi] = ZP2ADCCorrLG[gi]*equalCoeffZP2[gi];
955 // Energy calibration of detector responses
956 Float_t calibTowZN1HG[5], calibTowZN2HG[5], calibTowZP1HG[5], calibTowZP2HG[5];
957 Float_t calibSumZN1HG=0., calibSumZN2HG=0., calibSumZP1HG=0., calibSumZP2HG=0.;
958 Float_t calibTowZN1LG[5], calibTowZN2LG[5], calibTowZP1LG[5], calibTowZP2LG[5];
959 Float_t calibSumZN1LG=0., calibSumZN2LG=0., calibSumZ12LG=0., calibSumZP2LG=0.;
960 for(Int_t gi=0; gi<5; gi++){
961 calibTowZN1HG[gi] = equalTowZN1HG[gi]*calibEne[0];
962 calibTowZP1HG[gi] = equalTowZP1HG[gi]*calibEne[1];
963 calibTowZN2HG[gi] = equalTowZN2HG[gi]*calibEne[2];
964 calibTowZP2HG[gi] = equalTowZP2HG[gi]*calibEne[3];
965 calibSumZN1HG += calibTowZN1HG[gi];
966 calibSumZP1HG += calibTowZP1HG[gi];
967 calibSumZN2HG += calibTowZN2HG[gi];
968 calibSumZP2HG += calibTowZP2HG[gi];
970 calibTowZN1LG[gi] = equalTowZN1LG[gi]*calibEne[0];
971 calibTowZP1LG[gi] = equalTowZP1LG[gi]*calibEne[1];
972 calibTowZN2LG[gi] = equalTowZN2LG[gi]*calibEne[2];
973 calibTowZP2LG[gi] = equalTowZP2LG[gi]*calibEne[3];
974 calibSumZN1LG += calibTowZN1LG[gi];
975 calibSumZ12LG += calibTowZP1LG[gi];
976 calibSumZN2LG += calibTowZN2LG[gi];
977 calibSumZP2LG += calibTowZP2LG[gi];
980 // --- Number of detected spectator nucleons
981 // *** N.B. -> It works only in Pb-Pb
982 Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
983 nDetSpecNLeft = (Int_t) (calibSumZN1HG/2.760);
984 nDetSpecPLeft = (Int_t) (calibSumZP1HG/2.760);
985 nDetSpecNRight = (Int_t) (calibSumZN2HG/2.760);
986 nDetSpecPRight = (Int_t) (calibSumZP2HG/2.760);
987 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
988 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
989 nDetSpecNRight, nDetSpecPRight);*/
991 // --- Number of generated spectator nucleons (from HIJING parameterization)
992 Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
993 Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
996 // *** RECONSTRUCTION FROM SIMULATED DATA
997 // Cut value for Ezem (GeV)
998 // ### Results from production -> 0<b<18 fm (Apr 2002)
999 /*Float_t eZEMCut = 420.;
1000 Float_t deltaEZEMSup = 690.;
1001 Float_t deltaEZEMInf = 270.;
1002 if(zemenergy > (eZEMCut+deltaEZEMSup)){
1003 nGenSpecNLeft = (Int_t) (fZNCen->Eval(ZN1CalibSum));
1004 nGenSpecPLeft = (Int_t) (fZPCen->Eval(ZP1CalibSum));
1005 nGenSpecLeft = (Int_t) (fZDCCen->Eval(ZN1CalibSum+ZP1CalibSum));
1006 nGenSpecNRight = (Int_t) (fZNCen->Eval(ZN2CalibSum));
1007 nGenSpecPRight = (Int_t) (fZNCen->Eval(ZP2CalibSum));
1008 nGenSpecRight = (Int_t) (fZNCen->Eval(ZN2CalibSum+ZP2CalibSum));
1009 impPar = fbCen->Eval(ZN1CalibSum+ZP1CalibSum);
1011 else if(zemenergy < (eZEMCut-deltaEZEMInf)){
1012 nGenSpecNLeft = (Int_t) (fZNPer->Eval(ZN1CalibSum));
1013 nGenSpecPLeft = (Int_t) (fZPPer->Eval(ZP1CalibSum));
1014 nGenSpecLeft = (Int_t) (fZDCPer->Eval(ZN1CalibSum+ZP1CalibSum));
1015 impPar = fbPer->Eval(ZN1CalibSum+ZP1CalibSum);
1017 else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
1018 nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
1019 nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
1020 nGenSpecLeft = (Int_t)(fZEMsp->Eval(zemenergy));
1021 impPar = fZEMb->Eval(zemenergy);
1023 // ### Results from production -> 0<b<18 fm (Apr 2002)
1024 if(ZN1CalibSum>162.) nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
1025 if(ZP1CalibSum>59.75) nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
1026 if(ZN1CalibSum+ZP1CalibSum>221.5) nGenSpecLeft = (Int_t)(fZEMsp->Eval(zemenergy));
1027 if(ZN1CalibSum+ZP1CalibSum>220.) impPar = fZEMb->Eval(zemenergy);
1031 // *** RECONSTRUCTION FROM REAL DATA
1033 if(corrADCZEMHG > supValueZEM){
1034 nGenSpecNLeft = (Int_t) (fZNCen->Eval(calibSumZN1HG));
1035 nGenSpecPLeft = (Int_t) (fZPCen->Eval(calibSumZP1HG));
1036 nGenSpecLeft = (Int_t) (fZDCCen->Eval(calibSumZN1HG+calibSumZP1HG));
1037 nGenSpecNRight = (Int_t) (fZNCen->Eval(calibSumZN2HG));
1038 nGenSpecPRight = (Int_t) (fZNCen->Eval(calibSumZP2HG));
1039 nGenSpecRight = (Int_t) (fZNCen->Eval(calibSumZN2HG+calibSumZP2HG));
1040 impPar = fbCen->Eval(calibSumZN1HG+calibSumZP1HG);
1042 else if(corrADCZEMHG < infValueZEM){
1043 nGenSpecNLeft = (Int_t) (fZNPer->Eval(calibSumZN1HG));
1044 nGenSpecPLeft = (Int_t) (fZPPer->Eval(calibSumZP1HG));
1045 nGenSpecLeft = (Int_t) (fZDCPer->Eval(calibSumZN1HG+calibSumZP1HG));
1046 impPar = fbPer->Eval(calibSumZN1HG+calibSumZP1HG);
1048 else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
1049 nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
1050 nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
1051 nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
1052 impPar = fZEMb->Eval(corrADCZEMHG);
1055 if(calibSumZN1HG/maxValEZN1>1.) nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
1056 if(calibSumZP1HG/maxValEZP1>1.) nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
1057 if((calibSumZN1HG+calibSumZP1HG/maxValEZDC1)>1.){
1058 nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
1059 impPar = fZEMb->Eval(corrADCZEMHG);
1061 if(calibSumZN2HG/maxValEZN2>1.) nGenSpecNRight = (Int_t) (fZEMn->Eval(corrADCZEMHG));
1062 if(calibSumZP2HG/maxValEZP2>1.) nGenSpecPRight = (Int_t) (fZEMp->Eval(corrADCZEMHG));
1063 if((calibSumZN2HG+calibSumZP2HG/maxValEZDC2)>1.) nGenSpecRight = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
1065 if(nGenSpecNLeft>125) nGenSpecNLeft=125;
1066 else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
1067 if(nGenSpecPLeft>82) nGenSpecPLeft=82;
1068 else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
1069 if(nGenSpecLeft>207) nGenSpecLeft=207;
1070 else if(nGenSpecLeft<0) nGenSpecLeft=0;
1072 // --- Number of generated participants (from HIJING parameterization)
1073 Int_t nPart, nPartTotLeft, nPartTotRight;
1074 nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
1075 nPartTotLeft = 207-nGenSpecLeft;
1076 nPartTotRight = 207-nGenSpecRight;
1077 if(nPart<0) nPart=0;
1078 if(nPartTotLeft<0) nPartTotLeft=0;
1079 if(nPartTotRight<0) nPartTotRight=0;
1082 printf("\n\t AliZDCReconstructor -> calibSumZN1HG %1.0f, calibSumZP1HG %1.0f,"
1083 " calibSumZN2HG %1.0f, calibSumZP2HG %1.0f, corrADCZEMHG %1.0f\n",
1084 calibSumZN1HG,calibSumZP1HG,calibSumZN2HG,calibSumZP2HG,corrADCZEMHG);
1085 printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
1086 "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n",
1087 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft,
1088 nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
1089 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1091 // create the output tree
1092 AliZDCReco reco(calibSumZN1HG, calibSumZP1HG, calibSumZN2HG, calibSumZP2HG,
1093 calibTowZN1LG, calibTowZN2LG, calibTowZP1LG, calibTowZP2LG,
1095 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1096 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
1097 nGenSpecPRight, nGenSpecRight,
1098 nPartTotLeft, nPartTotRight, impPar);
1100 AliZDCReco* preco = &reco;
1101 const Int_t kBufferSize = 4000;
1102 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
1104 // write the output tree
1105 clustersTree->Fill();
1108 //_____________________________________________________________________________
1109 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1111 // fill energies and number of participants to the ESD
1114 AliZDCReco* preco = &reco;
1115 clustersTree->SetBranchAddress("ZDC", &preco);
1117 clustersTree->GetEntry(0);
1118 /*Double_t tZN1Ene[4], tZN2Ene[4];
1119 for(Int_t i=0; i<4; i++){
1120 tZN1Ene[i] = reco.GetZN1EnTow(i);
1121 tZN2Ene[i] = reco.GetZN2EnTow(i);
1123 esd->SetZDC(tZN1Ene, tZN2Ene, reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEMsignal(),
1124 reco.GetZN2Energy(), reco.GetZP2Energy(),
1125 reco.GetNPartLeft());
1127 esd->SetZDC(reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEMsignal(),
1128 reco.GetZN2Energy(), reco.GetZP2Energy(),
1129 reco.GetNPartLeft());
1131 /*Double_t tZN1Ene[4], tZN2Ene[4];
1132 for(Int_t i=0; i<4; i++){
1133 tZN1Ene[i] = reco.GetZN1EnTow(i);
1134 tZN2Ene[i] = reco.GetZN2EnTow(i);
1136 //esd->SetZN1TowerEnergy(tZN1Ene);
1137 //esd->SetZN2TowerEnergy(tZN2Ene);
1141 //_____________________________________________________________________________
1142 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1144 // Setting the storage
1146 Bool_t deleteManager = kFALSE;
1148 AliCDBManager *manager = AliCDBManager::Instance();
1149 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1151 if(!defstorage || !(defstorage->Contains("ZDC"))){
1152 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1153 manager->SetDefaultStorage(uri);
1154 deleteManager = kTRUE;
1157 AliCDBStorage *storage = manager->GetDefaultStorage();
1160 AliCDBManager::Instance()->UnsetDefaultStorage();
1161 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1167 //_____________________________________________________________________________
1168 AliZDCCalibData* AliZDCReconstructor::GetCalibData() const
1171 // Getting calibration object for ZDC set
1173 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Data");
1174 if(!entry) AliFatal("No calibration data loaded!");
1176 AliZDCCalibData *calibdata = dynamic_cast<AliZDCCalibData*> (entry->GetObject());
1177 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");