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