Setting "const" to member functions of AliPHOSEMCAGeometry
[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx
... / ...
CommitLineData
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#include <TMap.h>
27
28#include "AliRunLoader.h"
29#include "AliRawReader.h"
30#include "AliGRPObject.h"
31#include "AliESDEvent.h"
32#include "AliESDZDC.h"
33#include "AliZDCDigit.h"
34#include "AliZDCRawStream.h"
35#include "AliZDCReco.h"
36#include "AliZDCReconstructor.h"
37#include "AliZDCPedestals.h"
38#include "AliZDCCalib.h"
39#include "AliZDCRecoParam.h"
40#include "AliZDCRecoParampp.h"
41#include "AliZDCRecoParamPbPb.h"
42
43
44ClassImp(AliZDCReconstructor)
45AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0; //reconstruction parameters
46
47//_____________________________________________________________________________
48AliZDCReconstructor:: AliZDCReconstructor() :
49 fPedData(GetPedData()),
50 fECalibData(GetECalibData()),
51 fRecoMode(0),
52 fBeamEnergy(0.)
53{
54 // **** Default constructor
55 SetRecoMode();
56
57}
58
59
60//_____________________________________________________________________________
61AliZDCReconstructor::~AliZDCReconstructor()
62{
63// destructor
64 if(fRecoParam) delete fRecoParam;
65 if(fPedData) delete fPedData;
66 if(fECalibData) delete fECalibData;
67}
68
69//____________________________________________________________________________
70void AliZDCReconstructor::SetRecoMode()
71{
72 // Setting reconstruction mode
73
74 // Initialization of the GRP entry
75 AliGRPObject* grpData;
76 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
77 if(entry){
78 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
79 if(m){
80 m->Print();
81 grpData = new AliGRPObject();
82 grpData->ReadValuesFromMap(m);
83 }
84 else{
85 grpData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
86 entry->SetOwner(0);
87 }
88 AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
89 }
90 if(!grpData) AliError("No GRP entry found in OCDB!");
91
92 TString beamType = grpData->GetBeamType();
93 if(beamType==AliGRPObject::GetInvalidString()){
94 AliError("GRP/GRP/Data entry: missing value for the beam energy !");
95 AliError("\t ZDC does not reconstruct event 4 UNKNOWN beam type\n");
96 return;
97 }
98 //
99 if((beamType.CompareTo("p-p")) == 0){
100 fRecoMode=0;
101 fRecoParam = (AliZDCRecoParampp*) AliZDCRecoParampp::GetppRecoParam();
102 }
103 else if((beamType.CompareTo("A-A")) == 0){
104 fRecoMode=1;
105 fRecoParam = (AliZDCRecoParamPbPb*) AliZDCRecoParamPbPb::GetPbPbRecoParam();
106 }
107
108 fBeamEnergy = grpData->GetBeamEnergy();
109 if(fBeamEnergy==AliGRPObject::GetInvalidFloat()) {
110 AliError("GRP/GRP/Data entry: missing value for the beam energy ! Using 0");
111 fBeamEnergy = 0.;
112 }
113
114 printf("\n ***** ZDC reconstruction initialized for %s @ %1.3f GeV\n\n",beamType.Data(), fBeamEnergy);
115}
116
117//_____________________________________________________________________________
118void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
119{
120 // *** Local ZDC reconstruction for digits
121 // Works on the current event
122
123 // Retrieving calibration data
124 Float_t meanPed[48];
125 for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
126
127 // get digits
128 AliZDCDigit digit;
129 AliZDCDigit* pdigit = &digit;
130 digitsTree->SetBranchAddress("ZDC", &pdigit);
131 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
132
133 // loop over digits
134 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
135 Float_t dZEM1Corr[2], dZEM2Corr[2], PMRef1[2], PMRef2[2];
136 for(Int_t i=0; i<10; i++){
137 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
138 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = PMRef1[i] = PMRef2[i] = 0.;
139 }
140 //
141 for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
142 digitsTree->GetEntry(iDigit);
143 if (!pdigit) continue;
144 //
145 Int_t det = digit.GetSector(0);
146 Int_t quad = digit.GetSector(1);
147 Int_t pedindex = -1, kNch = 24;
148 //printf("\n\t Digit #%d det %d quad %d", iDigit, det, quad);
149 //
150 if(quad != 5){ // ZDC (not reference PTMs!)
151 if(det == 1){ // *** ZNC
152 pedindex = quad;
153 tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
154 if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
155 tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
156 if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
157 //printf("\t pedindex %d tZN1Corr[%d] = %1.0f tZN1Corr[%d] = %1.0f",
158 // pedindex, quad, tZN1Corr[quad], quad+5, tZN1Corr[quad+5]);
159 }
160 else if(det == 2){ // *** ZP1
161 pedindex = quad+5;
162 tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
163 if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
164 tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
165 if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
166 //printf("\t pedindex %d tZP1Corr[%d] = %1.0f tZP1Corr[%d] = %1.0f",
167 // pedindex, quad, tZP1Corr[quad], quad+5, tZP1Corr[quad+5]);
168 }
169 else if(det == 3){
170 pedindex = quad+9;
171 if(quad == 1){ // *** ZEM1
172 dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
173 if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
174 dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
175 if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
176 //printf("\t pedindex %d tZEM1Corr[%d] = %1.0f tZEM1Corr[%d] = %1.0f",
177 // pedindex, quad, tZEM1Corr[quad], quad+1, tZEM1Corr[quad+1]);
178 }
179 else if(quad == 2){ // *** ZEM2
180 dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
181 if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
182 dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
183 if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
184 //printf("\t pedindex %d tZEM2Corr[%d] = %1.0f tZEM2Corr[%d] = %1.0f",
185 // pedindex, quad, tZEM2Corr[quad], quad+1, tZEM2Corr[quad+1]);
186 }
187 }
188 else if(det == 4){ // *** ZN2
189 pedindex = quad+12;
190 tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
191 if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
192 tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
193 if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
194 //printf("\t pedindex %d tZN2Corr[%d] = %1.0f tZN2Corr[%d] = %1.0f\n",
195 // pedindex, quad, tZN2Corr[quad], quad+5, tZN2Corr[quad+5]);
196 }
197 else if(det == 5){ // *** ZP2
198 pedindex = quad+17;
199 tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
200 if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
201 tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
202 if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
203 //printf("\t pedindex %d tZP2Corr[%d] = %1.0f tZP2Corr[%d] = %1.0f\n",
204 // pedindex, quad, tZP2Corr[quad], quad+5, tZP2Corr[quad+5]);
205 }
206 }
207 else{ // Reference PMs
208 pedindex = (det-1)/3+22;
209 if(det == 1){
210 PMRef1[0] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
211 if(PMRef1[0]<0.) PMRef1[0] = 0.;
212 PMRef1[1] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
213 if(PMRef2[1]<0.) PMRef1[1] = 0.;
214 }
215 else if(det == 4){
216 PMRef2[0] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
217 if(PMRef2[0]<0.) PMRef2[0] = 0.;
218 PMRef2[1] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
219 if(PMRef2[1]<0.) PMRef2[1] = 0.;
220 }
221 }
222 }
223
224 // reconstruct the event
225 if(fRecoMode==0)
226 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
227 dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
228 else if(fRecoMode==1)
229 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
230 dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
231
232}
233
234//_____________________________________________________________________________
235void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
236{
237 // *** ZDC raw data reconstruction
238 // Works on the current event
239
240 // Retrieving calibration data
241 Float_t meanPed[48];
242 for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
243
244 rawReader->Reset();
245
246 // loop over raw data
247 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
248 Float_t dZEM1Corr[2], dZEM2Corr[2], PMRef1[2], PMRef2[2];
249 for(Int_t i=0; i<10; i++){
250 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
251 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = PMRef1[i] = PMRef2[i] = 0.;
252 }
253 //
254 AliZDCRawStream rawData(rawReader);
255 Int_t kNch = 24;
256 while (rawData.Next()) {
257 if(rawData.IsADCDataWord()){
258 Int_t det = rawData.GetSector(0);
259 Int_t quad = rawData.GetSector(1);
260 Int_t gain = rawData.GetADCGain();
261 Int_t pedindex=0;
262 //
263 if(quad !=5){ // ZDCs (not reference PTMs)
264 if(det == 1){
265 pedindex = quad;
266 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
267 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
268 }
269 else if(det == 2){
270 pedindex = quad+5;
271 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
272 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
273 }
274 else if(det == 3){
275 pedindex = quad+9;
276 if(quad==1){
277 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
278 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
279 }
280 else if(quad==2){
281 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
282 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
283 }
284 }
285 else if(det == 4){
286 pedindex = quad+12;
287 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
288 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
289 }
290 else if(det == 5){
291 pedindex = quad+17;
292 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
293 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
294 }
295 //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
296 // det,quad,gain, pedindex, meanPed[pedindex]);
297 }
298 else{ // reference PM
299 pedindex = (det-1)/3 + 22;
300 if(det == 1){
301 if(gain==0) PMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
302 else PMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
303 }
304 else if(det ==4){
305 if(gain==0) PMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
306 else PMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
307 }
308 }
309 }//IsADCDataWord
310 }
311
312 // reconstruct the event
313 if(fRecoMode==0)
314 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
315 dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
316 else if(fRecoMode==1)
317 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
318 dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
319
320}
321
322//_____________________________________________________________________________
323void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* ZN1ADCCorr,
324 Float_t* ZP1ADCCorr, Float_t* ZN2ADCCorr, Float_t* ZP2ADCCorr,
325 Float_t* ZEM1ADCCorr, Float_t* ZEM2ADCCorr, Float_t* PMRef1, Float_t* PMRef2) const
326{
327 // ***** Reconstruct one event
328
329 // *** RECONSTRUCTION FROM "REAL" DATA
330 //
331 // Retrieving calibration data
332 // --- Equalization coefficients ---------------------------------------------
333 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
334 for(Int_t ji=0; ji<5; ji++){
335 equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
336 equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
337 equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
338 equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
339 }
340 // --- Energy calibration factors ------------------------------------
341 Float_t calibEne[4];
342 // *********************************************************************
343 // **** Until the beam type info isn't known @ reconstruction level ****
344 // **** the energy calibration coefficient are manually set to 1 ****
345 // **** as it will be in real life for pp data taking ****
346 // *********************************************************************
347 //for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
348 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = 1.;
349
350 // Equalization of detector responses
351 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
352 for(Int_t gi=0; gi<5; gi++){
353 equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
354 equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
355 equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
356 equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
357 equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
358 equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
359 equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
360 equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
361 }
362
363 // Energy calibration of detector responses
364 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
365 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
366 for(Int_t gi=0; gi<10; gi++){
367 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
368 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
369 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
370 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
371 //
372 if(gi<5){
373 calibSumZN1[0] += calibTowZN1[gi];
374 calibSumZP1[0] += calibTowZP1[gi];
375 calibSumZN2[0] += calibTowZN2[gi];
376 calibSumZP2[0] += calibTowZP2[gi];
377 }
378 //
379 else{
380 calibSumZN1[1] += calibTowZN1[gi];
381 calibSumZP1[1] += calibTowZP1[gi];
382 calibSumZN2[1] += calibTowZN2[gi];
383 calibSumZP2[1] += calibTowZP2[gi];
384 }
385 }
386
387 // --- Number of detected spectator nucleons
388 // *** N.B. -> It works only in Pb-Pb
389 Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
390 if(fBeamEnergy!=0){
391 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
392 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
393 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
394 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
395 }
396 else AliWarning(" ATTENTION -> fBeamEnergy = 0\n");
397 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
398 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
399 nDetSpecNRight, nDetSpecPRight);*/
400
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;
404 Int_t nPartTotLeft=0, nPartTotRight=0;
405 Double_t impPar=0.;
406
407 // create the output tree
408 AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
409 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
410 ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
411 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
412 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
413 nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
414
415 AliZDCReco* preco = &reco;
416 const Int_t kBufferSize = 4000;
417 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
418
419 // write the output tree
420 clustersTree->Fill();
421}
422
423//_____________________________________________________________________________
424void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree, Float_t* ZN1ADCCorr,
425 Float_t* ZP1ADCCorr, Float_t* ZN2ADCCorr, Float_t* ZP2ADCCorr,
426 Float_t* ZEM1ADCCorr, Float_t* ZEM2ADCCorr, Float_t* PMRef1, Float_t* PMRef2) const
427{
428 // ***** Reconstruct one event
429
430 // *** RECONSTRUCTION FROM "REAL" DATA
431 //
432 // Retrieving calibration data
433 // --- Equalization coefficients ---------------------------------------------
434 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
435 for(Int_t ji=0; ji<5; ji++){
436 equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
437 equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
438 equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
439 equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
440 }
441 // --- Energy calibration factors ------------------------------------
442 Float_t calibEne[4];
443 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
444
445 // Equalization of detector responses
446 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
447 for(Int_t gi=0; gi<5; gi++){
448 equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
449 equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
450 equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
451 equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
452 equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
453 equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
454 equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
455 equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
456 }
457
458 // Energy calibration of detector responses
459 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
460 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
461 for(Int_t gi=0; gi<10; gi++){
462 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
463 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
464 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
465 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
466 //
467 if(gi<5){
468 calibSumZN1[0] += calibTowZN1[gi];
469 calibSumZP1[0] += calibTowZP1[gi];
470 calibSumZN2[0] += calibTowZN2[gi];
471 calibSumZP2[0] += calibTowZP2[gi];
472 }
473 //
474 else{
475 calibSumZN1[1] += calibTowZN1[gi];
476 calibSumZP1[1] += calibTowZP1[gi];
477 calibSumZN2[1] += calibTowZN2[gi];
478 calibSumZP2[1] += calibTowZP2[gi];
479 }
480 }
481
482 //
483 // --- Reconstruction parameters ------------------
484 if(!fRecoParam) fRecoParam = (AliZDCRecoParamPbPb*) AliZDCRecoParamPbPb::GetPbPbRecoParam();
485 //
486 Float_t endPointZEM = fRecoParam->GetZEMEndValue();
487 Float_t cutFractionZEM = fRecoParam->GetZEMCutFraction();
488 Float_t dZEMSup = fRecoParam->GetDZEMSup();
489 Float_t dZEMInf = fRecoParam->GetDZEMInf();
490 //
491 Float_t cutValueZEM = endPointZEM*cutFractionZEM;
492 Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
493 Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
494 //
495 Float_t maxValEZN1 = fRecoParam->GetEZN1MaxValue();
496 Float_t maxValEZP1 = fRecoParam->GetEZP1MaxValue();
497 Float_t maxValEZDC1 = fRecoParam->GetEZDC1MaxValue();
498 Float_t maxValEZN2 = fRecoParam->GetEZN2MaxValue();
499 Float_t maxValEZP2 = fRecoParam->GetEZP2MaxValue();
500 Float_t maxValEZDC2 = fRecoParam->GetEZDC2MaxValue();
501 //
502 //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
503 // " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
504
505 // --- Number of detected spectator nucleons
506 // *** N.B. -> It works only in Pb-Pb
507 Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
508 if(fBeamEnergy!=0){
509 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
510 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
511 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
512 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
513 }
514 else AliWarning(" ATTENTION -> fBeamEnergy = 0\n");
515 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
516 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
517 nDetSpecNRight, nDetSpecPRight);*/
518
519 // --- Number of generated spectator nucleons (from HIJING parameterization)
520 Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
521 Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
522 Double_t impPar=0.;
523 //
524 Float_t corrADCZEMHG = ZEM1ADCCorr[0] + ZEM2ADCCorr[0];
525 //
526 if(corrADCZEMHG > supValueZEM){
527 nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN1[0]));
528 nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZPCen())->Eval(calibSumZP1[0]));
529 nGenSpecLeft = (Int_t) ((fRecoParam->GetfZDCCen())->Eval(calibSumZN1[0]+calibSumZP1[0]));
530 nGenSpecNRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN2[0]));
531 nGenSpecPRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZP2[0]));
532 nGenSpecRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN2[0]+calibSumZP2[0]));
533 impPar = (fRecoParam->GetfbCen())->Eval(calibSumZN1[0]+calibSumZP1[0]);
534 }
535 else if(corrADCZEMHG < infValueZEM){
536 nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZNPer())->Eval(calibSumZN1[0]));
537 nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZPPer())->Eval(calibSumZP1[0]));
538 nGenSpecLeft = (Int_t) ((fRecoParam->GetfZDCPer())->Eval(calibSumZN1[0]+calibSumZP1[0]));
539 impPar = (fRecoParam->GetfbPer())->Eval(calibSumZN1[0]+calibSumZP1[0]);
540 }
541 else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
542 nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
543 nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
544 nGenSpecLeft = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
545 impPar = (fRecoParam->GetfZEMb())->Eval(corrADCZEMHG);
546 }
547 //
548 if(calibSumZN1[0]/maxValEZN1>1.) nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
549 if(calibSumZP1[0]/maxValEZP1>1.) nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
550 if((calibSumZN1[0]+calibSumZP1[0]/maxValEZDC1)>1.){
551 nGenSpecLeft = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
552 impPar = (fRecoParam->GetfZEMb())->Eval(corrADCZEMHG);
553 }
554 if(calibSumZN2[0]/maxValEZN2>1.) nGenSpecNRight = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
555 if(calibSumZP2[0]/maxValEZP2>1.) nGenSpecPRight = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
556 if((calibSumZN2[0]+calibSumZP2[0]/maxValEZDC2)>1.) nGenSpecRight = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
557 //
558 if(nGenSpecNLeft>125) nGenSpecNLeft=125;
559 else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
560 if(nGenSpecPLeft>82) nGenSpecPLeft=82;
561 else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
562 if(nGenSpecLeft>207) nGenSpecLeft=207;
563 else if(nGenSpecLeft<0) nGenSpecLeft=0;
564
565 // --- Number of generated participants (from HIJING parameterization)
566 Int_t nPart, nPartTotLeft, nPartTotRight;
567 nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
568 nPartTotLeft = 207-nGenSpecLeft;
569 nPartTotRight = 207-nGenSpecRight;
570 if(nPart<0) nPart=0;
571 if(nPartTotLeft<0) nPartTotLeft=0;
572 if(nPartTotRight<0) nPartTotRight=0;
573 //
574 // *** DEBUG ***
575 /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
576 " calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
577 calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
578 printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
579 "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n",
580 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft,
581 nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
582 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
583 */
584
585 // create the output tree
586 AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
587 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
588 ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
589 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
590 nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
591 nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
592
593 AliZDCReco* preco = &reco;
594 const Int_t kBufferSize = 4000;
595 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
596
597 // write the output tree
598 clustersTree->Fill();
599}
600
601//_____________________________________________________________________________
602void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
603{
604 // fill energies and number of participants to the ESD
605
606 AliZDCReco reco;
607 AliZDCReco* preco = &reco;
608 clustersTree->SetBranchAddress("ZDC", &preco);
609
610 clustersTree->GetEntry(0);
611 //
612 AliESDZDC * esdzdc = esd->GetESDZDC();
613 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
614 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
615 for(Int_t i=0; i<5; i++){
616 tZN1Ene[i] = reco.GetZN1HREnTow(i);
617 tZN2Ene[i] = reco.GetZN2HREnTow(i);
618 tZP1Ene[i] = reco.GetZP1HREnTow(i);
619 tZP2Ene[i] = reco.GetZP2HREnTow(i);
620 //
621 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
622 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
623 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
624 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
625 }
626 esdzdc->SetZN1TowerEnergy(tZN1Ene);
627 esdzdc->SetZN2TowerEnergy(tZN2Ene);
628 esdzdc->SetZP1TowerEnergy(tZP1Ene);
629 esdzdc->SetZP2TowerEnergy(tZP2Ene);
630 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
631 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
632 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
633 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
634 //
635 esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
636 reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
637 reco.GetNPartLeft(), reco.GetNPartRight());
638 //
639
640}
641
642//_____________________________________________________________________________
643AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
644{
645 // Setting the storage
646
647 Bool_t deleteManager = kFALSE;
648
649 AliCDBManager *manager = AliCDBManager::Instance();
650 AliCDBStorage *defstorage = manager->GetDefaultStorage();
651
652 if(!defstorage || !(defstorage->Contains("ZDC"))){
653 AliWarning("No default storage set or default storage doesn't contain ZDC!");
654 manager->SetDefaultStorage(uri);
655 deleteManager = kTRUE;
656 }
657
658 AliCDBStorage *storage = manager->GetDefaultStorage();
659
660 if(deleteManager){
661 AliCDBManager::Instance()->UnsetDefaultStorage();
662 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
663 }
664
665 return storage;
666}
667
668//_____________________________________________________________________________
669AliZDCPedestals* AliZDCReconstructor::GetPedData() const
670{
671
672 // Getting pedestal calibration object for ZDC set
673
674 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
675 if(!entry) AliFatal("No calibration data loaded!");
676
677 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
678 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
679
680 return calibdata;
681}
682
683//_____________________________________________________________________________
684AliZDCCalib* AliZDCReconstructor::GetECalibData() const
685{
686
687 // Getting energy and equalization calibration object for ZDC set
688
689 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EMDCalib");
690 if(!entry) AliFatal("No calibration data loaded!");
691
692 AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*> (entry->GetObject());
693 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
694
695 return calibdata;
696}
697