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