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