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