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