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