Removing warnings
[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx
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 //                  Author: Chiara.Oppedisano@to.infn.it                     //
22 //                                                                           //
23 // NOTATIONS ADOPTED TO IDENTIFY DETECTORS (used in different ages!):        //
24 //   (ZN1,ZP1) or (ZNC, ZPC) or RIGHT refers to side C (RB26)                //
25 //   (ZN2,ZP2) or (ZNA, ZPA) or LEFT refers to side A (RB24)                 //
26 //                                                                           //
27 ///////////////////////////////////////////////////////////////////////////////
28
29
30 #include <TH2F.h>
31 #include <TH1D.h>
32 #include <TAxis.h>
33 #include <TMap.h>
34
35 #include "AliRawReader.h"
36 #include "AliESDEvent.h"
37 #include "AliESDZDC.h"
38 #include "AliZDCDigit.h"
39 #include "AliZDCRawStream.h"
40 #include "AliZDCReco.h"
41 #include "AliZDCReconstructor.h"
42 #include "AliZDCPedestals.h"
43 #include "AliZDCEnCalib.h"
44 #include "AliZDCTowerCalib.h"
45 #include "AliZDCRecoParam.h"
46 #include "AliZDCRecoParampp.h"
47 #include "AliZDCRecoParamPbPb.h"
48 #include "AliRunInfo.h"
49
50
51 ClassImp(AliZDCReconstructor)
52 AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0;  //reconstruction parameters
53
54 //_____________________________________________________________________________
55 AliZDCReconstructor:: AliZDCReconstructor() :
56   fPedData(GetPedestalData()),
57   fEnCalibData(GetEnergyCalibData()),
58   fTowCalibData(GetTowerCalibData()),
59   fRecoMode(0),
60   fBeamEnergy(0.),
61   fNRun(0),
62   fIsCalibrationMB(kFALSE),
63   fPedSubMode(0),
64   fSignalThreshold(7)
65 {
66   // **** Default constructor
67 }
68
69
70 //_____________________________________________________________________________
71 AliZDCReconstructor::~AliZDCReconstructor()
72 {
73 // destructor
74 //   if(fRecoParam)    delete fRecoParam;
75    if(fPedData)      delete fPedData;    
76    if(fEnCalibData)  delete fEnCalibData;
77    if(fTowCalibData) delete fTowCalibData;
78 }
79
80 //____________________________________________________________________________
81 void AliZDCReconstructor::Init()
82 {
83   // Setting reconstruction mode
84   // Getting beam type and beam energy from GRP calibration object
85   
86   TString runType = GetRunInfo()->GetRunType();
87   if((runType.CompareTo("CALIBRATION_MB")) == 0){
88     fIsCalibrationMB = kTRUE;
89   }
90     
91   TString beamType = GetRunInfo()->GetBeamType();
92   // This is a temporary solution to allow reconstruction in tests without beam
93   if(((beamType.CompareTo("UNKNOWN"))==0) && ((runType.CompareTo("PHYSICS")) == 0)){
94     fRecoMode=1;
95   }
96   /*else if((beamType.CompareTo("UNKNOWN"))==0){
97     AliError("\t UNKNOWN beam type\n");
98     return;
99   }*/
100
101   if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
102      ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)){
103     fRecoMode=1;
104   }
105   else if((beamType.CompareTo("A-A")) == 0){
106     fRecoMode=2;
107   }
108     
109   fBeamEnergy = GetRunInfo()->GetBeamEnergy();
110   if(fBeamEnergy==0.){
111     AliWarning(" Beam energy value missing -> E_beam = 0");
112     fBeamEnergy = 0.;
113   }
114     
115   if(fIsCalibrationMB==kFALSE)  
116     printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f GeV *****\n\n",beamType.Data(), fBeamEnergy);
117   
118 }
119
120 //_____________________________________________________________________________
121 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
122 {
123   // *** Local ZDC reconstruction for digits
124   // Works on the current event
125     
126   // Retrieving calibration data  
127   // Parameters for mean value pedestal subtraction
128   int const kNch = 24;
129   Float_t meanPed[2*kNch];    
130   for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
131   // Parameters pedestal subtraction through correlation with out-of-time signals
132   Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
133   for(Int_t jj=0; jj<2*kNch; jj++){
134      corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
135      corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
136   }
137
138   // get digits
139   AliZDCDigit digit;
140   AliZDCDigit* pdigit = &digit;
141   digitsTree->SetBranchAddress("ZDC", &pdigit);
142   //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
143
144   // loop over digits
145   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10]; 
146   Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2]; 
147   for(Int_t i=0; i<10; i++){
148      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
149      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
150   }  
151   
152   Int_t digNentries = digitsTree->GetEntries();
153   Float_t ootDigi[kNch];
154   // -- Reading out-of-time signals (last kNch entries) for current event
155   if(fPedSubMode==1){
156     for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
157        ootDigi[iDigit] = digitsTree->GetEntry(iDigit);
158     }
159   }
160   
161   for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
162    digitsTree->GetEntry(iDigit);
163    if (!pdigit) continue;
164    //  
165    Int_t det = digit.GetSector(0);
166    Int_t quad = digit.GetSector(1);
167    Int_t pedindex = -1;
168    Float_t ped2SubHg=0., ped2SubLg=0.;
169    if(quad!=5){
170      if(det==1)      pedindex = quad;
171      else if(det==2) pedindex = quad+5;
172      else if(det==3) pedindex = quad+9;
173      else if(det==4) pedindex = quad+12;
174      else if(det==5) pedindex = quad+17;
175    }
176    else pedindex = (det-1)/3+22;
177    //
178    if(fPedSubMode==0){
179      ped2SubHg = meanPed[pedindex];
180      ped2SubLg = meanPed[pedindex+kNch];
181    }
182    else if(fPedSubMode==1){
183      ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
184      ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
185    }
186       
187    if(quad != 5){ // ZDC (not reference PTMs!)
188     if(det == 1){ // *** ZNC
189        tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
190        tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
191        if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
192        if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
193     }
194     else if(det == 2){ // *** ZP1
195        tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
196        tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
197        if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
198        if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
199     }
200     else if(det == 3){
201        if(quad == 1){       // *** ZEM1  
202          dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg); 
203          dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg); 
204          if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
205          if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
206        }
207        else if(quad == 2){  // *** ZEM2
208          dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg); 
209          dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg); 
210          if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
211          if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
212        }
213     }
214     else if(det == 4){  // *** ZN2
215        tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
216        tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
217        if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
218        if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
219    }
220     else if(det == 5){  // *** ZP2 
221        tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
222        tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
223        if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
224        if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
225     }
226    }
227    else{ // Reference PMs
228      if(det == 1){
229        sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
230        sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
231        // Ch. debug
232        if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
233        if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
234      }
235      else if(det == 4){
236        sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
237        sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
238        // Ch. debug
239        if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
240        if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
241      }
242    }
243
244    // Ch. debug
245    /*printf("AliZDCReconstructor: digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
246          iDigit, det, quad, ped2SubHg, ped2SubLg);
247    printf(" -> pedindex %d\n", pedindex);
248    printf("   HGChain -> RawDig %d DigCorr %1.2f", 
249         digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg); 
250    printf("   LGChain -> RawDig %d DigCorr %1.2f\n", 
251         digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);*/ 
252    
253   }//digits loop
254  
255   UInt_t counts[32];
256   for(Int_t jj=0; jj<32; jj++) counts[jj]=0;
257   
258   Int_t  evQualityBlock[4] = {1,0,0,0};
259   Int_t  triggerBlock[4] = {0,0,0,0};
260   Int_t  chBlock[3] = {0,0,0};
261   UInt_t puBits=0;
262   
263   // reconstruct the event
264   if(fRecoMode==1)
265     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr, 
266       dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, 
267       kFALSE, counts, 
268       evQualityBlock,  triggerBlock,  chBlock, puBits);
269   else if(fRecoMode==2)
270     ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr, 
271       dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, 
272       kFALSE, counts, 
273       evQualityBlock,  triggerBlock,  chBlock, puBits);
274 }
275
276 //_____________________________________________________________________________
277 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
278 {
279   // *** ZDC raw data reconstruction
280   // Works on the current event
281   
282   // Retrieving calibration data  
283   // Parameters for pedestal subtraction
284   int const kNch = 24;
285   Float_t meanPed[2*kNch];    
286   for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
287   // Parameters pedestal subtraction through correlation with out-of-time signals
288   Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
289   for(Int_t jj=0; jj<2*kNch; jj++){
290      corrCoeff0[jj] =  fPedData->GetPedCorrCoeff0(jj);
291      corrCoeff1[jj] =  fPedData->GetPedCorrCoeff1(jj);
292   }
293
294   Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
295   Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
296   Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
297   Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
298   Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
299   Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
300   for(Int_t ich=0; ich<5; ich++){
301     adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
302     adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
303     adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
304     adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
305     if(ich<2){
306       adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
307       pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
308     }
309   }
310   
311   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10]; 
312   Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2]; 
313   for(Int_t i=0; i<10; i++){
314      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
315      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
316   }  
317
318   Bool_t isScalerOn=kFALSE;
319   Int_t jsc=0;
320   UInt_t scalerData[32];        
321   for(Int_t k=0; k<32; k++) scalerData[k]=0;
322   
323   Int_t  evQualityBlock[4] = {1,0,0,0};
324   Int_t  triggerBlock[4] = {0,0,0,0};
325   Int_t  chBlock[3] = {0,0,0};
326   UInt_t puBits=0;
327
328   //fNRun = (Int_t) rawReader->GetRunNumber();
329   Int_t kFirstADCGeo=0, kLastADCGeo=3, kScalerGeo=8, kPUGeo=29;
330   //Int_t kTrigScales=30, kTrigHistory=31;
331
332   // loop over raw data
333   rawReader->Reset();
334   AliZDCRawStream rawData(rawReader);
335   while(rawData.Next()){
336    
337    // ***************************** Reading ADCs
338    if((rawData.GetADCModule()>=kFirstADCGeo) && (rawData.GetADCModule()<=kLastADCGeo)){    
339     //printf(" **** Reading ADC raw data from module %d **** \n",rawData.GetADCModule());
340     //
341     if((rawData.IsADCDataWord()) && (rawData.GetNChannelsOn()<48))    chBlock[0] = kTRUE;
342     if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE))  chBlock[1] = kTRUE;
343     if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)) chBlock[2] = kTRUE;
344     if((rawData.IsADCDataWord()) && (rawData.IsADCEventGood() == kTRUE)) evQualityBlock[0] = kTRUE;
345     
346     if((rawData.IsADCDataWord()) && (rawData.IsUnderflow()==kFALSE) 
347         && (rawData.IsOverflow()==kFALSE) && (rawData.IsADCEventGood()==kTRUE)){
348      
349       Int_t adcMod = rawData.GetADCModule();
350       Int_t det = rawData.GetSector(0);
351       Int_t quad = rawData.GetSector(1);
352       Int_t gain = rawData.GetADCGain();
353       Int_t pedindex=0;
354       //
355       // Mean pedestal value subtraction -------------------------------------------------------
356       if(fPedSubMode == 0){
357        // Not interested in o.o.t. signals (ADC modules 2, 3)
358        if(adcMod == 2 || adcMod == 3) continue;
359        //
360        if(quad != 5){ // ZDCs (not reference PTMs)
361         if(det == 1){    
362           pedindex = quad;
363           if(gain == 0) tZN1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
364           else tZN1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
365         }
366         else if(det == 2){ 
367           pedindex = quad+5;
368           if(gain == 0) tZP1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
369           else tZP1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
370         }
371         else if(det == 3){ 
372           pedindex = quad+9;
373           if(quad==1){     
374             if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
375             else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
376           }
377           else if(quad==2){ 
378             if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
379             else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
380           }
381         }
382         else if(det == 4){       
383           pedindex = quad+12;
384           if(gain == 0) tZN2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
385           else tZN2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
386         }
387         else if(det == 5){
388           pedindex = quad+17;
389           if(gain == 0) tZP2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
390           else tZP2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
391         }
392        }
393        else{ // reference PM
394          pedindex = (det-1)/3 + 22;
395          if(det == 1){
396            if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
397          else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
398          }
399          else if(det == 4){
400            if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
401            else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
402          }
403        }
404        // Ch. debug
405        /*if(gain==0){
406          printf(" AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f", 
407            det,quad,gain, pedindex, meanPed[pedindex]);
408          printf("   RawADC %d ADCCorr %1.0f\n", 
409            rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);
410        }*/ 
411       }// mean pedestal subtraction
412       // Pedestal subtraction from correlation ------------------------------------------------
413       else if(fPedSubMode == 1){
414        // In time signals
415        if(adcMod==0 || adcMod==1){
416          if(quad != 5){ // signals from ZDCs
417            if(det == 1){
418              if(gain==0) adcZN1[quad] = rawData.GetADCValue();
419              else adcZN1lg[quad] = rawData.GetADCValue();
420            }
421            else if(det == 2){
422              if(gain==0) adcZP1[quad] = rawData.GetADCValue();
423              else adcZP1lg[quad] = rawData.GetADCValue();
424            }
425            else if(det == 3){
426              if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
427              else adcZEMlg[quad-1] = rawData.GetADCValue();
428            }
429            else if(det == 4){
430              if(gain==0) adcZN2[quad] = rawData.GetADCValue();
431              else adcZN2lg[quad] = rawData.GetADCValue();
432            }
433            else if(det == 5){
434              if(gain==0) adcZP2[quad] = rawData.GetADCValue();
435              else adcZP2lg[quad] = rawData.GetADCValue();
436            }
437          }
438          else{ // signals from reference PM
439             if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
440             else pmReflg[quad-1] = rawData.GetADCValue();
441          }
442        }
443        // Out-of-time pedestals
444        else if(adcMod==2 || adcMod==3){
445          if(quad != 5){ // signals from ZDCs
446            if(det == 1){
447              if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
448              else adcZN1ootlg[quad] = rawData.GetADCValue();
449            }
450            else if(det == 2){
451              if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
452              else adcZP1ootlg[quad] = rawData.GetADCValue();
453            }
454            else if(det == 3){
455              if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
456              else adcZEMootlg[quad-1] = rawData.GetADCValue();
457            }
458            else if(det == 4){
459              if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
460              else adcZN2ootlg[quad] = rawData.GetADCValue();
461            }
462            else if(det == 5){
463              if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
464              else adcZP2ootlg[quad] = rawData.GetADCValue();
465            }
466          }
467          else{ // signals from reference PM
468             if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
469             else pmRefootlg[quad-1] = rawData.GetADCValue();
470          }
471        }
472       } // pedestal subtraction from correlation
473       // Ch. debug
474       /*printf("\t AliZDCReconstructor: det %d quad %d res %d -> Ped[%d] = %1.0f\n", 
475         det,quad,gain, pedindex, meanPed[pedindex]);*/
476     }//IsADCDataWord
477    }// ADC DATA
478    // ***************************** Reading Scaler
479    else if(rawData.GetADCModule()==kScalerGeo){
480      if(rawData.IsScHeaderRead()==kTRUE && rawData.IsScEventGood()==kTRUE){
481        isScalerOn = kTRUE;
482        scalerData[jsc] = rawData.GetTriggerCount();
483        jsc++;
484      }
485    }// VME SCALER DATA
486    // ***************************** Reading PU
487    else if(rawData.GetADCModule()==kPUGeo){
488      puBits = rawData.GetDetectorPattern();
489    }
490     // ***************************** Reading trigger history
491    else if(rawData.IstriggerHistoryWord()==kTRUE){
492      triggerBlock[0] = rawData.IsCPTInputEMDTrigger();
493      triggerBlock[1] = rawData.IsCPTInputSemiCentralTrigger();
494      triggerBlock[2] = rawData.IsCPTInputCentralTrigger();
495      triggerBlock[3] = rawData.IsCPTInputMBTrigger();
496    }
497   
498   }//loop on raw data
499   
500   if(fPedSubMode==1){
501     for(Int_t t=0; t<5; t++){
502        tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
503        tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
504        //
505        tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
506        tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
507        //
508        tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
509        tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
510        //
511        tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
512        tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
513        // 0---------0 Ch. debug 0---------0
514 /*       printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
515        printf("\tCorrCoeff0\tCorrCoeff1\n");
516        printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
517        printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
518        printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
519        printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
520        printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
521        printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
522        printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
523        printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
524        
525        printf("ZN1 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
526                 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
527        printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
528                 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
529        //
530        printf("ZP1 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
531                 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
532        printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
533                 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
534        //
535        printf("ZN2 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
536                 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
537        printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
538                 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
539        //
540        printf("ZP2 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
541                 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
542        printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
543                 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
544 */
545     }
546     dZEM1Corr[0] = adcZEM[0]   - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
547     dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
548     dZEM2Corr[0] = adcZEM[1]   - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
549     dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
550     //
551     sPMRef1[0] = pmRef[0]   - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
552     sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
553     sPMRef2[0] = pmRef[0]   - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
554     sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
555   }
556     
557   if(fRecoMode==1) // p-p data
558     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr, 
559       dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, 
560       isScalerOn, scalerData, 
561       evQualityBlock, triggerBlock, chBlock, puBits);
562   else if(fRecoMode==2) // Pb-Pb data
563     ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr, 
564       dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, 
565       isScalerOn, scalerData, 
566       evQualityBlock, triggerBlock, chBlock, puBits);
567 }
568
569 //_____________________________________________________________________________
570 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* corrADCZN1, 
571         Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
572         Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
573         Bool_t isScalerOn, UInt_t* scaler, 
574         Int_t* evQualityBlock, Int_t* triggerBlock, Int_t* chBlock, UInt_t puBits) const
575 {
576   // ****************** Reconstruct one event ******************
577   
578   // CH. debug
579   /*printf("\n*************************************************\n");
580   printf(" ReconstructEventpp -> values after pedestal subtraction:\n");
581   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
582         corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
583   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
584         corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
585   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
586         corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
587   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
588         corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
589   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
590   printf("*************************************************\n");*/
591     
592   // ---------------------- Setting reco flags for ESD
593   UInt_t rFlags[32];
594   for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
595   
596   if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
597   else rFlags[31] = 0x1;
598   //
599   if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
600   if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
601   if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
602
603   if(triggerBlock[0] == 1) rFlags[27] = 0x1;
604   if(triggerBlock[1] == 1) rFlags[26] = 0x1;
605   if(triggerBlock[2] == 1) rFlags[25] = 0x1;
606   if(triggerBlock[3] == 1) rFlags[24] = 0x1;
607   
608   if(chBlock[0] == 1) rFlags[18] = 0x1;
609   if(chBlock[1] == 1) rFlags[17] = 0x1;
610   if(chBlock[2] == 1) rFlags[16] = 0x1;
611   
612   
613   rFlags[13] = puBits & 0x00000020;
614   rFlags[12] = puBits & 0x00000010;
615   rFlags[11] = puBits & 0x00000080;
616   rFlags[10] = puBits & 0x00000040;
617   rFlags[9]  = puBits & 0x00000020;
618   rFlags[8]  = puBits & 0x00000010;
619   
620   if(corrADCZP1[0]>fSignalThreshold)  rFlags[5] = 0x1;
621   if(corrADCZN1[0]>fSignalThreshold)  rFlags[4] = 0x1;
622   if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
623   if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
624   if(corrADCZP2[0]>fSignalThreshold)  rFlags[1] = 0x1;
625   if(corrADCZN2[0]>fSignalThreshold)  rFlags[0] = 0x1;
626
627   UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
628              rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
629              0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
630              0x0 << 19 | rFlags[18] << 18 |  rFlags[17] << 17 |  rFlags[16] << 16 |
631              0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 | 
632              rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
633              0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 | 
634              rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
635   // --------------------------------------------------
636
637   // ******     Retrieving calibration data 
638   // --- Equalization coefficients ---------------------------------------------
639   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
640   for(Int_t ji=0; ji<5; ji++){
641      equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
642      equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji); 
643      equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji); 
644      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji); 
645   }
646   // --- Energy calibration factors ------------------------------------
647   Float_t calibEne[4];
648   // **** Energy calibration coefficient set to 1 
649   // **** (no trivial way to calibrate in p-p runs)
650   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
651   
652   // ******     Equalization of detector responses
653   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
654   for(Int_t gi=0; gi<10; gi++){
655      if(gi<5){
656        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
657        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
658        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
659        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
660      }
661      else{
662        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
663        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
664        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
665        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
666      }
667   }
668   // Ch. debug
669   /*printf("\n ------------- EQUALIZATION -------------\n");
670   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
671         equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
672   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
673         equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
674   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
675         equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
676   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
677         equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
678   printf(" ----------------------------------------\n");*/
679   
680   // ******     Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
681   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
682   for(Int_t gi=0; gi<5; gi++){
683        calibSumZN1[0] += equalTowZN1[gi];
684        calibSumZP1[0] += equalTowZP1[gi];
685        calibSumZN2[0] += equalTowZN2[gi];
686        calibSumZP2[0] += equalTowZP2[gi];
687        //
688        calibSumZN1[1] += equalTowZN1[gi+5];
689        calibSumZP1[1] += equalTowZP1[gi+5];
690        calibSumZN2[1] += equalTowZN2[gi+5];
691        calibSumZP2[1] += equalTowZP2[gi+5];
692   }
693   // High gain chain
694   calibSumZN1[0] = calibSumZN1[0]*calibEne[0];
695   calibSumZP1[0] = calibSumZP1[0]*calibEne[1];
696   calibSumZN2[0] = calibSumZN2[0]*calibEne[2];
697   calibSumZP2[0] = calibSumZP2[0]*calibEne[3];
698   // Low gain chain
699   calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
700   calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
701   calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
702   calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
703   
704   // ******     Energy calibration of detector responses
705   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
706   for(Int_t gi=0; gi<5; gi++){
707      // High gain chain
708      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
709      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
710      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
711      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
712      // Low gain chain
713      calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
714      calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
715      calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
716      calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
717   }
718   //
719   Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
720   calibZEM1[0] = corrADCZEM1[0]*calibEne[5];
721   calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
722   calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
723   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
724   for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
725   // Ch. debug
726   /*printf("\n ------------- CALIBRATION -------------\n");
727   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
728         calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
729   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
730         calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
731   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
732         calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
733   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
734         calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
735   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
736   printf(" ----------------------------------------\n");*/
737   
738   //  ******    No. of spectator and participants nucleons
739   //  Variables calculated to comply with ESD structure
740   //  *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
741   Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
742   Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
743   Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
744   Double_t impPar=0., impPar1=0., impPar2=0.;
745   
746   // create the output tree
747   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2, 
748                    calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2, 
749                    calibZEM1, calibZEM2, sPMRef1, sPMRef2,
750                    nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight, 
751                    nGenSpec, nGenSpecLeft, nGenSpecRight, 
752                    nPart, nPartTotLeft, nPartTotRight, 
753                    impPar, impPar1, impPar2,
754                    recoFlag, isScalerOn, scaler);
755                   
756   const Int_t kBufferSize = 4000;
757   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
758   // write the output tree
759   clustersTree->Fill();
760 }
761
762 //_____________________________________________________________________________
763 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree, 
764         Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
765         Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
766         Bool_t isScalerOn, UInt_t* scaler, 
767         Int_t* evQualityBlock, Int_t* triggerBlock, Int_t* chBlock, UInt_t puBits) const
768 {
769   // ****************** Reconstruct one event ******************
770   // ---------------------- Setting reco flags for ESD
771   UInt_t rFlags[32];
772   for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
773   
774   if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
775   else rFlags[31] = 0x1;
776   //
777   if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
778   if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
779   if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
780
781   if(triggerBlock[0] == 1) rFlags[27] = 0x1;
782   if(triggerBlock[1] == 1) rFlags[26] = 0x1;
783   if(triggerBlock[2] == 1) rFlags[25] = 0x1;
784   if(triggerBlock[3] == 1) rFlags[24] = 0x1;
785   
786   if(chBlock[0] == 1) rFlags[18] = 0x1;
787   if(chBlock[1] == 1) rFlags[17] = 0x1;
788   if(chBlock[2] == 1) rFlags[16] = 0x1;
789   
790   rFlags[13] = puBits & 0x00000020;
791   rFlags[12] = puBits & 0x00000010;
792   rFlags[11] = puBits & 0x00000080;
793   rFlags[10] = puBits & 0x00000040;
794   rFlags[9]  = puBits & 0x00000020;
795   rFlags[8]  = puBits & 0x00000010;  
796   
797   if(corrADCZP1[0]>fSignalThreshold)  rFlags[5] = 0x1;
798   if(corrADCZN1[0]>fSignalThreshold)  rFlags[4] = 0x1;
799   if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
800   if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
801   if(corrADCZP2[0]>fSignalThreshold)  rFlags[1] = 0x1;
802   if(corrADCZN2[0]>fSignalThreshold)  rFlags[0] = 0x1;
803
804   UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
805              rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
806              0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
807              0x0 << 19 | rFlags[18] << 18 |  rFlags[17] << 17 |  rFlags[16] << 16 |
808              0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 | 
809              rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
810              0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 | 
811              rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
812   // --------------------------------------------------
813   
814
815   // ******     Retrieving calibration data 
816   // --- Equalization coefficients ---------------------------------------------
817   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
818   for(Int_t ji=0; ji<5; ji++){
819      equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
820      equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji); 
821      equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji); 
822      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji); 
823   }
824   // --- Energy calibration factors ------------------------------------
825   Float_t valFromOCDB[6], calibEne[6];
826   for(Int_t ij=0; ij<6; ij++){
827     valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
828     if(ij<4){
829       if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
830       else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
831     } 
832     else calibEne[ij] = valFromOCDB[ij];
833   }
834   
835   // ******     Equalization of detector responses
836   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
837   for(Int_t gi=0; gi<10; gi++){
838      if(gi<5){
839        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
840        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
841        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
842        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
843      }
844      else{
845        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
846        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
847        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
848        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
849      }
850   }
851   
852   // ******     Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
853   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
854   for(Int_t gi=0; gi<5; gi++){
855        calibSumZN1[0] += equalTowZN1[gi];
856        calibSumZP1[0] += equalTowZP1[gi];
857        calibSumZN2[0] += equalTowZN2[gi];
858        calibSumZP2[0] += equalTowZP2[gi];
859        //
860        calibSumZN1[1] += equalTowZN1[gi+5];
861        calibSumZP1[1] += equalTowZP1[gi+5];
862        calibSumZN2[1] += equalTowZN2[gi+5];
863        calibSumZP2[1] += equalTowZP2[gi+5];
864   }
865   // High gain chain
866   calibSumZN1[0] = calibSumZN1[0]*calibEne[0]*8.;
867   calibSumZP1[0] = calibSumZP1[0]*calibEne[1]*8.;
868   calibSumZN2[0] = calibSumZN2[0]*calibEne[2]*8.;
869   calibSumZP2[0] = calibSumZP2[0]*calibEne[3]*8.;
870   // Low gain chain
871   calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
872   calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
873   calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
874   calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
875   //
876   Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
877   calibZEM1[0] = corrADCZEM1[0]*calibEne[5]*8.;
878   calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
879   calibZEM2[0] = corrADCZEM2[0]*calibEne[5]*8.;
880   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
881   for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
882     
883   // ******     Energy calibration of detector responses
884   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
885   for(Int_t gi=0; gi<5; gi++){
886      // High gain chain
887      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]*8.;
888      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]*8.;
889      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]*8.;
890      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]*8.;
891      // Low gain chain
892      calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
893      calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
894      calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
895      calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
896   }
897   
898   //  ******    Number of detected spectator nucleons
899   Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
900   if(fBeamEnergy!=0){
901     nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
902     nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
903     nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
904     nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
905   }
906   else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
907   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
908     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft, 
909     nDetSpecNRight, nDetSpecPRight);*/
910   
911   Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
912   Int_t nPart=0, nPartA=0, nPartC=0;
913   Double_t b=0., bA=0., bC=0.;
914   
915   if(fIsCalibrationMB == kFALSE){
916     // ******   Reconstruction parameters ------------------ 
917     // Ch. debug
918     //fRecoParam->Print("");
919     //
920  
921     if (!fRecoParam) fRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam()); 
922  
923     TH2F *hZDCvsZEM  = fRecoParam->GethZDCvsZEM();
924     TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
925     TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
926     TH1D *hNpartDist = fRecoParam->GethNpartDist();
927     TH1D *hbDist = fRecoParam->GethbDist();    
928     Float_t ClkCenter = fRecoParam->GetClkCenter();
929     //
930     Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
931     Double_t origin = xHighEdge*ClkCenter;
932     // Ch. debug
933     //printf("\n\n  xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
934     //
935     // ====> Summed ZDC info (sideA+side C)
936     TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
937     Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
938     Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
939     line->SetParameter(0, y/(x-origin));
940     line->SetParameter(1, -origin*y/(x-origin));
941     // Ch. debug
942     //printf("  ***************** Summed ZDC info (sideA+side C) \n");
943     //printf("  E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f   ", x, y,y/(x-origin),-origin*y/(x-origin));
944     //
945     Double_t countPerc=0;
946     Double_t xBinCenter=0, yBinCenter=0;
947     for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
948       for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
949          xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
950          yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
951          //
952          if(line->GetParameter(0)>0){
953            if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
954              countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
955              // Ch. debug
956              /*printf(" xBinCenter  %1.3f, yBinCenter %1.0f,  countPerc %1.0f\n", 
957                 xBinCenter, yBinCenter, countPerc);*/
958            }
959          }
960          else{
961            if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
962              countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
963              // Ch. debug
964              /*printf(" xBinCenter  %1.3f, yBinCenter %1.0f,  countPerc %1.0f\n", 
965                 xBinCenter, yBinCenter, countPerc);*/
966            }
967          }
968       }
969     }
970     //
971     Double_t xSecPerc = 0.;
972     if(hZDCvsZEM->GetEntries()!=0){ 
973       xSecPerc = countPerc/hZDCvsZEM->GetEntries();
974     }
975     else{
976       AliWarning("  Histogram hZDCvsZEM from OCDB has no entries!!!");
977     }
978     // Ch. debug
979     //printf("  xSecPerc %1.4f  \n", xSecPerc);
980
981     // ====> side C
982     TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
983     Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
984     lineC->SetParameter(0, yC/(x-origin));
985     lineC->SetParameter(1, -origin*yC/(x-origin));
986     // Ch. debug
987     //printf("  ***************** Side C \n");
988     //printf("  E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f   ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
989     //
990     Double_t countPercC=0;
991     Double_t xBinCenterC=0, yBinCenterC=0;
992     for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
993       for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
994          xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
995          yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
996          if(lineC->GetParameter(0)>0){
997            if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
998              countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
999            }
1000          }
1001          else{
1002            if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1003              countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1004            }
1005          }
1006       }
1007     }
1008     //
1009     Double_t xSecPercC = 0.;
1010     if(hZDCCvsZEM->GetEntries()!=0){ 
1011       xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
1012     }
1013     else{
1014       AliWarning("  Histogram hZDCCvsZEM from OCDB has no entries!!!");
1015     }
1016     // Ch. debug
1017     //printf("  xSecPercC %1.4f  \n", xSecPercC);
1018     
1019     // ====> side A
1020     TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
1021     Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
1022     lineA->SetParameter(0, yA/(x-origin));
1023     lineA->SetParameter(1, -origin*yA/(x-origin));
1024     //
1025     // Ch. debug
1026     //printf("  ***************** Side A \n");
1027     //printf("  E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f   ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
1028     //
1029     Double_t countPercA=0;
1030     Double_t xBinCenterA=0, yBinCenterA=0;
1031     for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
1032       for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
1033          xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
1034          yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
1035          if(lineA->GetParameter(0)>0){
1036            if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1037              countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1038            }
1039          }
1040          else{
1041            if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1042              countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1043            }
1044          }
1045       }
1046     }
1047     //
1048     Double_t xSecPercA = 0.;
1049     if(hZDCAvsZEM->GetEntries()!=0){ 
1050       xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
1051     }
1052     else{
1053       AliWarning("  Histogram hZDCAvsZEM from OCDB has no entries!!!");
1054     }
1055     // Ch. debug
1056     //printf("  xSecPercA %1.4f  \n", xSecPercA);
1057     
1058     //  ******    Number of participants (from E_ZDC vs. E_ZEM correlation)
1059     Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1060     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1061       nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1062       if((1.-nPartFrac) < xSecPerc){
1063         nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1064         // Ch. debug
1065         //printf("  ***************** Summed ZDC info (sideA+side C) \n");
1066         //printf("  nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1067         break;
1068       }
1069     }
1070     if(nPart<0) nPart=0;
1071     //
1072     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1073       nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1074       if((1.-nPartFracC) < xSecPercC){
1075         nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1076         // Ch. debug
1077         //printf("  ***************** Side C \n");
1078         //printf("  nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1079         break;
1080     }
1081     }
1082     if(nPartC<0) nPartC=0;
1083     //
1084     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1085       nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1086       if((1.-nPartFracA) < xSecPercA){
1087         nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1088         // Ch. debug
1089         //printf("  ***************** Side A \n");
1090         //printf("  nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1091         break;
1092       }
1093     }
1094     if(nPartA<0) nPartA=0;
1095     
1096     //  ******    Impact parameter (from E_ZDC vs. E_ZEM correlation)
1097     Double_t bFrac=0., bFracC=0., bFracA=0.;
1098     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1099       bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1100       if((1.-bFrac) < xSecPerc){
1101         b = hbDist->GetBinLowEdge(ibbin);
1102         break;
1103       }
1104     }
1105     //
1106     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1107       bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1108       if((1.-bFracC) < xSecPercC){
1109         bC = hbDist->GetBinLowEdge(ibbin);
1110         break;
1111       }
1112     }
1113     //
1114     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1115       bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1116       if((1.-bFracA) < xSecPercA){
1117         bA = hbDist->GetBinLowEdge(ibbin);
1118         break;
1119       }
1120     }
1121
1122     //  ******  Number of spectator nucleons 
1123     nGenSpec = 416 - nPart;
1124     nGenSpecC = 416 - nPartC;
1125     nGenSpecA = 416 - nPartA;
1126     if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1127     if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1128     if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1129   
1130     //  Ch. debug
1131     /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1132         "  calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n", 
1133         calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1134     printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n", 
1135         nGenSpecLeft, nGenSpecRight);
1136     printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1137     */
1138   
1139     delete lineC;  delete lineA;
1140
1141   } // ONLY IF fIsCalibrationMB==kFALSE
1142   
1143   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2, 
1144                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2, 
1145                   calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1146                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight, 
1147                   nGenSpec, nGenSpecA, nGenSpecC, 
1148                   nPart, nPartA, nPartC, b, bA, bC,
1149                   recoFlag, isScalerOn, scaler);
1150                     
1151   const Int_t kBufferSize = 4000;
1152   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1153   reco->Print("");
1154   // write the output tree
1155   clustersTree->Fill();
1156 }
1157
1158 //_____________________________________________________________________________
1159 void AliZDCReconstructor::BuildRecoParam(Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1160 {
1161   // Calculate RecoParam object for Pb-Pb data
1162   (fRecoParam->GethZDCvsZEM())->Fill(ZDCC+ZDCA, ZEM);
1163   (fRecoParam->GethZDCCvsZEM())->Fill(ZDCC, ZEM);
1164   (fRecoParam->GethZDCAvsZEM())->Fill(ZDCA, ZEM);
1165  
1166 }
1167
1168 //_____________________________________________________________________________
1169 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1170 {
1171   // fill energies and number of participants to the ESD
1172
1173   AliZDCReco reco;
1174   AliZDCReco* preco = &reco;
1175   clustersTree->SetBranchAddress("ZDC", &preco);
1176
1177   clustersTree->GetEntry(0);
1178   //
1179   AliESDZDC * esdzdc = esd->GetESDZDC();
1180   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1181   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1182   for(Int_t i=0; i<5; i++){
1183      tZN1Ene[i] = reco.GetZN1HREnTow(i);
1184      tZN2Ene[i] = reco.GetZN2HREnTow(i);
1185      tZP1Ene[i] = reco.GetZP1HREnTow(i);
1186      tZP2Ene[i] = reco.GetZP2HREnTow(i);
1187      //
1188      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1189      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1190      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1191      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1192   }
1193   //
1194   esdzdc->SetZN1TowerEnergy(tZN1Ene);
1195   esdzdc->SetZN2TowerEnergy(tZN2Ene);
1196   esdzdc->SetZP1TowerEnergy(tZP1Ene);
1197   esdzdc->SetZP2TowerEnergy(tZP2Ene);
1198   //
1199   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1200   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1201   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1202   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1203   // 
1204   Int_t nPart  = reco.GetNParticipants();
1205   Int_t nPartA = reco.GetNPartSideA();
1206   Int_t nPartC = reco.GetNPartSideC();
1207   Double_t b  = reco.GetImpParameter();
1208   Double_t bA = reco.GetImpParSideA();
1209   Double_t bC = reco.GetImpParSideC();
1210   UInt_t recoFlag = reco.GetRecoFlag();
1211   //
1212   esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(), 
1213               reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(), 
1214               nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1215   
1216   // Writing ZDC scaler for cross section calculation
1217   // ONLY IF the scaler has been read during the event
1218   if(reco.IsScalerOn()==kTRUE){
1219     UInt_t counts[32];
1220     for(Int_t jk=0; jk<32; jk++) counts[jk] = reco.GetZDCScaler(jk);
1221     esd->SetZDCScaler(counts);
1222   }
1223 }
1224
1225 //_____________________________________________________________________________
1226 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri) 
1227 {
1228   // Setting the storage
1229
1230   Bool_t deleteManager = kFALSE;
1231   
1232   AliCDBManager *manager = AliCDBManager::Instance();
1233   AliCDBStorage *defstorage = manager->GetDefaultStorage();
1234   
1235   if(!defstorage || !(defstorage->Contains("ZDC"))){ 
1236      AliWarning("No default storage set or default storage doesn't contain ZDC!");
1237      manager->SetDefaultStorage(uri);
1238      deleteManager = kTRUE;
1239   }
1240  
1241   AliCDBStorage *storage = manager->GetDefaultStorage();
1242
1243   if(deleteManager){
1244     AliCDBManager::Instance()->UnsetDefaultStorage();
1245     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
1246   }
1247
1248   return storage; 
1249 }
1250
1251 //_____________________________________________________________________________
1252 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1253 {
1254
1255   // Getting pedestal calibration object for ZDC set
1256
1257   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1258   if(!entry) AliFatal("No calibration data loaded!");  
1259
1260   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
1261   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1262
1263   return calibdata;
1264 }
1265
1266 //_____________________________________________________________________________
1267 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1268 {
1269
1270   // Getting energy and equalization calibration object for ZDC set
1271
1272   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1273   if(!entry) AliFatal("No calibration data loaded!");  
1274
1275   AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1276   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1277
1278   return calibdata;
1279 }
1280
1281 //_____________________________________________________________________________
1282 AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1283 {
1284
1285   // Getting energy and equalization calibration object for ZDC set
1286
1287   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1288   if(!entry) AliFatal("No calibration data loaded!");  
1289
1290   AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1291   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1292
1293   return calibdata;
1294 }
1295
1296 //_____________________________________________________________________________
1297 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1298 {
1299
1300   // Getting reconstruction parameters from OCDB
1301
1302   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1303   if(!entry) AliFatal("No RecoParam data found in OCDB!");  
1304   
1305   AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1306   if(!param)  AliFatal("No RecoParam object in OCDB entry!");
1307   
1308   return param;
1309
1310 }
1311
1312 //_____________________________________________________________________________
1313 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1314 {
1315
1316   // Getting reconstruction parameters from OCDB
1317
1318   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1319   if(!entry) AliFatal("No RecoParam data found in OCDB!");  
1320   
1321   AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1322   if(!param)  AliFatal("No RecoParam object in OCDB entry!");
1323   
1324   return param;
1325
1326 }