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