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