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