ZDC/AliZDCReconstructor.cxx

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