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==2 && quad==2) || (det==3)) && (quad != 5)){
 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 == 2){
 441            if(adcMod==0 || adcMod==1){
 442              if(gain==0) adcZP2[quad] = rawData.GetADCValue();
 443              else adcZP2lg[quad] = rawData.GetADCValue();
 444            }
 445            else if(adcMod==2 || adcMod==3){
 446              if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
 447              else adcZP2ootlg[quad] = rawData.GetADCValue();
 448            }
 449          }
 450          else if(det == 3){
 451            if(adcMod==0 || adcMod==1){
 452              if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
 453              else adcZEMlg[quad-1] = rawData.GetADCValue();
 454            }
 455            else if(adcMod==2 || adcMod==3){
 456              if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
 457              else adcZEMootlg[quad-1] = rawData.GetADCValue();
 458            }
 459          }
 460        }
 461        // When oot values are read the ADC modules 2, 3 can be skipped!!!
 462        if(adcMod == 2 || adcMod == 3) continue;
 463
 464        // *************************************************************************
 465        if(quad != 5){ // ZDCs (not reference PTMs)
 466         if(det==1 && quad!=0){
 467           pedindex = quad;
 468           if(gain == 0) tZN1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 469           else tZN1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 470         }
 471         else if(det==2 && quad!=2){
 472           pedindex = quad+5;
 473           if(gain == 0) tZP1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 474           else tZP1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 475         }
 476         /*else if(det == 3){
 477           pedindex = quad+9;
 478           if(quad==1){
 479             if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 480             else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 481           }
 482           else if(quad==2){
 483             if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 484             else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 485           }
 486         }*/
 487         else if(det == 4){
 488           pedindex = quad+12;
 489           if(gain == 0) tZN2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 490           else tZN2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 491         }
 492         else if(det == 5){
 493           pedindex = quad+17;
 494           if(gain == 0) tZP2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 495           else tZP2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 496         }
 497        }
 498        else{ // reference PM
 499          pedindex = (det-1)/3 + 22;
 500          if(det == 1){
 501            if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 502          else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 503          }
 504          else if(det == 4){
 505            if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 506            else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 507          }
 508        }
 509        // Ch. debug
 510        /*if(gain==0){
 511          printf(" AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f",
 512            det,quad,gain, pedindex, meanPed[pedindex]);
 513          printf("   RawADC %d ADCCorr %1.0f\n",
 514            rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);
 515        }*/
 516       }// mean pedestal subtraction
 517       // Pedestal subtraction from correlation ------------------------------------------------
 518       else if(fPedSubMode == 1){
 519        // In time signals
 520        if(adcMod==0 || adcMod==1){
 521          if(quad != 5){ // signals from ZDCs
 522            if(det == 1){
 523              if(gain==0) adcZN1[quad] = rawData.GetADCValue();
 524              else adcZN1lg[quad] = rawData.GetADCValue();
 525            }
 526            else if(det == 2){
 527              if(gain==0) adcZP1[quad] = rawData.GetADCValue();
 528              else adcZP1lg[quad] = rawData.GetADCValue();
 529            }
 530            else if(det == 3){
 531              if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
 532              else adcZEMlg[quad-1] = rawData.GetADCValue();
 533            }
 534            else if(det == 4){
 535              if(gain==0) adcZN2[quad] = rawData.GetADCValue();
 536              else adcZN2lg[quad] = rawData.GetADCValue();
 537            }
 538            else if(det == 5){
 539              if(gain==0) adcZP2[quad] = rawData.GetADCValue();
 540              else adcZP2lg[quad] = rawData.GetADCValue();
 541            }
 542          }
 543          else{ // signals from reference PM
 544             if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
 545             else pmReflg[quad-1] = rawData.GetADCValue();
 546          }
 547        }
 548        // Out-of-time pedestals
 549        else if(adcMod==2 || adcMod==3){
 550          if(quad != 5){ // signals from ZDCs
 551            if(det == 1){
 552              if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
 553              else adcZN1ootlg[quad] = rawData.GetADCValue();
 554            }
 555            else if(det == 2){
 556              if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
 557              else adcZP1ootlg[quad] = rawData.GetADCValue();
 558            }
 559            else if(det == 3){
 560              if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
 561              else adcZEMootlg[quad-1] = rawData.GetADCValue();
 562            }
 563            else if(det == 4){
 564              if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
 565              else adcZN2ootlg[quad] = rawData.GetADCValue();
 566            }
 567            else if(det == 5){
 568              if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
 569              else adcZP2ootlg[quad] = rawData.GetADCValue();
 570            }
 571          }
 572          else{ // signals from reference PM
 573             if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
 574             else pmRefootlg[quad-1] = rawData.GetADCValue();
 575          }
 576        }
 577       } // pedestal subtraction from correlation
 578       // Ch. debug
 579       /*printf("\t AliZDCReconstructor: det %d quad %d res %d -> Ped[%d] = %1.0f\n",
 580         det,quad,gain, pedindex, meanPed[pedindex]);*/
 581     }//IsADCDataWord
 582    }// ADC DATA
 583    // ***************************** Reading Scaler
 584    else if(rawData.GetADCModule()==kScalerGeo){
 585      if(rawData.IsScalerWord()==kTRUE && rawData.IsScEventGood()==kTRUE){
 586        isScalerOn = kTRUE;
 587        scalerData[jsc] = rawData.GetTriggerCount();
 588        // Ch. debug
 589        //printf("   Reconstructed VME Scaler: %d %d  ",jsc,scalerData[jsc]);
 590        //
 591        jsc++;
 592      }
 593    }// VME SCALER DATA
 594    // ***************************** Reading ZDC TDC
 595    else if(rawData.GetADCModule()==kZDCTDCGeo && rawData.IsZDCTDCDatum()==kTRUE){
 596        itdc = rawData.GetChannel();
 597        if(itdc==iprevtdc) ihittdc++;
 598        else ihittdc=0;
 599        iprevtdc=itdc;
 600        tdcData[itdc][ihittdc] = rawData.GetZDCTDCDatum();
 601        // Ch. debug
 602        //printf("   Reconstructed TDC[%d, %d] %d  ",itdc, ihittdc, tdcData[itdc][ihittdc]);
 603    }// ZDC TDC DATA
 604    // ***************************** Reading PU
 605    else if(rawData.GetADCModule()==kPUGeo){
 606      puBits = rawData.GetDetectorPattern();
 607    }
 608     // ***************************** Reading trigger history
 609    else if(rawData.IstriggerHistoryWord()==kTRUE){
 610      triggerBlock[0] = rawData.IsCPTInputEMDTrigger();
 611      triggerBlock[1] = rawData.IsCPTInputSemiCentralTrigger();
 612      triggerBlock[2] = rawData.IsCPTInputCentralTrigger();
 613      triggerBlock[3] = rawData.IsCPTInputMBTrigger();
 614    }
 615
 616   }//loop on raw data
 617
 618   if(fPedSubMode==1){
 619     for(Int_t t=0; t<5; t++){
 620        tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
 621        tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
 622        //
 623        tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
 624        tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
 625        //
 626        tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
 627        tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
 628        //
 629        tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
 630        tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
 631     }
 632     dZEM1Corr[0] = adcZEM[0]   - (corrCoeff1[10]*adcZEMoot[0]+corrCoeff0[10]);
 633     dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[10+kNch]*adcZEMootlg[0]+corrCoeff0[10+kNch]);
 634     dZEM2Corr[0] = adcZEM[1]   - (corrCoeff1[11]*adcZEMoot[1]+corrCoeff0[11]);
 635     dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[11+kNch]*adcZEMootlg[1]+corrCoeff0[11+kNch]);
 636     //
 637     sPMRef1[0] = pmRef[0]   - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
 638     sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
 639     sPMRef2[0] = pmRef[0]   - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
 640     sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
 641   }
 642   else{
 643     //  **** Pb-Pb data taking 2010 -> subtracting some ch. from correlation ****
 644     tZN1Corr[0] = adcZN1[0] - (corrCoeff1[0]*adcZN1oot[0]+corrCoeff0[0]);
 645     tZN1Corr[5] = adcZN1lg[0] - (corrCoeff1[kNch]*adcZN1ootlg[0]+corrCoeff0[kNch]);
 646     // Ch. debug
 647     //printf(" adcZN1 %d  adcZN1oot %d tZN1Corr %1.2f \n", adcZN1[0],adcZN1oot[0],tZN1Corr[0]);
 648     //printf(" adcZN1lg %d  adcZN1ootlg %d tZN1Corrlg %1.2f \n", adcZN1lg[0],adcZN1ootlg[0],tZN1Corr[5]);
 649     //
 650     tZP1Corr[2] = adcZP1[2] - (corrCoeff1[2+5]*adcZP1oot[2]+corrCoeff0[2+5]);
 651     tZP1Corr[2+5] = adcZP1lg[2] - (corrCoeff1[2+5+kNch]*adcZP1ootlg[2]+corrCoeff0[2+5+kNch]);
 652     //
 653     dZEM1Corr[0] = adcZEM[0]   - (corrCoeff1[10]*adcZEMoot[0]+corrCoeff0[10]);
 654     dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[10+kNch]*adcZEMootlg[0]+corrCoeff0[10+kNch]);
 655     dZEM2Corr[0] = adcZEM[1]   - (corrCoeff1[11]*adcZEMoot[1]+corrCoeff0[11]);
 656     dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[11+kNch]*adcZEMootlg[1]+corrCoeff0[11+kNch]);
 657     // *************************************************************************
 658   }
 659
 660   if(fRecoMode==1) // p-p data
 661     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 662       dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
 663       isScalerOn, scalerData, tdcData,
 664       evQualityBlock, triggerBlock, chBlock, puBits);
 665   else if(fRecoMode==2) // Pb-Pb data
 666       ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 667       dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
 668       isScalerOn, scalerData,  tdcData,
 669       evQualityBlock, triggerBlock, chBlock, puBits);
 670 }
 671
 672 //_____________________________________________________________________________
 673 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree,
 674         const Float_t* const corrADCZN1, const Float_t* const corrADCZP1,
 675         const Float_t* const corrADCZN2, const Float_t* const corrADCZP2,
 676         const Float_t* const corrADCZEM1, const Float_t* const corrADCZEM2,
 677         Float_t* sPMRef1, Float_t* sPMRef2, Bool_t isScalerOn, UInt_t* scaler,
 678         Int_t tdcData[32][4], const Int_t* const evQualityBlock,
 679         const Int_t* const triggerBlock, const Int_t* const chBlock, UInt_t puBits) const
 680 {
 681   // ****************** Reconstruct one event ******************
 682
 683   // CH. debug
 684   /*printf("\n*************************************************\n");
 685   printf(" ReconstructEventpp -> values after pedestal subtraction:\n");
 686   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 687         corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
 688   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 689         corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
 690   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 691         corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
 692   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 693         corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
 694   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
 695   printf("*************************************************\n");*/
 696
 697   // ---------------------- Setting reco flags for ESD
 698   UInt_t rFlags[32];
 699   for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
 700
 701   if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
 702   else rFlags[31] = 0x1;
 703   //
 704   if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
 705   if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
 706   if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
 707
 708   if(triggerBlock[0] == 1) rFlags[27] = 0x1;
 709   if(triggerBlock[1] == 1) rFlags[26] = 0x1;
 710   if(triggerBlock[2] == 1) rFlags[25] = 0x1;
 711   if(triggerBlock[3] == 1) rFlags[24] = 0x1;
 712
 713   if(chBlock[0] == 1) rFlags[18] = 0x1;
 714   if(chBlock[1] == 1) rFlags[17] = 0x1;
 715   if(chBlock[2] == 1) rFlags[16] = 0x1;
 716
 717
 718   rFlags[13] = puBits & 0x00000020;
 719   rFlags[12] = puBits & 0x00000010;
 720   rFlags[11] = puBits & 0x00000080;
 721   rFlags[10] = puBits & 0x00000040;
 722   rFlags[9]  = puBits & 0x00000020;
 723   rFlags[8]  = puBits & 0x00000010;
 724
 725   if(corrADCZP1[0]>fSignalThreshold)  rFlags[5] = 0x1;
 726   if(corrADCZN1[0]>fSignalThreshold)  rFlags[4] = 0x1;
 727   if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
 728   if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
 729   if(corrADCZP2[0]>fSignalThreshold)  rFlags[1] = 0x1;
 730   if(corrADCZN2[0]>fSignalThreshold)  rFlags[0] = 0x1;
 731
 732   UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
 733              rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
 734              0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
 735              0x0 << 19 | rFlags[18] << 18 |  rFlags[17] << 17 |  rFlags[16] << 16 |
 736              0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
 737              rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
 738              0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
 739              rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
 740   // --------------------------------------------------
 741
 742   // ******     Retrieving calibration data
 743   // --- Equalization coefficients ---------------------------------------------
 744   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 745   for(Int_t ji=0; ji<5; ji++){
 746      equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
 747      equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
 748      equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
 749      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
 750   }
 751   // --- Energy calibration factors ------------------------------------
 752   Float_t calibEne[6];
 753   // **** Energy calibration coefficient set to 1
 754   // **** (no trivial way to calibrate in p-p runs)
 755   for(Int_t ij=0; ij<6; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
 756
 757   // ******     Equalization of detector responses
 758   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 759   for(Int_t gi=0; gi<10; gi++){
 760      if(gi<5){
 761        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
 762        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
 763        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
 764        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
 765      }
 766      else{
 767        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
 768        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
 769        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
 770        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
 771      }
 772   }
 773   // Ch. debug
 774   /*printf("\n ------------- EQUALIZATION -------------\n");
 775   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 776         equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
 777   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 778         equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
 779   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 780         equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
 781   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 782         equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
 783   printf(" ----------------------------------------\n");*/
 784
 785   // ******     Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
 786   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 787   for(Int_t gi=0; gi<5; gi++){
 788        calibSumZN1[0] += equalTowZN1[gi];
 789        calibSumZP1[0] += equalTowZP1[gi];
 790        calibSumZN2[0] += equalTowZN2[gi];
 791        calibSumZP2[0] += equalTowZP2[gi];
 792        //
 793        calibSumZN1[1] += equalTowZN1[gi+5];
 794        calibSumZP1[1] += equalTowZP1[gi+5];
 795        calibSumZN2[1] += equalTowZN2[gi+5];
 796        calibSumZP2[1] += equalTowZP2[gi+5];
 797   }
 798   // High gain chain
 799   calibSumZN1[0] = calibSumZN1[0]*calibEne[0];
 800   calibSumZP1[0] = calibSumZP1[0]*calibEne[1];
 801   calibSumZN2[0] = calibSumZN2[0]*calibEne[2];
 802   calibSumZP2[0] = calibSumZP2[0]*calibEne[3];
 803   // Low gain chain
 804   calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
 805   calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
 806   calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
 807   calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
 808
 809   // ******     Energy calibration of detector responses
 810   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 811   for(Int_t gi=0; gi<5; gi++){
 812      // High gain chain
 813      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
 814      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
 815      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
 816      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
 817      // Low gain chain
 818      calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
 819      calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
 820      calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
 821      calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
 822   }
 823   //
 824   Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
 825   calibZEM1[0] = corrADCZEM1[0]*calibEne[4];
 826   calibZEM1[1] = corrADCZEM1[1]*calibEne[4];
 827   calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
 828   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
 829   for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
 830   // Ch. debug
 831   /*printf("\n ------------- CALIBRATION -------------\n");
 832   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 833         calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
 834   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 835         calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
 836   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 837         calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
 838   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 839         calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
 840   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
 841   printf(" ----------------------------------------\n");*/
 842
 843   //  ******    No. of spectator and participants nucleons
 844   //  Variables calculated to comply with ESD structure
 845   //  *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
 846   Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
 847   Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
 848   Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
 849   Double_t impPar=0., impPar1=0., impPar2=0.;
 850
 851   // create the output tree
 852   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 853                    calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
 854                    calibZEM1, calibZEM2, sPMRef1, sPMRef2,
 855                    nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 856                    nGenSpec, nGenSpecLeft, nGenSpecRight,
 857                    nPart, nPartTotLeft, nPartTotRight,
 858                    impPar, impPar1, impPar2,
 859                    recoFlag, isScalerOn, scaler, tdcData);
 860
 861   const Int_t kBufferSize = 4000;
 862   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
 863   // write the output tree
 864   clustersTree->Fill();
 865   delete reco;
 866 }
 867
 868 //_____________________________________________________________________________
 869 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
 870         const Float_t* const corrADCZN1, const Float_t* const corrADCZP1,
 871         const Float_t* const corrADCZN2, const Float_t* const corrADCZP2,
 872         const Float_t* const corrADCZEM1, const Float_t* const corrADCZEM2,
 873         Float_t* sPMRef1, Float_t* sPMRef2, Bool_t isScalerOn, UInt_t* scaler,
 874         Int_t tdcData[32][4], const Int_t* const evQualityBlock,
 875         const Int_t* const triggerBlock, const Int_t* const chBlock, UInt_t puBits) const
 876 {
 877   // ****************** Reconstruct one event ******************
 878   // ---------------------- Setting reco flags for ESD
 879   UInt_t rFlags[32];
 880   for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
 881
 882   if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
 883   else rFlags[31] = 0x1;
 884   //
 885   if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
 886   if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
 887   if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
 888
 889   if(triggerBlock[0] == 1) rFlags[27] = 0x1;
 890   if(triggerBlock[1] == 1) rFlags[26] = 0x1;
 891   if(triggerBlock[2] == 1) rFlags[25] = 0x1;
 892   if(triggerBlock[3] == 1) rFlags[24] = 0x1;
 893
 894   if(chBlock[0] == 1) rFlags[18] = 0x1;
 895   if(chBlock[1] == 1) rFlags[17] = 0x1;
 896   if(chBlock[2] == 1) rFlags[16] = 0x1;
 897
 898   rFlags[13] = puBits & 0x00000020;
 899   rFlags[12] = puBits & 0x00000010;
 900   rFlags[11] = puBits & 0x00000080;
 901   rFlags[10] = puBits & 0x00000040;
 902   rFlags[9]  = puBits & 0x00000020;
 903   rFlags[8]  = puBits & 0x00000010;
 904
 905   if(corrADCZP1[0]>fSignalThreshold)  rFlags[5] = 0x1;
 906   if(corrADCZN1[0]>fSignalThreshold)  rFlags[4] = 0x1;
 907   if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
 908   if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
 909   if(corrADCZP2[0]>fSignalThreshold)  rFlags[1] = 0x1;
 910   if(corrADCZN2[0]>fSignalThreshold)  rFlags[0] = 0x1;
 911
 912   UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
 913              rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
 914              0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
 915              0x0 << 19 | rFlags[18] << 18 |  rFlags[17] << 17 |  rFlags[16] << 16 |
 916              0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
 917              rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
 918              0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
 919              rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
 920   // --------------------------------------------------
 921
 922
 923   // CH. debug
 924 /*  printf("\n*************************************************\n");
 925   printf(" ReconstructEventPbPb -> values after pedestal subtraction:\n");
 926   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 927         corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
 928   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 929         corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
 930   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 931         corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
 932   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 933         corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
 934   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
 935   printf("*************************************************\n");
 936 */
 937   // ******     Retrieving calibration data
 938   // --- Equalization coefficients ---------------------------------------------
 939   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 940   for(Int_t ji=0; ji<5; ji++){
 941      equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
 942      equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
 943      equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
 944      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
 945   }
 946   // --- Energy calibration factors ------------------------------------
 947   Float_t calibEne[6];
 948   // The energy calibration object already takes into account of E_beam
 949   // -> the value from the OCDB can be directly used (Jul 2010)
 950   for(Int_t ij=0; ij<6; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
 951
 952   // ******     Equalization of detector responses
 953   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 954   for(Int_t gi=0; gi<10; gi++){
 955      if(gi<5){
 956        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
 957        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
 958        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
 959        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
 960      }
 961      else{
 962        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
 963        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
 964        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
 965        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
 966      }
 967   }
 968
 969   // Ch. debug
 970 /*  printf("\n ------------- EQUALIZATION -------------\n");
 971   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 972         equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
 973   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 974         equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
 975   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 976         equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
 977   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 978         equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
 979   printf(" ----------------------------------------\n");
 980 */
 981
 982   // ******     Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
 983   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 984   for(Int_t gi=0; gi<5; gi++){
 985        calibSumZN1[0] += equalTowZN1[gi];
 986        calibSumZP1[0] += equalTowZP1[gi];
 987        calibSumZN2[0] += equalTowZN2[gi];
 988        calibSumZP2[0] += equalTowZP2[gi];
 989        //
 990        calibSumZN1[1] += equalTowZN1[gi+5];
 991        calibSumZP1[1] += equalTowZP1[gi+5];
 992        calibSumZN2[1] += equalTowZN2[gi+5];
 993        calibSumZP2[1] += equalTowZP2[gi+5];
 994   }
 995   //
 996   //fEnCalibData->Print("");
 997
 998   // High gain chain
 999   calibSumZN1[0] = calibSumZN1[0]*calibEne[0]*8.;
1000   calibSumZP1[0] = calibSumZP1[0]*calibEne[1]*8.;
1001   calibSumZN2[0] = calibSumZN2[0]*calibEne[2]*8.;
1002   calibSumZP2[0] = calibSumZP2[0]*calibEne[3]*8.;
1003   // Low gain chain
1004   calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
1005   calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
1006   calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
1007   calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
1008   //
1009   Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
1010   calibZEM1[0] = corrADCZEM1[0]*calibEne[4]*8.;
1011   calibZEM1[1] = corrADCZEM1[1]*calibEne[4];
1012   calibZEM2[0] = corrADCZEM2[0]*calibEne[5]*8.;
1013   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
1014   for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
1015
1016   // ******     Energy calibration of detector responses
1017   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
1018   for(Int_t gi=0; gi<5; gi++){
1019      // High gain chain
1020      calibTowZN1[gi] = equalTowZN1[gi]*2*calibEne[0]*8.;
1021      calibTowZP1[gi] = equalTowZP1[gi]*2*calibEne[1]*8.;
1022      calibTowZN2[gi] = equalTowZN2[gi]*2*calibEne[2]*8.;
1023      calibTowZP2[gi] = equalTowZP2[gi]*2*calibEne[3]*8.;
1024      // Low gain chain
1025      calibTowZN1[gi+5] = equalTowZN1[gi+5]*2*calibEne[0];
1026      calibTowZP1[gi+5] = equalTowZP1[gi+5]*2*calibEne[1];
1027      calibTowZN2[gi+5] = equalTowZN2[gi+5]*2*calibEne[2];
1028      calibTowZP2[gi+5] = equalTowZP2[gi+5]*2*calibEne[3];
1029   }
1030
1031   // Ch. debug
1032 /*  printf("\n ------------- CALIBRATION -------------\n");
1033   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1034         calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
1035   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1036         calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
1037   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1038         calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
1039   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1040         calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
1041   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
1042   printf(" ----------------------------------------\n");
1043 */
1044   //  ******    Number of detected spectator nucleons
1045   Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
1046   if(fBeamEnergy>0.01){
1047     nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
1048     nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
1049     nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
1050     nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
1051   }
1052   else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
1053   /*printf("\n\t AliZDCReconstructor -> fBeamEnergy %1.0f: nDetSpecNsideA %d, nDetSpecPsideA %d,"
1054     " nDetSpecNsideC %d, nDetSpecPsideC %d\n",fBeamEnergy,nDetSpecNLeft, nDetSpecPLeft,
1055     nDetSpecNRight, nDetSpecPRight);*/
1056
1057   Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
1058   Int_t nPart=0, nPartA=0, nPartC=0;
1059   Double_t b=0., bA=0., bC=0.;
1060
1061   if(fIsCalibrationMB == kFALSE){
1062     // ******   Reconstruction parameters ------------------
1063     if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
1064     if(!fgRecoParam){
1065       AliError("  RecoParam object not retrieved correctly: not reconstructing event!!!");
1066       return;
1067     }
1068     TH1D* hNpartDist = fgRecoParam->GethNpartDist();
1069     TH1D*   hbDist = fgRecoParam->GethbDist();
1070     Float_t  fClkCenter = fgRecoParam->GetClkCenter();
1071     if(!hNpartDist || !hbDist){
1072        AliError("Something wrong in Glauber MC histos got from AliZDCREcoParamPbPb: NO EVENT RECO FOR ZDC DATA!!!\n\n");
1073        return;
1074     }
1075
1076     if(!fgMBCalibData) fgMBCalibData = const_cast<AliZDCMBCalib*>(GetMBCalibData());
1077     TH2F *hZDCvsZEM  = fgMBCalibData->GethZDCvsZEM();
1078     TH2F *hZDCCvsZEM = fgMBCalibData->GethZDCCvsZEM();
1079     TH2F *hZDCAvsZEM = fgMBCalibData->GethZDCAvsZEM();
1080     //
1081     Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
1082     Double_t origin = xHighEdge*fClkCenter;
1083     // Ch. debug
1084     //printf("\n\n  xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
1085     //
1086     // ====> Summed ZDC info (sideA+side C)
1087     TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
1088     Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
1089     Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
1090     line->SetParameter(0, y/(x-origin));
1091     line->SetParameter(1, -origin*y/(x-origin));
1092     // Ch. debug
1093     //printf("  ***************** Summed ZDC info (sideA+side C) \n");
1094     //printf("  E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f   ", x, y,y/(x-origin),-origin*y/(x-origin));
1095     //
1096     Double_t countPerc=0;
1097     Double_t xBinCenter=0, yBinCenter=0;
1098     for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
1099       for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
1100          xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
1101          yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
1102          //
1103          if(line->GetParameter(0)>0){
1104            if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
1105              countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
1106              // Ch. debug
1107              //printf(" xBinCenter  %1.3f, yBinCenter %1.0f,  countPerc %1.0f\n",
1108                 //xBinCenter, yBinCenter, countPerc);
1109            }
1110          }
1111          else{
1112            if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
1113              countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
1114              // Ch. debug
1115              //printf(" xBinCenter  %1.3f, yBinCenter %1.0f,  countPerc %1.0f\n",
1116                 //xBinCenter, yBinCenter, countPerc);
1117            }
1118          }
1119       }
1120     }
1121     //
1122     Double_t xSecPerc = 0.;
1123     if(hZDCvsZEM->GetEntries()!=0){
1124       xSecPerc = countPerc/hZDCvsZEM->GetEntries();
1125     }
1126     else{
1127       AliWarning("  Histogram hZDCvsZEM from OCDB has no entries!!!");
1128     }
1129     // Ch. debug
1130     //printf("  xSecPerc %1.4f  \n", xSecPerc);
1131
1132     // ====> side C
1133     TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
1134     Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
1135     lineC->SetParameter(0, yC/(x-origin));
1136     lineC->SetParameter(1, -origin*yC/(x-origin));
1137     // Ch. debug
1138     //printf("  ***************** Side C \n");
1139     //printf("  E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f   ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
1140     //
1141     Double_t countPercC=0;
1142     Double_t xBinCenterC=0, yBinCenterC=0;
1143     for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
1144       for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
1145          xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
1146          yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
1147          if(lineC->GetParameter(0)>0){
1148            if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1149              countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1150            }
1151          }
1152          else{
1153            if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1154              countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1155            }
1156          }
1157       }
1158     }
1159     //
1160     Double_t xSecPercC = 0.;
1161     if(hZDCCvsZEM->GetEntries()!=0){
1162       xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
1163     }
1164     else{
1165       AliWarning("  Histogram hZDCCvsZEM from OCDB has no entries!!!");
1166     }
1167     // Ch. debug
1168     //printf("  xSecPercC %1.4f  \n", xSecPercC);
1169
1170     // ====> side A
1171     TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
1172     Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
1173     lineA->SetParameter(0, yA/(x-origin));
1174     lineA->SetParameter(1, -origin*yA/(x-origin));
1175     //
1176     // Ch. debug
1177     //printf("  ***************** Side A \n");
1178     //printf("  E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f   ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
1179     //
1180     Double_t countPercA=0;
1181     Double_t xBinCenterA=0, yBinCenterA=0;
1182     for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
1183       for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
1184          xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
1185          yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
1186          if(lineA->GetParameter(0)>0){
1187            if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1188              countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1189            }
1190          }
1191          else{
1192            if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1193              countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1194            }
1195          }
1196       }
1197     }
1198     //
1199     Double_t xSecPercA = 0.;
1200     if(hZDCAvsZEM->GetEntries()!=0){
1201       xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
1202     }
1203     else{
1204       AliWarning("  Histogram hZDCAvsZEM from OCDB has no entries!!!");
1205     }
1206     // Ch. debug
1207     //printf("  xSecPercA %1.4f  \n", xSecPercA);
1208
1209     //  ******    Number of participants (from E_ZDC vs. E_ZEM correlation)
1210     Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1211     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1212       nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1213       if((1.-nPartFrac) < xSecPerc){
1214         nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1215         // Ch. debug
1216         //printf("  ***************** Summed ZDC info (sideA+side C) \n");
1217         //printf("  nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1218         break;
1219       }
1220     }
1221     if(nPart<0) nPart=0;
1222     //
1223     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1224       nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1225       if((1.-nPartFracC) < xSecPercC){
1226         nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1227         // Ch. debug
1228         //printf("  ***************** Side C \n");
1229         //printf("  nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1230         break;
1231     }
1232     }
1233     if(nPartC<0) nPartC=0;
1234     //
1235     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1236       nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1237       if((1.-nPartFracA) < xSecPercA){
1238         nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1239         // Ch. debug
1240         //printf("  ***************** Side A \n");
1241         //printf("  nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1242         break;
1243       }
1244     }
1245     if(nPartA<0) nPartA=0;
1246
1247     //  ******    Impact parameter (from E_ZDC vs. E_ZEM correlation)
1248     Double_t bFrac=0., bFracC=0., bFracA=0.;
1249     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1250       bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1251       if(bFrac > xSecPerc){
1252         b = hbDist->GetBinLowEdge(ibbin);
1253         break;
1254       }
1255     }
1256     //
1257     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1258       bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1259       if(bFracC > xSecPercC){
1260         bC = hbDist->GetBinLowEdge(ibbin);
1261         break;
1262       }
1263     }
1264     //
1265     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1266       bFracA += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1267       if(bFracA > xSecPercA){
1268         bA = hbDist->GetBinLowEdge(ibbin);
1269         break;
1270       }
1271     }
1272
1273     //  ******  Number of spectator nucleons
1274     nGenSpec = 416 - nPart;
1275     nGenSpecC = 416 - nPartC;
1276     nGenSpecA = 416 - nPartA;
1277     if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1278     if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1279     if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1280
1281     delete line;
1282     delete lineC;  delete lineA;
1283
1284   } // ONLY IF fIsCalibrationMB==kFALSE
1285
1286   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1287                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1288                   calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1289                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1290                   nGenSpec, nGenSpecA, nGenSpecC,
1291                   nPart, nPartA, nPartC, b, bA, bC,
1292                   recoFlag, isScalerOn, scaler, tdcData);
1293
1294   const Int_t kBufferSize = 4000;
1295   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1296   //reco->Print("");
1297   // write the output tree
1298   clustersTree->Fill();
1299   delete reco;
1300 }
1301
1302
1303 //_____________________________________________________________________________
1304 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1305 {
1306   // fill energies and number of participants to the ESD
1307
1308   AliZDCReco reco;
1309   AliZDCReco* preco = &reco;
1310   clustersTree->SetBranchAddress("ZDC", &preco);
1311   clustersTree->GetEntry(0);
1312   //
1313   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1314   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1315   for(Int_t i=0; i<5; i++){
1316      tZN1Ene[i] = reco.GetZN1HREnTow(i);
1317      tZN2Ene[i] = reco.GetZN2HREnTow(i);
1318      tZP1Ene[i] = reco.GetZP1HREnTow(i);
1319      tZP2Ene[i] = reco.GetZP2HREnTow(i);
1320      //
1321      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1322      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1323      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1324      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1325   }
1326   //
1327   fESDZDC->SetZN1TowerEnergy(tZN1Ene);
1328   fESDZDC->SetZN2TowerEnergy(tZN2Ene);
1329   fESDZDC->SetZP1TowerEnergy(tZP1Ene);
1330   fESDZDC->SetZP2TowerEnergy(tZP2Ene);
1331   //
1332   fESDZDC->SetZN1TowerEnergyLR(tZN1EneLR);
1333   fESDZDC->SetZN2TowerEnergyLR(tZN2EneLR);
1334   fESDZDC->SetZP1TowerEnergyLR(tZP1EneLR);
1335   fESDZDC->SetZP2TowerEnergyLR(tZP2EneLR);
1336   //
1337   Int_t nPart  = reco.GetNParticipants();
1338   Int_t nPartA = reco.GetNPartSideA();
1339   Int_t nPartC = reco.GetNPartSideC();
1340   Double_t b  = reco.GetImpParameter();
1341   Double_t bA = reco.GetImpParSideA();
1342   Double_t bC = reco.GetImpParSideC();
1343   UInt_t recoFlag = reco.GetRecoFlag();
1344
1345   fESDZDC->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(),
1346         reco.GetZEM1HRsignal(), reco.GetZEM2HRsignal(),
1347         reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1348         nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1349
1350   // Writing ZDC scaler for cross section calculation
1351   // ONLY IF the scaler has been read during the event
1352   if(reco.IsScalerOn()==kTRUE){
1353     UInt_t counts[32];
1354     for(Int_t jk=0; jk<32; jk++) counts[jk] = reco.GetZDCScaler(jk);
1355     fESDZDC->SetZDCScaler(counts);
1356   }
1357
1358   // Writing TDC data into ZDC ESDs
1359   Int_t tdcValues[32][4];
1360   Float_t tdcCorrected[32][4];
1361   for(Int_t jk=0; jk<32; jk++){
1362     for(Int_t lk=0; lk<4; lk++){
1363       tdcValues[jk][lk] = reco.GetZDCTDCData(jk, lk);
1364     }
1365   }
1366   // 4/2/2011 -> Subtracting L0 (tdcValues[15]) instead of ADC gate
1367   // we try to keep the TDC oscillations as low as possible!
1368   for(Int_t jk=0; jk<32; jk++){
1369     for(Int_t lk=0; lk<4; lk++){
1370       if(tdcValues[jk][lk]!=0.) tdcCorrected[jk][lk] = 0.025*(tdcValues[jk][lk]-tdcValues[15][0])+fMeanPhase;
1371     }
1372   }
1373   fESDZDC->SetZDCTDCData(tdcValues);
1374   fESDZDC->SetZDCTDCCorrected(tdcCorrected);
1375   fESDZDC->AliESDZDC::SetBit(AliESDZDC::kCorrectedTDCFilled, kTRUE);
1376   fESDZDC->AliESDZDC::SetBit(AliESDZDC::kEnergyCalibratedSignal, kTRUE);
1377
1378   if(esd) esd->SetZDCData(fESDZDC);
1379 }
1380
1381 //_____________________________________________________________________________
1382 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1383 {
1384   // Setting the storage
1385
1386   Bool_t deleteManager = kFALSE;
1387
1388   AliCDBManager *manager = AliCDBManager::Instance();
1389   AliCDBStorage *defstorage = manager->GetDefaultStorage();
1390
1391   if(!defstorage || !(defstorage->Contains("ZDC"))){
1392      AliWarning("No default storage set or default storage doesn't contain ZDC!");
1393      manager->SetDefaultStorage(uri);
1394      deleteManager = kTRUE;
1395   }
1396
1397   AliCDBStorage *storage = manager->GetDefaultStorage();
1398
1399   if(deleteManager){
1400     AliCDBManager::Instance()->UnsetDefaultStorage();
1401     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
1402   }
1403
1404   return storage;
1405 }
1406
1407 //_____________________________________________________________________________
1408 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1409 {
1410
1411   // Getting pedestal calibration object for ZDC set
1412
1413   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1414   if(!entry) AliFatal("No calibration data loaded!");
1415   entry->SetOwner(kFALSE);
1416
1417   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
1418   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1419
1420   return calibdata;
1421 }
1422
1423 //_____________________________________________________________________________
1424 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1425 {
1426
1427   // Getting energy and equalization calibration object for ZDC set
1428
1429   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1430   if(!entry) AliFatal("No calibration data loaded!");
1431   entry->SetOwner(kFALSE);
1432
1433   AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1434   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1435
1436   return calibdata;
1437 }
1438
1439 //_____________________________________________________________________________
1440 AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1441 {
1442
1443   // Getting energy and equalization calibration object for ZDC set
1444
1445   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1446   if(!entry) AliFatal("No calibration data loaded!");
1447   entry->SetOwner(kFALSE);
1448
1449   AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1450   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1451
1452   return calibdata;
1453 }
1454
1455 //_____________________________________________________________________________
1456 AliZDCMBCalib* AliZDCReconstructor::GetMBCalibData() const
1457 {
1458
1459   // Getting energy and equalization calibration object for ZDC set
1460
1461   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/MBCalib");
1462   if(!entry) AliFatal("No calibration data loaded!");
1463   entry->SetOwner(kFALSE);
1464
1465   AliZDCMBCalib *calibdata = dynamic_cast<AliZDCMBCalib*> (entry->GetObject());
1466   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1467
1468   return calibdata;
1469 }