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