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