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