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