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