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