ZDC/AliZDCReconstructor.cxx

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