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