ZDC/AliZDCReconstructor.cxx

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