]> git.uio.no Git - u/mrichter/AliRoot.git/blob - ZDC/AliZDCReconstructor.cxx
08e7add4f6af85a531314472e1537fbd506a7e01
[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     Int_t nPart=0, nPartC=0, nPartA=0;
1060     Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1061     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1062       nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1063       if((1.-nPartFrac) < xSecPerc){
1064         nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1065         // Ch. debug
1066         //printf("  ***************** Summed ZDC info (sideA+side C) \n");
1067         //printf("  nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1068         break;
1069       }
1070     }
1071     if(nPart<0) nPart=0;
1072     //
1073     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1074       nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1075       if((1.-nPartFracC) < xSecPercC){
1076         nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1077         // Ch. debug
1078         //printf("  ***************** Side C \n");
1079         //printf("  nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1080         break;
1081     }
1082     }
1083     if(nPartC<0) nPartC=0;
1084     //
1085     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1086       nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1087       if((1.-nPartFracA) < xSecPercA){
1088         nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1089         // Ch. debug
1090         //printf("  ***************** Side A \n");
1091         //printf("  nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1092         break;
1093       }
1094     }
1095     if(nPartA<0) nPartA=0;
1096     
1097     //  ******    Impact parameter (from E_ZDC vs. E_ZEM correlation)
1098     Float_t b=0, bC=0, bA=0;
1099     Double_t bFrac=0., bFracC=0., bFracA=0.;
1100     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1101       bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1102       if((1.-bFrac) < xSecPerc){
1103         b = hbDist->GetBinLowEdge(ibbin);
1104         break;
1105       }
1106     }
1107     //
1108     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1109       bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1110       if((1.-bFracC) < xSecPercC){
1111         bC = hbDist->GetBinLowEdge(ibbin);
1112         break;
1113       }
1114     }
1115     //
1116     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1117       bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1118       if((1.-bFracA) < xSecPercA){
1119         bA = hbDist->GetBinLowEdge(ibbin);
1120         break;
1121       }
1122     }
1123
1124     //  ******  Number of spectator nucleons 
1125     Int_t nGenSpec=0, nGenSpecC=0, nGenSpecA=0;
1126     //
1127     nGenSpec = 416 - nPart;
1128     nGenSpecC = 416 - nPartC;
1129     nGenSpecA = 416 - nPartA;
1130     if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1131     if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1132     if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1133   
1134     //  Ch. debug
1135     /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1136         "  calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n", 
1137         calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1138     printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n", 
1139         nGenSpecLeft, nGenSpecRight);
1140     printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1141     */
1142   
1143     delete lineC;  delete lineA;
1144
1145   } // ONLY IF fIsCalibrationMB==kFALSE
1146   
1147   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2, 
1148                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2, 
1149                   calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1150                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight, 
1151                   nGenSpec, nGenSpecA, nGenSpecC, 
1152                   nPart, nPartA, nPartC, b, bA, bC,
1153                   recoFlag, isScalerOn, scaler);
1154                     
1155   const Int_t kBufferSize = 4000;
1156   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1157   reco->Print("");
1158   // write the output tree
1159   clustersTree->Fill();
1160 }
1161
1162 //_____________________________________________________________________________
1163 void AliZDCReconstructor::BuildRecoParam(Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1164 {
1165   // Calculate RecoParam object for Pb-Pb data
1166   (fRecoParam->GethZDCvsZEM())->Fill(ZDCC+ZDCA, ZEM);
1167   (fRecoParam->GethZDCCvsZEM())->Fill(ZDCC, ZEM);
1168   (fRecoParam->GethZDCAvsZEM())->Fill(ZDCA, ZEM);
1169  
1170 }
1171
1172 //_____________________________________________________________________________
1173 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1174 {
1175   // fill energies and number of participants to the ESD
1176
1177   AliZDCReco reco;
1178   AliZDCReco* preco = &reco;
1179   clustersTree->SetBranchAddress("ZDC", &preco);
1180
1181   clustersTree->GetEntry(0);
1182   //
1183   AliESDZDC * esdzdc = esd->GetESDZDC();
1184   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1185   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1186   for(Int_t i=0; i<5; i++){
1187      tZN1Ene[i] = reco.GetZN1HREnTow(i);
1188      tZN2Ene[i] = reco.GetZN2HREnTow(i);
1189      tZP1Ene[i] = reco.GetZP1HREnTow(i);
1190      tZP2Ene[i] = reco.GetZP2HREnTow(i);
1191      //
1192      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1193      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1194      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1195      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1196   }
1197   //
1198   esdzdc->SetZN1TowerEnergy(tZN1Ene);
1199   esdzdc->SetZN2TowerEnergy(tZN2Ene);
1200   esdzdc->SetZP1TowerEnergy(tZP1Ene);
1201   esdzdc->SetZP2TowerEnergy(tZP2Ene);
1202   //
1203   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1204   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1205   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1206   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1207   // 
1208   Int_t nPart  = reco.GetNParticipants();
1209   Int_t nPartA = reco.GetNPartSideA();
1210   Int_t nPartC = reco.GetNPartSideC();
1211   Double_t b  = reco.GetImpParameter();
1212   Double_t bA = reco.GetImpParSideA();
1213   Double_t bC = reco.GetImpParSideC();
1214   UInt_t recoFlag = reco.GetRecoFlag();
1215   //
1216   esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(), 
1217               reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(), 
1218               nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1219   
1220   // Writing ZDC scaler for cross section calculation
1221   // ONLY IF the scaler has been read during the event
1222   if(reco.IsScalerOn()==kTRUE){
1223     UInt_t counts[32];
1224     for(Int_t jk=0; jk<32; jk++) counts[jk] = reco.GetZDCScaler(jk);
1225     esd->SetZDCScaler(counts);
1226   }
1227 }
1228
1229 //_____________________________________________________________________________
1230 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri) 
1231 {
1232   // Setting the storage
1233
1234   Bool_t deleteManager = kFALSE;
1235   
1236   AliCDBManager *manager = AliCDBManager::Instance();
1237   AliCDBStorage *defstorage = manager->GetDefaultStorage();
1238   
1239   if(!defstorage || !(defstorage->Contains("ZDC"))){ 
1240      AliWarning("No default storage set or default storage doesn't contain ZDC!");
1241      manager->SetDefaultStorage(uri);
1242      deleteManager = kTRUE;
1243   }
1244  
1245   AliCDBStorage *storage = manager->GetDefaultStorage();
1246
1247   if(deleteManager){
1248     AliCDBManager::Instance()->UnsetDefaultStorage();
1249     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
1250   }
1251
1252   return storage; 
1253 }
1254
1255 //_____________________________________________________________________________
1256 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1257 {
1258
1259   // Getting pedestal calibration object for ZDC set
1260
1261   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1262   if(!entry) AliFatal("No calibration data loaded!");  
1263
1264   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
1265   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1266
1267   return calibdata;
1268 }
1269
1270 //_____________________________________________________________________________
1271 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1272 {
1273
1274   // Getting energy and equalization calibration object for ZDC set
1275
1276   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1277   if(!entry) AliFatal("No calibration data loaded!");  
1278
1279   AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1280   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1281
1282   return calibdata;
1283 }
1284
1285 //_____________________________________________________________________________
1286 AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1287 {
1288
1289   // Getting energy and equalization calibration object for ZDC set
1290
1291   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1292   if(!entry) AliFatal("No calibration data loaded!");  
1293
1294   AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1295   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1296
1297   return calibdata;
1298 }
1299
1300 //_____________________________________________________________________________
1301 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1302 {
1303
1304   // Getting reconstruction parameters from OCDB
1305
1306   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1307   if(!entry) AliFatal("No RecoParam data found in OCDB!");  
1308   
1309   AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1310   if(!param)  AliFatal("No RecoParam object in OCDB entry!");
1311   
1312   return param;
1313
1314 }
1315
1316 //_____________________________________________________________________________
1317 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1318 {
1319
1320   // Getting reconstruction parameters from OCDB
1321
1322   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1323   if(!entry) AliFatal("No RecoParam data found in OCDB!");  
1324   
1325   AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1326   if(!param)  AliFatal("No RecoParam object in OCDB entry!");
1327   
1328   return param;
1329
1330 }