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 sumZEM[]={0,0}, 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   for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
 839   // Ch. debug
 840   /*printf("\n ------------- CALIBRATION -------------\n");
 841   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 842         calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
 843   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 844         calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
 845   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 846         calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
 847   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 848         calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
 849   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
 850   printf(" ----------------------------------------\n");*/
 851
 852   //  ******    No. of spectator and participants nucleons
 853   //  Variables calculated to comply with ESD structure
 854   //  *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
 855   Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
 856   Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
 857   Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
 858   Double_t impPar=0., impPar1=0., impPar2=0.;
 859
 860   Bool_t energyFlag = kFALSE;
 861   // create the output tree
 862   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 863                    calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
 864                    calibZEM1, calibZEM2, sPMRef1, sPMRef2,
 865                    nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 866                    nGenSpec, nGenSpecLeft, nGenSpecRight,
 867                    nPart, nPartTotLeft, nPartTotRight,
 868                    impPar, impPar1, impPar2,
 869                    recoFlag, energyFlag, isScalerOn, scaler, tdcData);
 870
 871   const Int_t kBufferSize = 4000;
 872   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
 873   // write the output tree
 874   clustersTree->Fill();
 875   delete reco;
 876 }
 877
 878 //_____________________________________________________________________________
 879 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
 880         const Float_t* const corrADCZN1, const Float_t* const corrADCZP1,
 881         const Float_t* const corrADCZN2, const Float_t* const corrADCZP2,
 882         const Float_t* const corrADCZEM1, const Float_t* const corrADCZEM2,
 883         Float_t* sPMRef1, Float_t* sPMRef2, Bool_t isScalerOn, UInt_t* scaler,
 884         Int_t tdcData[32][4], const Int_t* const evQualityBlock,
 885         const Int_t* const triggerBlock, const Int_t* const chBlock, UInt_t puBits) const
 886 {
 887   // ****************** Reconstruct one event ******************
 888   // ---------------------- Setting reco flags for ESD
 889   UInt_t rFlags[32];
 890   for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
 891
 892   if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
 893   else rFlags[31] = 0x1;
 894   //
 895   if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
 896   if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
 897   if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
 898
 899   if(triggerBlock[0] == 1) rFlags[27] = 0x1;
 900   if(triggerBlock[1] == 1) rFlags[26] = 0x1;
 901   if(triggerBlock[2] == 1) rFlags[25] = 0x1;
 902   if(triggerBlock[3] == 1) rFlags[24] = 0x1;
 903
 904   if(chBlock[0] == 1) rFlags[18] = 0x1;
 905   if(chBlock[1] == 1) rFlags[17] = 0x1;
 906   if(chBlock[2] == 1) rFlags[16] = 0x1;
 907
 908   rFlags[13] = puBits & 0x00000020;
 909   rFlags[12] = puBits & 0x00000010;
 910   rFlags[11] = puBits & 0x00000080;
 911   rFlags[10] = puBits & 0x00000040;
 912   rFlags[9]  = puBits & 0x00000020;
 913   rFlags[8]  = puBits & 0x00000010;
 914
 915   if(corrADCZP1[0]>fSignalThreshold)  rFlags[5] = 0x1;
 916   if(corrADCZN1[0]>fSignalThreshold)  rFlags[4] = 0x1;
 917   if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
 918   if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
 919   if(corrADCZP2[0]>fSignalThreshold)  rFlags[1] = 0x1;
 920   if(corrADCZN2[0]>fSignalThreshold)  rFlags[0] = 0x1;
 921
 922   UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
 923              rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
 924              0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
 925              0x0 << 19 | rFlags[18] << 18 |  rFlags[17] << 17 |  rFlags[16] << 16 |
 926              0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
 927              rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
 928              0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
 929              rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
 930   // --------------------------------------------------
 931
 932
 933   // CH. debug
 934 /*  printf("\n*************************************************\n");
 935   printf(" ReconstructEventPbPb -> values after pedestal subtraction:\n");
 936   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 937         corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
 938   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 939         corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
 940   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 941         corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
 942   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 943         corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
 944   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
 945   printf("*************************************************\n");
 946 */
 947   // ******     Retrieving calibration data
 948   // --- Equalization coefficients ---------------------------------------------
 949   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 950   for(Int_t ji=0; ji<5; ji++){
 951      equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
 952      equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
 953      equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
 954      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
 955   }
 956   // --- Energy calibration factors ------------------------------------
 957   Float_t calibEne[6], calibSatZNA[4], calibSatZNC[4];
 958   // **** Energy calibration coefficient set to 1
 959   // **** (no trivial way to calibrate in p-p runs)
 960   for(Int_t ij=0; ij<6; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
 961   for(Int_t ij=0; ij<4; ij++){
 962     calibSatZNA[ij] = fSatCalibData->GetZNASatCalib(ij);
 963     calibSatZNC[ij] = fSatCalibData->GetZNCSatCalib(ij);
 964   }
 965
 966   // ******     Equalization of detector responses
 967   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 968   for(Int_t gi=0; gi<10; gi++){
 969      if(gi<5){
 970        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
 971        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
 972        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
 973        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
 974      }
 975      else{
 976        equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
 977        equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
 978        equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
 979        equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
 980      }
 981   }
 982
 983   // Ch. debug
 984 /*  printf("\n ------------- EQUALIZATION -------------\n");
 985   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 986         equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
 987   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 988         equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
 989   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 990         equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
 991   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
 992         equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
 993   printf(" ----------------------------------------\n");
 994 */
 995
 996   //  *** p-A RUN 2013 -> new calibration object
 997   //      to take into account saturation in ZN PMC
 998   //   -> 5th order pol. fun. to be applied BEFORE en. calibration
 999   equalTowZN1[0] = equalTowZN1[0] + calibSatZNC[0]*equalTowZN1[0]*equalTowZN1[0] +
1000         calibSatZNC[1]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0] +
1001         calibSatZNC[2]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0] +
1002         calibSatZNC[3]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0];
1003   equalTowZN2[0] = equalTowZN2[0] + calibSatZNA[0]*equalTowZN2[0]*equalTowZN2[0] +
1004         calibSatZNA[1]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0] +
1005         calibSatZNA[2]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0] +
1006         calibSatZNA[3]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0];
1007
1008   // ******     Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
1009   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
1010   for(Int_t gi=0; gi<5; gi++){
1011        calibSumZN1[0] += equalTowZN1[gi];
1012        calibSumZP1[0] += equalTowZP1[gi];
1013        calibSumZN2[0] += equalTowZN2[gi];
1014        calibSumZP2[0] += equalTowZP2[gi];
1015        //
1016        calibSumZN1[1] += equalTowZN1[gi+5];
1017        calibSumZP1[1] += equalTowZP1[gi+5];
1018        calibSumZN2[1] += equalTowZN2[gi+5];
1019        calibSumZP2[1] += equalTowZP2[gi+5];
1020   }
1021   //
1022   //fEnCalibData->Print("");
1023
1024   // High gain chain
1025   calibSumZN1[0] = calibSumZN1[0]*calibEne[0]*8.;
1026   calibSumZP1[0] = calibSumZP1[0]*calibEne[1]*8.;
1027   calibSumZN2[0] = calibSumZN2[0]*calibEne[2]*8.;
1028   calibSumZP2[0] = calibSumZP2[0]*calibEne[3]*8.;
1029   // Low gain chain
1030   calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
1031   calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
1032   calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
1033   calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
1034   //
1035   Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
1036   calibZEM1[0] = corrADCZEM1[0]*calibEne[4];
1037   calibZEM1[1] = corrADCZEM1[1]*calibEne[4];
1038   calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
1039   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
1040   for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
1041
1042   // ******     Energy calibration of detector responses
1043   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
1044   for(Int_t gi=0; gi<5; gi++){
1045      // High gain chain
1046      calibTowZN1[gi] = equalTowZN1[gi]*2*calibEne[0]*8.;
1047      calibTowZP1[gi] = equalTowZP1[gi]*2*calibEne[1]*8.;
1048      calibTowZN2[gi] = equalTowZN2[gi]*2*calibEne[2]*8.;
1049      calibTowZP2[gi] = equalTowZP2[gi]*2*calibEne[3]*8.;
1050      // Low gain chain
1051      calibTowZN1[gi+5] = equalTowZN1[gi+5]*2*calibEne[0];
1052      calibTowZP1[gi+5] = equalTowZP1[gi+5]*2*calibEne[1];
1053      calibTowZN2[gi+5] = equalTowZN2[gi+5]*2*calibEne[2];
1054      calibTowZP2[gi+5] = equalTowZP2[gi+5]*2*calibEne[3];
1055   }
1056
1057   // Ch. debug
1058 /*  printf("\n ------------- CALIBRATION -------------\n");
1059   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1060         calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
1061   printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1062         calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
1063   printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1064         calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
1065   printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1066         calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
1067   printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
1068   printf(" ----------------------------------------\n");
1069 */
1070   //  ******    Number of detected spectator nucleons
1071   Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
1072   if(fBeamEnergy>0.01){
1073     nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
1074     nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
1075     nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
1076     nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
1077   }
1078   else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
1079   /*printf("\n\t AliZDCReconstructor -> fBeamEnergy %1.0f: nDetSpecNsideA %d, nDetSpecPsideA %d,"
1080     " nDetSpecNsideC %d, nDetSpecPsideC %d\n",fBeamEnergy,nDetSpecNLeft, nDetSpecPLeft,
1081     nDetSpecNRight, nDetSpecPRight);*/
1082
1083   Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
1084   Int_t nPart=0, nPartA=0, nPartC=0;
1085   Double_t b=0., bA=0., bC=0.;
1086
1087   if(fIsCalibrationMB == kFALSE){
1088    // ******   Reconstruction parameters ------------------
1089    if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
1090    if(!fgRecoParam){
1091      AliError("  RecoParam object not retrieved correctly: not reconstructing ZDC event!!!");
1092      return;
1093    }
1094    TH1D* hNpartDist = fgRecoParam->GethNpartDist();
1095    TH1D*   hbDist = fgRecoParam->GethbDist();
1096    Float_t  fClkCenter = fgRecoParam->GetClkCenter();
1097    if(!hNpartDist || !hbDist){
1098       AliError("Something wrong in Glauber MC histos got from AliZDCREcoParamPbPb: NO EVENT RECO FOR ZDC DATA!!!\n\n");
1099       //return;
1100    }
1101    else{
1102     if(!fgMBCalibData) fgMBCalibData = const_cast<AliZDCMBCalib*>(GetMBCalibData());
1103     TH2F *hZDCvsZEM  = fgMBCalibData->GethZDCvsZEM();
1104     TH2F *hZDCCvsZEM = fgMBCalibData->GethZDCCvsZEM();
1105     TH2F *hZDCAvsZEM = fgMBCalibData->GethZDCAvsZEM();
1106     //
1107     Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
1108     Double_t origin = xHighEdge*fClkCenter;
1109     // Ch. debug
1110     //printf("\n\n  xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
1111     //
1112     // ====> Summed ZDC info (sideA+side C)
1113     TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
1114     Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
1115     Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
1116     line->SetParameter(0, y/(x-origin));
1117     line->SetParameter(1, -origin*y/(x-origin));
1118     // Ch. debug
1119     //printf("  ***************** Summed ZDC info (sideA+side C) \n");
1120     //printf("  E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f   ", x, y,y/(x-origin),-origin*y/(x-origin));
1121     //
1122     Double_t countPerc=0;
1123     Double_t xBinCenter=0, yBinCenter=0;
1124     for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
1125       for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
1126          xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
1127          yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
1128          //
1129          if(line->GetParameter(0)>0){
1130            if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
1131              countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
1132              // Ch. debug
1133              //printf(" xBinCenter  %1.3f, yBinCenter %1.0f,  countPerc %1.0f\n",
1134                 //xBinCenter, yBinCenter, countPerc);
1135            }
1136          }
1137          else{
1138            if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
1139              countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
1140              // Ch. debug
1141              //printf(" xBinCenter  %1.3f, yBinCenter %1.0f,  countPerc %1.0f\n",
1142                 //xBinCenter, yBinCenter, countPerc);
1143            }
1144          }
1145       }
1146     }
1147     //
1148     Double_t xSecPerc = 0.;
1149     if(hZDCvsZEM->GetEntries()!=0){
1150       xSecPerc = countPerc/hZDCvsZEM->GetEntries();
1151     }
1152     else{
1153       AliWarning("  Histogram hZDCvsZEM from OCDB has no entries!!!");
1154     }
1155     // Ch. debug
1156     //printf("  xSecPerc %1.4f  \n", xSecPerc);
1157
1158     // ====> side C
1159     TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
1160     Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
1161     lineC->SetParameter(0, yC/(x-origin));
1162     lineC->SetParameter(1, -origin*yC/(x-origin));
1163     // Ch. debug
1164     //printf("  ***************** Side C \n");
1165     //printf("  E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f   ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
1166     //
1167     Double_t countPercC=0;
1168     Double_t xBinCenterC=0, yBinCenterC=0;
1169     for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
1170       for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
1171          xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
1172          yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
1173          if(lineC->GetParameter(0)>0){
1174            if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1175              countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1176            }
1177          }
1178          else{
1179            if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1180              countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1181            }
1182          }
1183       }
1184     }
1185     //
1186     Double_t xSecPercC = 0.;
1187     if(hZDCCvsZEM->GetEntries()!=0){
1188       xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
1189     }
1190     else{
1191       AliWarning("  Histogram hZDCCvsZEM from OCDB has no entries!!!");
1192     }
1193     // Ch. debug
1194     //printf("  xSecPercC %1.4f  \n", xSecPercC);
1195
1196     // ====> side A
1197     TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
1198     Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
1199     lineA->SetParameter(0, yA/(x-origin));
1200     lineA->SetParameter(1, -origin*yA/(x-origin));
1201     //
1202     // Ch. debug
1203     //printf("  ***************** Side A \n");
1204     //printf("  E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f   ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
1205     //
1206     Double_t countPercA=0;
1207     Double_t xBinCenterA=0, yBinCenterA=0;
1208     for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
1209       for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
1210          xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
1211          yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
1212          if(lineA->GetParameter(0)>0){
1213            if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1214              countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1215            }
1216          }
1217          else{
1218            if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1219              countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1220            }
1221          }
1222       }
1223     }
1224     //
1225     Double_t xSecPercA = 0.;
1226     if(hZDCAvsZEM->GetEntries()!=0){
1227       xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
1228     }
1229     else{
1230       AliWarning("  Histogram hZDCAvsZEM from OCDB has no entries!!!");
1231     }
1232     // Ch. debug
1233     //printf("  xSecPercA %1.4f  \n", xSecPercA);
1234
1235     //  ******    Number of participants (from E_ZDC vs. E_ZEM correlation)
1236     Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1237     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1238       nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1239       if((1.-nPartFrac) < xSecPerc){
1240         nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1241         // Ch. debug
1242         //printf("  ***************** Summed ZDC info (sideA+side C) \n");
1243         //printf("  nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1244         break;
1245       }
1246     }
1247     if(nPart<0) nPart=0;
1248     //
1249     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1250       nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1251       if((1.-nPartFracC) < xSecPercC){
1252         nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1253         // Ch. debug
1254         //printf("  ***************** Side C \n");
1255         //printf("  nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1256         break;
1257     }
1258     }
1259     if(nPartC<0) nPartC=0;
1260     //
1261     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1262       nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1263       if((1.-nPartFracA) < xSecPercA){
1264         nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1265         // Ch. debug
1266         //printf("  ***************** Side A \n");
1267         //printf("  nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1268         break;
1269       }
1270     }
1271     if(nPartA<0) nPartA=0;
1272
1273     //  ******    Impact parameter (from E_ZDC vs. E_ZEM correlation)
1274     Double_t bFrac=0., bFracC=0., bFracA=0.;
1275     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1276       bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1277       if(bFrac > xSecPerc){
1278         b = hbDist->GetBinLowEdge(ibbin);
1279         break;
1280       }
1281     }
1282     //
1283     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1284       bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1285       if(bFracC > xSecPercC){
1286         bC = hbDist->GetBinLowEdge(ibbin);
1287         break;
1288       }
1289     }
1290     //
1291     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1292       bFracA += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1293       if(bFracA > xSecPercA){
1294         bA = hbDist->GetBinLowEdge(ibbin);
1295         break;
1296       }
1297     }
1298
1299     //  ******  Number of spectator nucleons
1300     nGenSpec = 416 - nPart;
1301     nGenSpecC = 416 - nPartC;
1302     nGenSpecA = 416 - nPartA;
1303     if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1304     if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1305     if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1306
1307     delete line;
1308     delete lineC;  delete lineA;
1309    }
1310   } // ONLY IF fIsCalibrationMB==kFALSE
1311
1312   Bool_t energyFlag = kTRUE;
1313   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1314                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1315                   calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1316                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1317                   nGenSpec, nGenSpecA, nGenSpecC,
1318                   nPart, nPartA, nPartC, b, bA, bC,
1319                   recoFlag, energyFlag, isScalerOn, scaler, tdcData);
1320
1321   const Int_t kBufferSize = 4000;
1322   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1323   //reco->Print("");
1324   // write the output tree
1325   clustersTree->Fill();
1326   delete reco;
1327 }
1328
1329
1330 //_____________________________________________________________________________
1331 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1332 {
1333   // fill energies and number of participants to the ESD
1334
1335   // Retrieving TDC calibration data
1336   // Parameters for TDC centering around zero
1337   int const knTDC = 6;
1338   Float_t tdcOffset[knTDC];
1339   for(Int_t jj=0; jj<knTDC; jj++) tdcOffset[jj] = fTDCCalibData->GetMeanTDC(jj);
1340   //fTDCCalibData->Print("");
1341
1342   AliZDCReco reco;
1343   AliZDCReco* preco = &reco;
1344   clustersTree->SetBranchAddress("ZDC", &preco);
1345   clustersTree->GetEntry(0);
1346   //
1347   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1348   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1349   for(Int_t i=0; i<5; i++){
1350      tZN1Ene[i] = reco.GetZN1HREnTow(i);
1351      tZN2Ene[i] = reco.GetZN2HREnTow(i);
1352      tZP1Ene[i] = reco.GetZP1HREnTow(i);
1353      tZP2Ene[i] = reco.GetZP2HREnTow(i);
1354      //
1355      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1356      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1357      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1358      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1359   }
1360   //
1361   fESDZDC->SetZN1TowerEnergy(tZN1Ene);
1362   fESDZDC->SetZN2TowerEnergy(tZN2Ene);
1363   fESDZDC->SetZP1TowerEnergy(tZP1Ene);
1364   fESDZDC->SetZP2TowerEnergy(tZP2Ene);
1365   //
1366   fESDZDC->SetZN1TowerEnergyLR(tZN1EneLR);
1367   fESDZDC->SetZN2TowerEnergyLR(tZN2EneLR);
1368   fESDZDC->SetZP1TowerEnergyLR(tZP1EneLR);
1369   fESDZDC->SetZP2TowerEnergyLR(tZP2EneLR);
1370   //
1371   Int_t nPart  = reco.GetNParticipants();
1372   Int_t nPartA = reco.GetNPartSideA();
1373   Int_t nPartC = reco.GetNPartSideC();
1374   Double_t b  = reco.GetImpParameter();
1375   Double_t bA = reco.GetImpParSideA();
1376   Double_t bC = reco.GetImpParSideC();
1377   UInt_t recoFlag = reco.GetRecoFlag();
1378
1379   fESDZDC->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(),
1380         reco.GetZEM1HRsignal(), reco.GetZEM2HRsignal(),
1381         reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1382         nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1383
1384   // Writing ZDC scaler for cross section calculation
1385   // ONLY IF the scaler has been read during the event
1386   if(reco.IsScalerOn()==kTRUE){
1387     UInt_t counts[32];
1388     for(Int_t jk=0; jk<32; jk++) counts[jk] = reco.GetZDCScaler(jk);
1389     fESDZDC->SetZDCScaler(counts);
1390   }
1391
1392   Int_t tdcValues[32][4] = {{0,}};
1393   Float_t tdcCorrected[32][4] = {{9999.,}};
1394   for(Int_t jk=0; jk<32; jk++){
1395     for(Int_t lk=0; lk<4; lk++){
1396       tdcValues[jk][lk] = reco.GetZDCTDCData(jk, lk);
1397       //
1398       if(jk==8 && TMath::Abs(tdcValues[jk][lk])>1e-09)      fESDZDC->SetZEM1TDChit(kTRUE);
1399       else if(jk==9 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZEM2TDChit(kTRUE);
1400       else if(jk==10 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZNCTDChit(kTRUE);
1401       else if(jk==11 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZPCTDChit(kTRUE);
1402       else if(jk==12 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZNATDChit(kTRUE);
1403       else if(jk==13 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZPATDChit(kTRUE);
1404       //Ch debug
1405       //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]);
1406     }
1407   }
1408
1409   // Writing TDC data into ZDC ESDs
1410   // 4/2/2011 -> Subtracting L0 (tdcValues[15]) instead of ADC gate
1411   // we try to keep the TDC oscillations as low as possible!
1412   for(Int_t jk=0; jk<32; jk++){
1413     for(Int_t lk=0; lk<4; lk++){
1414       if(tdcValues[jk][lk]!=0.){
1415         // Feb2013 _-> TDC correct entry is there ONLY IF tdc has a hit!
1416         if(TMath::Abs(tdcValues[jk][lk])>1e-09){
1417            tdcCorrected[jk][lk] = 0.025*(tdcValues[jk][lk]-tdcValues[15][0])+fMeanPhase;
1418            // Sep 2011: TDC ch. from 8 to 13 centered around 0 using OCDB
1419            if(jk>=8 && jk<=13) tdcCorrected[jk][lk] =  tdcCorrected[jk][lk] - tdcOffset[jk-8];
1420            //Ch. debug
1421            //if(jk>=8 && jk<=13) printf(" *** tdcOffset%d %f  tdcCorr%d %f \n",jk,tdcOffset[jk-8],tdcCorrected[jk][lk]);
1422         }
1423       }
1424     }
1425   }
1426
1427   fESDZDC->SetZDCTDCData(tdcValues);
1428   fESDZDC->SetZDCTDCCorrected(tdcCorrected);
1429   fESDZDC->AliESDZDC::SetBit(AliESDZDC::kCorrectedTDCFilled, reco.GetEnergyFlag());
1430   fESDZDC->AliESDZDC::SetBit(AliESDZDC::kEnergyCalibratedSignal, kTRUE);
1431
1432   if(esd) esd->SetZDCData(fESDZDC);
1433 }
1434
1435 //_____________________________________________________________________________
1436 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1437 {
1438   // Setting the storage
1439
1440   Bool_t deleteManager = kFALSE;
1441
1442   AliCDBManager *manager = AliCDBManager::Instance();
1443   AliCDBStorage *defstorage = manager->GetDefaultStorage();
1444
1445   if(!defstorage || !(defstorage->Contains("ZDC"))){
1446      AliWarning("No default storage set or default storage doesn't contain ZDC!");
1447      manager->SetDefaultStorage(uri);
1448      deleteManager = kTRUE;
1449   }
1450
1451   AliCDBStorage *storage = manager->GetDefaultStorage();
1452
1453   if(deleteManager){
1454     AliCDBManager::Instance()->UnsetDefaultStorage();
1455     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
1456   }
1457
1458   return storage;
1459 }
1460
1461 //_____________________________________________________________________________
1462 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1463 {
1464
1465   // Getting pedestal calibration object for ZDC set
1466
1467   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1468   if(!entry) AliFatal("No calibration data loaded!");
1469   entry->SetOwner(kFALSE);
1470
1471   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
1472   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1473
1474   return calibdata;
1475 }
1476
1477 //_____________________________________________________________________________
1478 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1479 {
1480
1481   // Getting energy and equalization calibration object for ZDC set
1482
1483   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1484   if(!entry) AliFatal("No calibration data loaded!");
1485   entry->SetOwner(kFALSE);
1486
1487   AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1488   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1489
1490   return calibdata;
1491 }
1492
1493 //_____________________________________________________________________________
1494 AliZDCSaturationCalib* AliZDCReconstructor::GetSaturationCalibData() const
1495 {
1496
1497   // Getting energy and equalization calibration object for ZDC set
1498
1499   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/SaturationCalib");
1500   if(!entry) AliFatal("No calibration data loaded!");
1501   entry->SetOwner(kFALSE);
1502
1503   AliZDCSaturationCalib *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
1515   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1516   if(!entry) AliFatal("No calibration data loaded!");
1517   entry->SetOwner(kFALSE);
1518
1519   AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1520   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1521
1522   return calibdata;
1523 }
1524
1525 //_____________________________________________________________________________
1526 AliZDCMBCalib* AliZDCReconstructor::GetMBCalibData() const
1527 {
1528
1529   // Getting energy and equalization calibration object for ZDC set
1530
1531   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/MBCalib");
1532   if(!entry) AliFatal("No calibration data loaded!");
1533   entry->SetOwner(kFALSE);
1534
1535   AliZDCMBCalib *calibdata = dynamic_cast<AliZDCMBCalib*> (entry->GetObject());
1536   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1537
1538   return calibdata;
1539 }
1540
1541 //_____________________________________________________________________________
1542 AliZDCTDCCalib* AliZDCReconstructor::GetTDCCalibData() const
1543 {
1544
1545   // Getting TDC object for ZDC
1546
1547   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TDCCalib");
1548   if(!entry) AliFatal("No calibration data loaded!");
1549   entry->SetOwner(kFALSE);
1550
1551   AliZDCTDCCalib *calibdata = dynamic_cast<AliZDCTDCCalib*> (entry->GetObject());
1552   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1553
1554   return calibdata;
1555 }