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