AliCutTask uses global cut initialization
[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx
CommitLineData
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 *
f5d41205 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 "AliESDEvent.h"
a85132e7 30#include "AliESDZDC.h"
f5d41205 31#include "AliZDCDigit.h"
32#include "AliZDCRawStream.h"
33#include "AliZDCReco.h"
34#include "AliZDCReconstructor.h"
6024ec85 35#include "AliZDCPedestals.h"
36#include "AliZDCCalib.h"
37#include "AliZDCRecParam.h"
f5d41205 38
39
40ClassImp(AliZDCReconstructor)
41
42
43//_____________________________________________________________________________
44AliZDCReconstructor:: AliZDCReconstructor() :
45
46 fZNCen(new TF1("fZNCen",
47 "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.)),
48 fZNPer(new TF1("fZNPer",
49 "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.)),
50 fZPCen(new TF1("fZPCen",
51 "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.)),
52 fZPPer(new TF1("fZPPer",
53 "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.)),
54 fZDCCen(new TF1("fZDCCen",
55 "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.)),
56 fZDCPer(new TF1("fZDCPer",
57 "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.)),
58 fbCen(new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.)),
59 fbPer(new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.)),
60 //
61 fZEMn(new TF1("fZEMn","121.7-0.1934*x+0.00007565*x*x",0.,1200.)),
62 fZEMp(new TF1("fZEMp","80.05-0.1315*x+0.00005327*x*x",0.,1200.)),
63 fZEMsp(new TF1("fZEMsp","201.7-0.325*x+0.0001292*x*x",0.,1200.)),
64 fZEMb(new TF1("fZEMb",
65 "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.)),
66 //
6024ec85 67 fPedData(GetPedData()),
68 fECalibData(GetECalibData()),
69 fRecParam(GetRecParams())
f5d41205 70{
71 // **** Default constructor
72
73}
74
75
76//_____________________________________________________________________________
77AliZDCReconstructor::~AliZDCReconstructor()
78{
79// destructor
80
81 delete fZNCen;
82 delete fZNPer;
83 delete fZPCen;
84 delete fZPPer;
85 delete fZDCCen;
86 delete fZDCPer;
87 delete fbCen;
88 delete fbPer;
89 delete fZEMn;
90 delete fZEMp;
91 delete fZEMsp;
92 delete fZEMb;
93
94}
95
96
97//_____________________________________________________________________________
98void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
99{
100 // *** Local ZDC reconstruction for digits
101 // Works on the current event
102
103 // Retrieving calibration data
a85132e7 104 Float_t meanPed[48];
105 for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
f5d41205 106
107 // get digits
108 AliZDCDigit digit;
109 AliZDCDigit* pdigit = &digit;
110 digitsTree->SetBranchAddress("ZDC", &pdigit);
111
112 // loop over digits
113 Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
a85132e7 114 Float_t dZEM1CorrHG=0., dZEM2CorrHG=0.;
f5d41205 115 Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
116 Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
a85132e7 117 Float_t dZEM1CorrLG=0., dZEM2CorrLG=0.;
f5d41205 118 Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
119
120 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
121 for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
a85132e7 122 digitsTree->GetEntry(iDigit);
123 if (!pdigit) continue;
124 //pdigit->Print("");
125 //
126 Int_t det = digit.GetSector(0);
127 Int_t quad = digit.GetSector(1);
128 Int_t pedindex = -1;
129 //printf("\n\t Digit #%d det %d quad %d", iDigit, det, quad);
130 //
131 if(quad != 5){ // ZDC (not reference PTMs!)
f5d41205 132 if(det == 1){ // *** ZN1
133 pedindex = quad;
134 tZN1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
135 if(tZN1CorrHG[quad]<0.) tZN1CorrHG[quad] = 0.;
a85132e7 136 tZN1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
f5d41205 137 if(tZN1CorrLG[quad]<0.) tZN1CorrLG[quad] = 0.;
138 //printf("\t pedindex %d tZN1CorrHG[%d] = %1.0f tZN1CorrLG[%d] = %1.0f",
139 // pedindex, quad, tZN1CorrHG[quad], quad, tZN1CorrLG[quad]);
140 }
141 else if(det == 2){ // *** ZP1
a85132e7 142 pedindex = quad+5;
f5d41205 143 tZP1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
144 if(tZP1CorrLG[quad]<0.) tZP1CorrLG[quad] = 0.;
a85132e7 145 tZP1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
f5d41205 146 if(tZP1CorrHG[quad]<0.) tZP1CorrHG[quad] = 0.;
147 //printf("\t pedindex %d tZP1CorrHG[%d] = %1.0f tZP1CorrLG[%d] = %1.0f",
148 // pedindex, quad, tZP1CorrHG[quad], quad, tZP1CorrLG[quad]);
149 }
150 else if(det == 3){
151 if(quad == 1){ // *** ZEM1
a85132e7 152 pedindex = quad+9;
153 dZEM1CorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
154 if(dZEM1CorrHG<0.) dZEM1CorrHG = 0.;
155 dZEM1CorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
156 if(dZEM1CorrLG<0.) dZEM1CorrLG = 0.;
f5d41205 157 //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n",
a85132e7 158 // pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEM1CorrHG);
f5d41205 159 //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n",
a85132e7 160 // pedindex+24, digit.GetADCValue(1), meanPed[pedindex+24], dZEM1CorrLG);
f5d41205 161 }
a85132e7 162 else if(quad == 2){ // *** ZEM2
163 pedindex = quad+9;
164 dZEM2CorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
165 if(dZEM2CorrHG<0.) dZEM2CorrHG = 0.;
166 dZEM2CorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
167 if(dZEM2CorrLG<0.) dZEM2CorrLG = 0.;
f5d41205 168 //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n",
a85132e7 169 // pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEM2CorrHG);
f5d41205 170 //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n",
a85132e7 171 // pedindex+2, digit.GetADCValue(1),meanPed[pedindex+2], dZEM2CorrLG);
f5d41205 172 }
173 }
174 else if(det == 4){ // *** ZN2
a85132e7 175 pedindex = quad+12;
f5d41205 176 tZN2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
177 if(tZN2CorrHG[quad]<0.) tZN2CorrHG[quad] = 0.;
a85132e7 178 tZN2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
f5d41205 179 if(tZN2CorrLG[quad]<0.) tZN2CorrLG[quad] = 0.;
180 //printf("\t pedindex %d tZN2CorrHG[%d] = %1.0f tZN2CorrLG[%d] = %1.0f\n",
181 // pedindex, quad, tZN2CorrHG[quad], quad, tZN2CorrLG[quad]);
182 }
183 else if(det == 5){ // *** ZP2
a85132e7 184 pedindex = quad+17;
f5d41205 185 tZP2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
186 if(tZP2CorrHG[quad]<0.) tZP2CorrHG[quad] = 0.;
a85132e7 187 tZP2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
f5d41205 188 if(tZP2CorrLG[quad]<0.) tZP2CorrLG[quad] = 0.;
189 //printf("\t pedindex %d tZP2CorrHG[%d] = %1.0f tZP2CorrLG[%d] = %1.0f\n",
190 // pedindex, quad, tZP2CorrHG[quad], quad, tZP2CorrLG[quad]);
191 }
a85132e7 192 }
f5d41205 193 }
194
195 // reconstruct the event
196 ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG,
197 tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG,
a85132e7 198 tZP2CorrLG, dZEM1CorrHG, dZEM2CorrHG);
f5d41205 199
200}
201
202//_____________________________________________________________________________
203void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
204{
205 // *** ZDC raw data reconstruction
206 // Works on the current event
207
208 // Retrieving calibration data
a85132e7 209 Float_t meanPed[48];
210 for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
f5d41205 211
212 rawReader->Reset();
213
214 // loop over raw data rawDatas
215 Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
a85132e7 216 Float_t dZEM1CorrHG=0., dZEM2CorrHG=0.;
f5d41205 217 Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
218 Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
a85132e7 219 Float_t dZEM1CorrLG=0., dZEM2CorrLG=0.;
f5d41205 220 Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
221 //
222 AliZDCRawStream rawData(rawReader);
223 while (rawData.Next()) {
224 if(rawData.IsADCDataWord()){
a85132e7 225 Int_t det = rawData.GetSector(0);
226 Int_t quad = rawData.GetSector(1);
227 Int_t gain = rawData.GetADCGain();
228 Int_t pedindex=0;
229 //
230 if(quad !=5){ // ZDCs (not reference PTMs)
f5d41205 231 if(det == 1){
232 pedindex = quad;
233 if(gain == 0) tZN1CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
a85132e7 234 else tZN1CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
f5d41205 235 }
236 else if(det == 2){
a85132e7 237 pedindex = quad+5;
f5d41205 238 if(gain == 0) tZP1CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
a85132e7 239 else tZP1CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
f5d41205 240 }
241 else if(det == 3){
242 if(quad==1){
a85132e7 243 pedindex = quad+9;
244 if(gain == 0) dZEM1CorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
245 else dZEM1CorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
f5d41205 246 }
247 else if(quad==2){
a85132e7 248 pedindex = quad+9;
249 if(gain == 0) dZEM2CorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
250 else dZEM2CorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
f5d41205 251 }
252 }
253 else if(det == 4){
a85132e7 254 pedindex = quad+12;
f5d41205 255 if(gain == 0) tZN2CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
a85132e7 256 else tZN2CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
f5d41205 257 }
258 else if(det == 5){
a85132e7 259 pedindex = quad+17;
f5d41205 260 if(gain == 0) tZP2CorrHG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
a85132e7 261 else tZP2CorrLG[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
f5d41205 262 }
a85132e7 263 printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
264 det,quad,gain, pedindex, meanPed[pedindex]);
265 }
266 }//IsADCDataWord
f5d41205 267 }
268
269 // reconstruct the event
270 ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG,
271 tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG,
a85132e7 272 tZP2CorrLG, dZEM1CorrHG, dZEM2CorrHG);
f5d41205 273
274}
275
276//_____________________________________________________________________________
277void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree,
278 Float_t* ZN1ADCCorrHG, Float_t* ZP1ADCCorrHG,
279 Float_t* ZN2ADCCorrHG, Float_t* ZP2ADCCorrHG,
280 Float_t* ZN1ADCCorrLG, Float_t* ZP1ADCCorrLG,
281 Float_t* ZN2ADCCorrLG, Float_t* ZP2ADCCorrLG,
a85132e7 282 Float_t corrADCZEM1HG, Float_t corrADCZEM2HG) const
f5d41205 283{
284 // ***** Reconstruct one event
285
286 // *** RECONSTRUCTION FROM SIMULATED DATA
287 // It passes trhough the no. of phe which is known from simulations
288 // --- ADCchannel -> photoelectrons
289 // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
290 // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
291 //Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
292 //zn1phe = ZN1Corr/convFactor;
293 //zp1phe = ZP1Corr/convFactor;
294 //zemphe = ZEMCorr/convFactor;
295 //zn2phe = ZN2Corr/convFactor;
296 //zp2phe = ZP2Corr/convFactor;
297 ////if AliDebug(1,Form("\n znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
298 //
299 //// --- Energy calibration
300 //// Conversion factors for hadronic ZDCs goes from phe yield to TRUE
301 //// incident energy (conversion from GeV to TeV is included); while for EM
302 //// calos conversion is from light yield to detected energy calculated by
303 //// GEANT NB -> ZN and ZP conversion factors are constant since incident
304 //// spectators have all the same energy, ZEM energy is obtained through a
305 //// fit over the whole range of incident particle energies
306 //// (obtained with full HIJING simulations)
307 //Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
308 //Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
309 //zn1energy = zn1phe/zn1phexTeV;
310 //zp1energy = zp1phe/zp1phexTeV;
311 //zdc1energy = zn1energy+zp1energy;
312 //zn2energy = zn2phe/zn2phexTeV;
313 //zp2energy = zp2phe/zp2phexTeV;
314 //zdc2energy = zn2energy+zp2energy;
315 //zemenergy = -4.81+0.3238*zemphe;
316 //if(zemenergy<0) zemenergy=0;
317 //// if AliDebug(1,Form(" znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
318 //// "\n zemenergy = %f TeV\n", znenergy, zpenergy,
319 //// zdcenergy, zemenergy);
320 //// if(zdcenergy==0)
321 //// if AliDebug(1,Form("\n\n ### ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
322 //// " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy);
323
324 //
325 // *** RECONSTRUCTION FROM "REAL" DATA
326 //
327 // Retrieving calibration data
84d6255e 328 // --- Equalization coefficients ---------------------------------------------
f5d41205 329 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
330 for(Int_t ji=0; ji<5; ji++){
6024ec85 331 equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
332 equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
333 equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
334 equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
f5d41205 335 }
84d6255e 336 // --- Energy calibration factors ------------------------------------
f5d41205 337 Float_t calibEne[4];
6024ec85 338 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
f5d41205 339 //
84d6255e 340 // --- Reconstruction parameters ------------------
6024ec85 341 Float_t endPointZEM = fRecParam->GetZEMEndValue();
342 Float_t cutFractionZEM = fRecParam->GetZEMCutFraction();
343 Float_t dZEMSup = fRecParam->GetDZEMSup();
344 Float_t dZEMInf = fRecParam->GetDZEMInf();
f5d41205 345 //
346 Float_t cutValueZEM = endPointZEM*cutFractionZEM;
347 Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
348 Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
349 //
6024ec85 350 Float_t maxValEZN1 = fRecParam->GetEZN1MaxValue();
351 Float_t maxValEZP1 = fRecParam->GetEZP1MaxValue();
352 Float_t maxValEZDC1 = fRecParam->GetEZDC1MaxValue();
353 Float_t maxValEZN2 = fRecParam->GetEZN2MaxValue();
354 Float_t maxValEZP2 = fRecParam->GetEZP2MaxValue();
355 Float_t maxValEZDC2 = fRecParam->GetEZDC2MaxValue();
f5d41205 356 //
357 //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
358 // " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
359
360 // Equalization of detector responses
361 Float_t equalTowZN1HG[5], equalTowZN2HG[5], equalTowZP1HG[5], equalTowZP2HG[5];
362 Float_t equalTowZN1LG[5], equalTowZN2LG[5], equalTowZP1LG[5], equalTowZP2LG[5];
363 for(Int_t gi=0; gi<5; gi++){
364 equalTowZN1HG[gi] = ZN1ADCCorrHG[gi]*equalCoeffZN1[gi];
365 equalTowZP1HG[gi] = ZP1ADCCorrHG[gi]*equalCoeffZP1[gi];
366 equalTowZN2HG[gi] = ZN2ADCCorrHG[gi]*equalCoeffZN2[gi];
367 equalTowZP2HG[gi] = ZP2ADCCorrHG[gi]*equalCoeffZP2[gi];
368 //
369 equalTowZN1LG[gi] = ZN1ADCCorrLG[gi]*equalCoeffZN1[gi];
370 equalTowZP1LG[gi] = ZP1ADCCorrLG[gi]*equalCoeffZP1[gi];
371 equalTowZN2LG[gi] = ZN2ADCCorrLG[gi]*equalCoeffZN2[gi];
372 equalTowZP2LG[gi] = ZP2ADCCorrLG[gi]*equalCoeffZP2[gi];
373 }
374
375 // Energy calibration of detector responses
376 Float_t calibTowZN1HG[5], calibTowZN2HG[5], calibTowZP1HG[5], calibTowZP2HG[5];
377 Float_t calibSumZN1HG=0., calibSumZN2HG=0., calibSumZP1HG=0., calibSumZP2HG=0.;
378 Float_t calibTowZN1LG[5], calibTowZN2LG[5], calibTowZP1LG[5], calibTowZP2LG[5];
379 Float_t calibSumZN1LG=0., calibSumZN2LG=0., calibSumZ12LG=0., calibSumZP2LG=0.;
380 for(Int_t gi=0; gi<5; gi++){
381 calibTowZN1HG[gi] = equalTowZN1HG[gi]*calibEne[0];
382 calibTowZP1HG[gi] = equalTowZP1HG[gi]*calibEne[1];
383 calibTowZN2HG[gi] = equalTowZN2HG[gi]*calibEne[2];
384 calibTowZP2HG[gi] = equalTowZP2HG[gi]*calibEne[3];
385 calibSumZN1HG += calibTowZN1HG[gi];
386 calibSumZP1HG += calibTowZP1HG[gi];
387 calibSumZN2HG += calibTowZN2HG[gi];
388 calibSumZP2HG += calibTowZP2HG[gi];
389 //
390 calibTowZN1LG[gi] = equalTowZN1LG[gi]*calibEne[0];
391 calibTowZP1LG[gi] = equalTowZP1LG[gi]*calibEne[1];
392 calibTowZN2LG[gi] = equalTowZN2LG[gi]*calibEne[2];
393 calibTowZP2LG[gi] = equalTowZP2LG[gi]*calibEne[3];
394 calibSumZN1LG += calibTowZN1LG[gi];
395 calibSumZ12LG += calibTowZP1LG[gi];
396 calibSumZN2LG += calibTowZN2LG[gi];
397 calibSumZP2LG += calibTowZP2LG[gi];
398 }
399
400 // --- Number of detected spectator nucleons
401 // *** N.B. -> It works only in Pb-Pb
402 Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
403 nDetSpecNLeft = (Int_t) (calibSumZN1HG/2.760);
404 nDetSpecPLeft = (Int_t) (calibSumZP1HG/2.760);
405 nDetSpecNRight = (Int_t) (calibSumZN2HG/2.760);
406 nDetSpecPRight = (Int_t) (calibSumZP2HG/2.760);
407 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
408 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
409 nDetSpecNRight, nDetSpecPRight);*/
410
411 // --- Number of generated spectator nucleons (from HIJING parameterization)
412 Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
413 Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
414 Double_t impPar=0.;
415 //
416 // *** RECONSTRUCTION FROM SIMULATED DATA
417 // Cut value for Ezem (GeV)
418 // ### Results from production -> 0<b<18 fm (Apr 2002)
419 /*Float_t eZEMCut = 420.;
420 Float_t deltaEZEMSup = 690.;
421 Float_t deltaEZEMInf = 270.;
422 if(zemenergy > (eZEMCut+deltaEZEMSup)){
423 nGenSpecNLeft = (Int_t) (fZNCen->Eval(ZN1CalibSum));
424 nGenSpecPLeft = (Int_t) (fZPCen->Eval(ZP1CalibSum));
425 nGenSpecLeft = (Int_t) (fZDCCen->Eval(ZN1CalibSum+ZP1CalibSum));
426 nGenSpecNRight = (Int_t) (fZNCen->Eval(ZN2CalibSum));
427 nGenSpecPRight = (Int_t) (fZNCen->Eval(ZP2CalibSum));
428 nGenSpecRight = (Int_t) (fZNCen->Eval(ZN2CalibSum+ZP2CalibSum));
429 impPar = fbCen->Eval(ZN1CalibSum+ZP1CalibSum);
430 }
431 else if(zemenergy < (eZEMCut-deltaEZEMInf)){
432 nGenSpecNLeft = (Int_t) (fZNPer->Eval(ZN1CalibSum));
433 nGenSpecPLeft = (Int_t) (fZPPer->Eval(ZP1CalibSum));
434 nGenSpecLeft = (Int_t) (fZDCPer->Eval(ZN1CalibSum+ZP1CalibSum));
435 impPar = fbPer->Eval(ZN1CalibSum+ZP1CalibSum);
436 }
437 else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
438 nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
439 nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
440 nGenSpecLeft = (Int_t)(fZEMsp->Eval(zemenergy));
441 impPar = fZEMb->Eval(zemenergy);
442 }
443 // ### Results from production -> 0<b<18 fm (Apr 2002)
444 if(ZN1CalibSum>162.) nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
445 if(ZP1CalibSum>59.75) nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
446 if(ZN1CalibSum+ZP1CalibSum>221.5) nGenSpecLeft = (Int_t)(fZEMsp->Eval(zemenergy));
447 if(ZN1CalibSum+ZP1CalibSum>220.) impPar = fZEMb->Eval(zemenergy);
448 */
449 //
450 //
451 // *** RECONSTRUCTION FROM REAL DATA
c0aec6f6 452 //
a85132e7 453 Float_t corrADCZEMHG = corrADCZEM1HG + corrADCZEM2HG;
454 //
a4cab348 455 if(corrADCZEMHG > supValueZEM){
456 nGenSpecNLeft = (Int_t) (fZNCen->Eval(calibSumZN1HG));
457 nGenSpecPLeft = (Int_t) (fZPCen->Eval(calibSumZP1HG));
458 nGenSpecLeft = (Int_t) (fZDCCen->Eval(calibSumZN1HG+calibSumZP1HG));
459 nGenSpecNRight = (Int_t) (fZNCen->Eval(calibSumZN2HG));
460 nGenSpecPRight = (Int_t) (fZNCen->Eval(calibSumZP2HG));
461 nGenSpecRight = (Int_t) (fZNCen->Eval(calibSumZN2HG+calibSumZP2HG));
462 impPar = fbCen->Eval(calibSumZN1HG+calibSumZP1HG);
646f1679 463 }
a4cab348 464 else if(corrADCZEMHG < infValueZEM){
465 nGenSpecNLeft = (Int_t) (fZNPer->Eval(calibSumZN1HG));
466 nGenSpecPLeft = (Int_t) (fZPPer->Eval(calibSumZP1HG));
467 nGenSpecLeft = (Int_t) (fZDCPer->Eval(calibSumZN1HG+calibSumZP1HG));
468 impPar = fbPer->Eval(calibSumZN1HG+calibSumZP1HG);
646f1679 469 }
a4cab348 470 else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
471 nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
472 nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
473 nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
474 impPar = fZEMb->Eval(corrADCZEMHG);
646f1679 475 }
476 //
a4cab348 477 if(calibSumZN1HG/maxValEZN1>1.) nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
478 if(calibSumZP1HG/maxValEZP1>1.) nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
479 if((calibSumZN1HG+calibSumZP1HG/maxValEZDC1)>1.){
480 nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
481 impPar = fZEMb->Eval(corrADCZEMHG);
646f1679 482 }
a4cab348 483 if(calibSumZN2HG/maxValEZN2>1.) nGenSpecNRight = (Int_t) (fZEMn->Eval(corrADCZEMHG));
484 if(calibSumZP2HG/maxValEZP2>1.) nGenSpecPRight = (Int_t) (fZEMp->Eval(corrADCZEMHG));
485 if((calibSumZN2HG+calibSumZP2HG/maxValEZDC2)>1.) nGenSpecRight = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
646f1679 486 //
487 if(nGenSpecNLeft>125) nGenSpecNLeft=125;
488 else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
489 if(nGenSpecPLeft>82) nGenSpecPLeft=82;
490 else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
491 if(nGenSpecLeft>207) nGenSpecLeft=207;
492 else if(nGenSpecLeft<0) nGenSpecLeft=0;
8309c1ab 493
980685f2 494 // --- Number of generated participants (from HIJING parameterization)
646f1679 495 Int_t nPart, nPartTotLeft, nPartTotRight;
496 nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
497 nPartTotLeft = 207-nGenSpecLeft;
498 nPartTotRight = 207-nGenSpecRight;
e90a5fef 499 if(nPart<0) nPart=0;
500 if(nPartTotLeft<0) nPartTotLeft=0;
501 if(nPartTotRight<0) nPartTotRight=0;
e97af564 502 //
e90a5fef 503 // *** DEBUG ***
422ed0d5 504 printf("\n\t AliZDCReconstructor -> calibSumZN1HG %1.0f, calibSumZP1HG %1.0f,"
e90a5fef 505 " calibSumZN2HG %1.0f, calibSumZP2HG %1.0f, corrADCZEMHG %1.0f\n",
506 calibSumZN1HG,calibSumZP1HG,calibSumZN2HG,calibSumZP2HG,corrADCZEMHG);
507 printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
508 "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n",
509 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft,
510 nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
511 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
422ed0d5 512
646f1679 513 // create the output tree
a4cab348 514 AliZDCReco reco(calibSumZN1HG, calibSumZP1HG, calibSumZN2HG, calibSumZP2HG,
515 calibTowZN1LG, calibTowZN2LG, calibTowZP1LG, calibTowZP2LG,
84d6255e 516 calibTowZN1LG, calibTowZP1LG, calibTowZN2LG, calibTowZP2LG,
a85132e7 517 corrADCZEM1HG, corrADCZEM2HG,
646f1679 518 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
519 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
520 nGenSpecPRight, nGenSpecRight,
521 nPartTotLeft, nPartTotRight, impPar);
522
8309c1ab 523 AliZDCReco* preco = &reco;
524 const Int_t kBufferSize = 4000;
70f04f6d 525 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
8309c1ab 526
527 // write the output tree
70f04f6d 528 clustersTree->Fill();
8309c1ab 529}
530
531//_____________________________________________________________________________
70f04f6d 532void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
8309c1ab 533{
70f04f6d 534 // fill energies and number of participants to the ESD
8309c1ab 535
8309c1ab 536 AliZDCReco reco;
537 AliZDCReco* preco = &reco;
70f04f6d 538 clustersTree->SetBranchAddress("ZDC", &preco);
8309c1ab 539
70f04f6d 540 clustersTree->GetEntry(0);
84d6255e 541 //
a85132e7 542 AliESDZDC * esdzdc = esd->GetESDZDC();
543 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
544 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
545 for(Int_t i=0; i<5; i++){
e90a5fef 546 tZN1Ene[i] = reco.GetZN1EnTow(i);
547 tZN2Ene[i] = reco.GetZN2EnTow(i);
a85132e7 548 tZP1Ene[i] = reco.GetZP1EnTow(i);
549 tZP2Ene[i] = reco.GetZP2EnTow(i);
550 tZN1EneLR[i] = reco.GetZN1SigLowRes(i);
551 tZN2EneLR[i] = reco.GetZP1SigLowRes(i);
552 tZP1EneLR[i] = reco.GetZN2SigLowRes(i);
553 tZP2EneLR[i] = reco.GetZP2SigLowRes(i);
e90a5fef 554 }
a85132e7 555 esdzdc->SetZN1TowerEnergy(tZN1Ene);
556 esdzdc->SetZN2TowerEnergy(tZN2Ene);
557 esdzdc->SetZP1TowerEnergy(tZP1Ene);
558 esdzdc->SetZP2TowerEnergy(tZP2Ene);
559 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
560 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
561 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
562 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
563 //
564 esd->SetZDC(reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEM1signal(),
565 reco.GetZEM2signal(), reco.GetZN2Energy(), reco.GetZP2Energy(),
646f1679 566 reco.GetNPartLeft());
a85132e7 567 //
a4cab348 568
8309c1ab 569}
48642b09 570
571//_____________________________________________________________________________
78d18275 572AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
48642b09 573{
cc2abffd 574 // Setting the storage
48642b09 575
78d18275 576 Bool_t deleteManager = kFALSE;
48642b09 577
78d18275 578 AliCDBManager *manager = AliCDBManager::Instance();
579 AliCDBStorage *defstorage = manager->GetDefaultStorage();
48642b09 580
78d18275 581 if(!defstorage || !(defstorage->Contains("ZDC"))){
582 AliWarning("No default storage set or default storage doesn't contain ZDC!");
583 manager->SetDefaultStorage(uri);
584 deleteManager = kTRUE;
585 }
586
587 AliCDBStorage *storage = manager->GetDefaultStorage();
588
589 if(deleteManager){
590 AliCDBManager::Instance()->UnsetDefaultStorage();
591 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
592 }
593
594 return storage;
595}
48642b09 596
78d18275 597//_____________________________________________________________________________
6024ec85 598AliZDCPedestals* AliZDCReconstructor::GetPedData() const
599{
600
601 // Getting pedestal calibration object for ZDC set
602
603 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
604 if(!entry) AliFatal("No calibration data loaded!");
605
606 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
607 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
608
609 return calibdata;
610}
611
612//_____________________________________________________________________________
613AliZDCCalib* AliZDCReconstructor::GetECalibData() const
614{
615
616 // Getting energy and equalization calibration object for ZDC set
617
618 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Calib");
619 if(!entry) AliFatal("No calibration data loaded!");
620
621 AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*> (entry->GetObject());
622 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
623
624 return calibdata;
625}
626
627//_____________________________________________________________________________
628AliZDCRecParam* AliZDCReconstructor::GetRecParams() const
78d18275 629{
48642b09 630
6024ec85 631 // Getting energy and equalization calibration object for ZDC set
4fda3ba1 632
6024ec85 633 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecParam");
457a440d 634 if(!entry) AliFatal("No calibration data loaded!");
4fda3ba1 635
6024ec85 636 AliZDCRecParam *calibdata = dynamic_cast<AliZDCRecParam*> (entry->GetObject());
457a440d 637 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
48642b09 638
78d18275 639 return calibdata;
48642b09 640}