ZDC/AliZDCReconstructor.cxx

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