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 "AliGRPObject.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDZDC.h"
  39 #include "AliZDCDigit.h"
  40 #include "AliZDCRawStream.h"
  41 #include "AliZDCReco.h"
  42 #include "AliZDCReconstructor.h"
  43 #include "AliZDCPedestals.h"
  44 #include "AliZDCEnCalib.h"
  45 #include "AliZDCTowerCalib.h"
  46 #include "AliZDCRecoParam.h"
  47 #include "AliZDCRecoParampp.h"
  48 #include "AliZDCRecoParamPbPb.h"
  49
  50
  51 ClassImp(AliZDCReconstructor)
  52 AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0;  //reconstruction parameters
  53
  54 //_____________________________________________________________________________
  55 AliZDCReconstructor:: AliZDCReconstructor() :
  56   fPedData(GetPedData()),
  57   fEnCalibData(GetEnCalibData()),
  58   fTowCalibData(GetTowCalibData()),
  59   fRecoMode(0),
  60   fBeamEnergy(0.),
  61   fNRun(0),
  62   fIsCalibrationMB(kFALSE),
  63   fPedSubMode(0),
  64   fRecoFlag(0x0)
  65 {
  66   // **** Default constructor
  67
  68 }
  69
  70
  71 //_____________________________________________________________________________
  72 AliZDCReconstructor::~AliZDCReconstructor()
  73 {
  74 // destructor
  75    if(fRecoParam)    delete fRecoParam;
  76    if(fPedData)      delete fPedData;
  77    if(fEnCalibData)  delete fEnCalibData;
  78    if(fTowCalibData) delete fTowCalibData;
  79 }
  80
  81 //____________________________________________________________________________
  82 void AliZDCReconstructor::Init()
  83 {
  84   // Setting reconstruction mode
  85   // Getting beam type and beam energy from GRP calibration object
  86
  87   if(fRecoMode==0 && fBeamEnergy==0.){
  88     // Initialization of the GRP entry
  89     AliCDBEntry*  entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
  90     AliGRPObject* grpData = 0x0;
  91     if(entry){
  92       TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
  93       if(m){
  94         m->Print();
  95         grpData = new AliGRPObject();
  96         grpData->ReadValuesFromMap(m);
  97       }
  98       else{
  99         grpData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
 100       }
 101       entry->SetOwner(0);
 102       AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
 103     }
 104     if(!grpData) AliError("No GRP entry found in OCDB!");
 105
 106     TString runType = grpData->GetRunType();
 107     if(runType==AliGRPObject::GetInvalidString()){
 108       AliError("GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN");
 109       runType = "UNKNOWN";
 110     }
 111     if((runType.CompareTo("CALIBRATION_MB")) == 0){
 112       fIsCalibrationMB = kTRUE;
 113       //
 114       fRecoParam = new AliZDCRecoParamPbPb();
 115       //
 116       TH2F* hZDCvsZEM = new TH2F("hZDCvsZEM","hZDCvsZEM",100,0.,10.,100,0.,1000.);
 117       hZDCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCvsZEM->SetYTitle("E_{ZDC} (TeV)");
 118       fRecoParam->SetZDCvsZEM(hZDCvsZEM);
 119       //
 120       TH2F* hZDCCvsZEM = new TH2F("hZDCCvsZEM","hZDCCvsZEM",100,0.,10.,100,0.,500.);
 121       hZDCCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCCvsZEM->SetYTitle("E_{ZDCC} (TeV)");
 122       fRecoParam->SetZDCCvsZEM(hZDCCvsZEM);
 123       //
 124       TH2F* hZDCAvsZEM = new TH2F("hZDCAvsZEM","hZDCAvsZEM",100,0.,10.,100,0.,500.);
 125       hZDCAvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCAvsZEM->SetYTitle("E_{ZDCA} (TeV)");
 126       fRecoParam->SetZDCAvsZEM(hZDCAvsZEM);
 127     }
 128
 129     TString beamType = grpData->GetBeamType();
 130     if(beamType==AliGRPObject::GetInvalidString()){
 131       AliError("GRP/GRP/Data entry:  missing value for the beam energy !");
 132       AliError("\t ZDC does not reconstruct event 4 UNKNOWN beam type\n");
 133       return;
 134     }
 135     if((beamType.CompareTo("p-p")) == 0){
 136       fRecoMode=1;
 137       fRecoParam = (AliZDCRecoParampp*) GetppRecoParamFromOCDB();
 138     }
 139     else if((beamType.CompareTo("A-A")) == 0){
 140       fRecoMode=2;
 141       if(fIsCalibrationMB == kFALSE)
 142          fRecoParam = (AliZDCRecoParamPbPb*) GetPbPbRecoParamFromOCDB();
 143     }
 144
 145     fBeamEnergy = grpData->GetBeamEnergy();
 146     if(fBeamEnergy==AliGRPObject::GetInvalidFloat()){
 147       AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0.");
 148       fBeamEnergy = 0.;
 149     }
 150
 151     if(fIsCalibrationMB==kTRUE){
 152       AliInfo("\n ***** CALIBRATION_MB data -> building AliZDCRecoParamPbPb object *****");
 153     }
 154     else{
 155       AliInfo(Form("\n\n ***** ZDC reconstruction initialized for %s @ %1.3f GeV *****\n",beamType.Data(), fBeamEnergy));
 156     }
 157   }
 158   else{
 159     AliError(" ATTENTION!!!!!! No beam type nor beam energy has been set!!!!!!\n");
 160   }
 161
 162 }
 163
 164 //_____________________________________________________________________________
 165 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
 166 {
 167   // *** Local ZDC reconstruction for digits
 168   // Works on the current event
 169
 170   // Retrieving calibration data
 171   // Parameters for mean value pedestal subtraction
 172   int const kNch = 24;
 173   Float_t meanPed[2*kNch];
 174   for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
 175   // Parameters pedestal subtraction through correlation with out-of-time signals
 176   Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
 177   for(Int_t jj=0; jj<2*kNch; jj++){
 178      corrCoeff0[jj] =  fPedData->GetPedCorrCoeff0(jj);
 179      corrCoeff1[jj] =  fPedData->GetPedCorrCoeff1(jj);
 180   }
 181
 182   // get digits
 183   AliZDCDigit digit;
 184   AliZDCDigit* pdigit = &digit;
 185   digitsTree->SetBranchAddress("ZDC", &pdigit);
 186   //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
 187
 188   // loop over digits
 189   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
 190   Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
 191   for(Int_t i=0; i<10; i++){
 192      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
 193      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
 194   }
 195
 196   Int_t digNentries = digitsTree->GetEntries();
 197   Float_t ootDigi[kNch];
 198   // -- Reading out-of-time signals (last kNch entries) for current event
 199   if(fPedSubMode==1){
 200     for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
 201        ootDigi[iDigit] = digitsTree->GetEntry(iDigit);
 202     }
 203   }
 204
 205   for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
 206    digitsTree->GetEntry(iDigit);
 207    if (!pdigit) continue;
 208    //
 209    Int_t det = digit.GetSector(0);
 210    Int_t quad = digit.GetSector(1);
 211    Int_t pedindex = -1;
 212    Float_t ped2SubHg=0., ped2SubLg=0.;
 213    if(quad!=5){
 214      if(det==1)      pedindex = quad;
 215      else if(det==2) pedindex = quad+5;
 216      else if(det==3) pedindex = quad+9;
 217      else if(det==4) pedindex = quad+12;
 218      else if(det==5) pedindex = quad+17;
 219    }
 220    else pedindex = (det-1)/3+22;
 221    //
 222    if(fPedSubMode==0){
 223      ped2SubHg = meanPed[pedindex];
 224      ped2SubLg = meanPed[pedindex+kNch];
 225    }
 226    else if(fPedSubMode==1){
 227      ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
 228      ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
 229    }
 230
 231    if(quad != 5){ // ZDC (not reference PTMs!)
 232     if(det == 1){ // *** ZNC
 233        tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
 234        tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
 235        if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
 236        if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
 237        // Ch. debug
 238        //printf("\t pedindex %d tZN1Corr[%d] = %1.0f tZN1Corr[%d] = %1.0f",
 239        //       pedindex, quad, tZN1Corr[quad], quad+5, tZN1Corr[quad+5]);
 240     }
 241     else if(det == 2){ // *** ZP1
 242        tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
 243        tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
 244        if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
 245        if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
 246        // Ch. debug
 247        //printf("\t pedindex %d tZP1Corr[%d] = %1.0f tZP1Corr[%d] = %1.0f",
 248        //       pedindex, quad, tZP1Corr[quad], quad+5, tZP1Corr[quad+5]);
 249     }
 250     else if(det == 3){
 251        if(quad == 1){       // *** ZEM1
 252          dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
 253          dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
 254          if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
 255          if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
 256          // Ch. debug
 257          //printf("\t pedindex %d tZEM1Corr[%d] = %1.0f tZEM1Corr[%d] = %1.0f",
 258          //     pedindex, quad, tZEM1Corr[quad], quad+1, tZEM1Corr[quad+1]);
 259        }
 260        else if(quad == 2){  // *** ZEM2
 261          dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
 262          dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
 263          if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
 264          if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
 265          // Ch. debug
 266          //printf("\t pedindex %d tZEM2Corr[%d] = %1.0f tZEM2Corr[%d] = %1.0f",
 267          //     pedindex, quad, tZEM2Corr[quad], quad+1, tZEM2Corr[quad+1]);
 268        }
 269     }
 270     else if(det == 4){  // *** ZN2
 271        tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
 272        tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
 273        if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
 274        if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
 275        // Ch. debug
 276        //printf("\t pedindex %d tZN2Corr[%d] = %1.0f tZN2Corr[%d] = %1.0f\n",
 277        //       pedindex, quad, tZN2Corr[quad], quad+5, tZN2Corr[quad+5]);
 278     }
 279     else if(det == 5){  // *** ZP2
 280        tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
 281        tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
 282        if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
 283        if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
 284        // Ch. debug
 285        //printf("\t pedindex %d tZP2Corr[%d] = %1.0f tZP2Corr[%d] = %1.0f\n",
 286        //       pedindex, quad, tZP2Corr[quad], quad+5, tZP2Corr[quad+5]);
 287     }
 288    }
 289    else{ // Reference PMs
 290      if(det == 1){
 291        sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
 292        sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
 293        // Ch. debug
 294        if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
 295        if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
 296      }
 297      else if(det == 4){
 298        sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
 299        sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
 300        // Ch. debug
 301        if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
 302        if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
 303      }
 304    }
 305
 306    // Ch. debug
 307    /*printf(" - AliZDCReconstructor -> digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
 308          iDigit, det, quad, ped2SubHg, ped2SubLg);
 309    printf("   HGChain -> RawDig %d DigCorr %1.2f\n", digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg);
 310    printf("   LGChain -> RawDig %d DigCorr %1.2f\n", digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);
 311    */
 312   }//digits loop
 313
 314   // If CALIBRATION_MB run build the RecoParam object
 315   if(fIsCalibrationMB){
 316     Float_t ZDCC=0., ZDCA=0., ZEM=0;
 317     ZEM += dZEM1Corr[0] + dZEM2Corr[0];
 318     for(Int_t jkl=0; jkl<5; jkl++){
 319        ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
 320        ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
 321     }
 322     // Using energies in TeV in fRecoParam object
 323     BuildRecoParam(fRecoParam->GethZDCvsZEM(), fRecoParam->GethZDCCvsZEM(),
 324                    fRecoParam->GethZDCAvsZEM(), ZDCC/1000., ZDCA/1000., ZEM/1000.);
 325   }
 326   // reconstruct the event
 327   if(fRecoMode==1)
 328     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 329       dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2);
 330   else if(fRecoMode==2)
 331     ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 332       dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2);
 333 }
 334
 335 //_____________________________________________________________________________
 336 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree)
 337 {
 338   // *** ZDC raw data reconstruction
 339   // Works on the current event
 340
 341   // Retrieving calibration data
 342   // Parameters for pedestal subtraction
 343   int const kNch = 24;
 344   Float_t meanPed[2*kNch];
 345   for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
 346   // Parameters pedestal subtraction through correlation with out-of-time signals
 347   Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
 348   for(Int_t jj=0; jj<2*kNch; jj++){
 349      corrCoeff0[jj] =  fPedData->GetPedCorrCoeff0(jj);
 350      corrCoeff1[jj] =  fPedData->GetPedCorrCoeff1(jj);
 351   }
 352
 353   Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
 354   Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
 355   Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
 356   Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
 357   Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
 358   Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
 359   for(Int_t ich=0; ich<5; ich++){
 360     adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
 361     adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
 362     adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
 363     adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
 364     if(ich<2){
 365       adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
 366       pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
 367     }
 368   }
 369
 370   // loop over raw data
 371   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
 372   Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
 373   for(Int_t i=0; i<10; i++){
 374      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
 375      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
 376   }
 377   //
 378   rawReader->Reset();
 379   fNRun = (Int_t) rawReader->GetRunNumber();
 380   AliZDCRawStream rawData(rawReader);
 381   while(rawData.Next()){
 382    // Do
 383    Bool_t ch2process = kTRUE;
 384    //
 385    if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)){
 386      fRecoFlag = 0x1<< 9;
 387      ch2process = kFALSE;
 388    }
 389    if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)){
 390      fRecoFlag = 0x1 << 8;
 391      ch2process = kFALSE;
 392    }
 393    if(rawData.GetNChannelsOn() < 48 ) fRecoFlag = 0x1 << 2;
 394
 395    if((rawData.IsADCDataWord()) && (ch2process == kTRUE)){
 396
 397      Int_t adcMod = rawData.GetADCModule();
 398      Int_t det = rawData.GetSector(0);
 399      Int_t quad = rawData.GetSector(1);
 400      Int_t gain = rawData.GetADCGain();
 401      Int_t pedindex=0;
 402      //
 403      // Mean pedestal value subtraction -------------------------------------------------------
 404      if(fPedSubMode == 0){
 405        // Not interested in o.o.t. signals (ADC modules 2, 3)
 406        if(adcMod == 2 || adcMod == 3) return;
 407        //
 408        if(quad != 5){ // ZDCs (not reference PTMs)
 409         if(det == 1){
 410           pedindex = quad;
 411           if(gain == 0) tZN1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 412           else tZN1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 413         }
 414         else if(det == 2){
 415           pedindex = quad+5;
 416           if(gain == 0) tZP1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 417           else tZP1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 418         }
 419         else if(det == 3){
 420           pedindex = quad+9;
 421           if(quad==1){
 422             if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 423             else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 424           }
 425           else if(quad==2){
 426             if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 427             else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 428           }
 429         }
 430         else if(det == 4){
 431           pedindex = quad+12;
 432           if(gain == 0) tZN2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 433           else tZN2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 434         }
 435         else if(det == 5){
 436           pedindex = quad+17;
 437           if(gain == 0) tZP2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 438           else tZP2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 439         }
 440        }
 441        else{ // reference PM
 442          pedindex = (det-1)/3 + 22;
 443          if(det == 1){
 444            if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 445          else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 446          }
 447          else if(det == 4){
 448            if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 449            else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 450          }
 451        }
 452        // Ch. debug
 453        /*printf(" -> AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f\n",
 454          det,quad,gain, pedindex, meanPed[pedindex]);
 455        printf(" -> AliZDCReconstructor: RawADC %1.0f ADCCorr %1.0f\n",
 456          rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);*/
 457
 458      }// mean pedestal subtraction
 459      // Pedestal subtraction from correlation ------------------------------------------------
 460      else if(fPedSubMode == 1){
 461        // In time signals
 462        if(adcMod==0 || adcMod==1){
 463          if(quad != 5){ // signals from ZDCs
 464            if(det == 1){
 465              if(gain==0) adcZN1[quad] = rawData.GetADCValue();
 466              else adcZN1lg[quad] = rawData.GetADCValue();
 467            }
 468            else if(det == 2){
 469              if(gain==0) adcZP1[quad] = rawData.GetADCValue();
 470              else adcZP1lg[quad] = rawData.GetADCValue();
 471            }
 472            else if(det == 3){
 473              if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
 474              else adcZEMlg[quad-1] = rawData.GetADCValue();
 475            }
 476            else if(det == 4){
 477              if(gain==0) adcZN2[quad] = rawData.GetADCValue();
 478              else adcZN2lg[quad] = rawData.GetADCValue();
 479            }
 480            else if(det == 5){
 481              if(gain==0) adcZP2[quad] = rawData.GetADCValue();
 482              else adcZP2lg[quad] = rawData.GetADCValue();
 483            }
 484          }
 485          else{ // signals from reference PM
 486             if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
 487             else pmReflg[quad-1] = rawData.GetADCValue();
 488          }
 489        }
 490        // Out-of-time pedestals
 491        else if(adcMod==2 || adcMod==3){
 492          if(quad != 5){ // signals from ZDCs
 493            if(det == 1){
 494              if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
 495              else adcZN1ootlg[quad] = rawData.GetADCValue();
 496            }
 497            else if(det == 2){
 498              if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
 499              else adcZP1ootlg[quad] = rawData.GetADCValue();
 500            }
 501            else if(det == 3){
 502              if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
 503              else adcZEMootlg[quad-1] = rawData.GetADCValue();
 504            }
 505            else if(det == 4){
 506              if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
 507              else adcZN2ootlg[quad] = rawData.GetADCValue();
 508            }
 509            else if(det == 5){
 510              if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
 511              else adcZP2ootlg[quad] = rawData.GetADCValue();
 512            }
 513          }
 514          else{ // signals from reference PM
 515             if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
 516             else pmRefootlg[quad-1] = rawData.GetADCValue();
 517          }
 518        }
 519      } // pedestal subtraction from correlation
 520         // Ch. debug
 521         //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
 522         //  det,quad,gain, pedindex, meanPed[pedindex]);
 523    }//IsADCDataWord
 524   }//loop on raw data
 525
 526   if(fPedSubMode==1){
 527     for(Int_t t=0; t<5; t++){
 528        tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
 529        tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
 530        //
 531        tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
 532        tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
 533        //
 534        tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
 535        tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
 536        //
 537        tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
 538        tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
 539        // 0---------0 Ch. debug 0---------0
 540 /*       printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
 541        printf("\tCorrCoeff0\tCorrCoeff1\n");
 542        printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
 543        printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
 544        printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
 545        printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
 546        printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
 547        printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
 548        printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
 549        printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
 550
 551        printf("ZN1 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
 552                 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
 553        printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
 554                 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
 555        //
 556        printf("ZP1 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
 557                 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
 558        printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
 559                 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
 560        //
 561        printf("ZN2 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
 562                 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
 563        printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
 564                 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
 565        //
 566        printf("ZP2 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
 567                 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
 568        printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
 569                 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
 570 */
 571     }
 572     dZEM1Corr[0] = adcZEM[0]   - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
 573     dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
 574     dZEM2Corr[0] = adcZEM[1]   - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
 575     dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
 576     //
 577     sPMRef1[0] = pmRef[0]   - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
 578     sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
 579     sPMRef2[0] = pmRef[0]   - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
 580     sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
 581   }
 582
 583   // If CALIBRATION_MB run build the RecoParam object
 584   if(fIsCalibrationMB){
 585     Float_t ZDCC=0., ZDCA=0., ZEM=0;
 586     ZEM += dZEM1Corr[0] + dZEM2Corr[0];
 587     for(Int_t jkl=0; jkl<5; jkl++){
 588        ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
 589        ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
 590     }
 591     BuildRecoParam(fRecoParam->GethZDCvsZEM(), fRecoParam->GethZDCCvsZEM(),
 592                    fRecoParam->GethZDCAvsZEM(), ZDCC/100., ZDCA/100., ZEM/100.);
 593   }
 594   // reconstruct the event
 595   else{
 596     if(fRecoMode==1) // p-p data
 597       ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 598         dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2);
 599     else if(fRecoMode==2) // Pb-Pb data
 600       ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 601         dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2);
 602   }
 603 }
 604
 605 //_____________________________________________________________________________
 606 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* corrADCZN1,
 607         Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
 608         Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2) const
 609 {
 610   // ****************** Reconstruct one event ******************
 611
 612   // ******     Retrieving calibration data
 613   // --- Equalization coefficients ---------------------------------------------
 614   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 615   for(Int_t ji=0; ji<5; ji++){
 616      equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
 617      equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
 618      equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
 619      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
 620   }
 621   // --- Energy calibration factors ------------------------------------
 622   Float_t calibEne[4];
 623   // **** Energy calibration coefficient set to 1
 624   // **** (no trivial way to calibrate in p-p runs)
 625   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
 626
 627   // ******     Equalization of detector responses
 628   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 629   for(Int_t gi=0; gi<10; gi++){
 630      equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
 631      equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
 632      equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
 633      equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
 634   }
 635
 636   // ******     Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
 637   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 638   for(Int_t gi=0; gi<5; gi++){
 639        calibSumZN1[0] += equalTowZN1[gi];
 640        calibSumZP1[0] += equalTowZP1[gi];
 641        calibSumZN2[0] += equalTowZN2[gi];
 642        calibSumZP2[0] += equalTowZP2[gi];
 643        //
 644        calibSumZN1[1] += equalTowZN1[gi+5];
 645        calibSumZP1[1] += equalTowZP1[gi+5];
 646        calibSumZN2[1] += equalTowZN2[gi+5];
 647        calibSumZP2[1] += equalTowZP2[gi+5];
 648   }
 649   // High gain chain
 650   calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
 651   calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
 652   calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
 653   calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
 654   // Low gain chain
 655   calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
 656   calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
 657   calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
 658   calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
 659
 660   // ******     Energy calibration of detector responses
 661   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 662   for(Int_t gi=0; gi<5; gi++){
 663      // High gain chain
 664      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
 665      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
 666      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
 667      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
 668      // Low gain chain
 669      calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
 670      calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
 671      calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
 672      calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
 673   }
 674   //
 675   Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
 676   calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
 677   calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
 678   calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
 679   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
 680   for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
 681
 682   //  ******    No. of spectator and participants nucleons
 683   //  Variables calculated to comply with ESD structure
 684   //  *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
 685   Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
 686   Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
 687   Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
 688   Double_t impPar=0., impPar1=0., impPar2=0.;
 689
 690   // create the output tree
 691   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 692                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
 693                   calibZEM1, calibZEM2, sPMRef1, sPMRef2,
 694                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 695                   nGenSpec, nGenSpecLeft, nGenSpecRight,
 696                   nPart, nPartTotLeft, nPartTotRight,
 697                   impPar, impPar1, impPar2);
 698
 699   AliZDCReco* preco = &reco;
 700   const Int_t kBufferSize = 4000;
 701   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
 702
 703   // write the output tree
 704   clustersTree->Fill();
 705 }
 706
 707 //_____________________________________________________________________________
 708 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
 709         Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
 710         Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2) const
 711 {
 712   // ****************** Reconstruct one event ******************
 713
 714   // ******     Retrieving calibration data
 715   // --- Equalization coefficients ---------------------------------------------
 716   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 717   for(Int_t ji=0; ji<5; ji++){
 718      equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
 719      equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
 720      equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
 721      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
 722   }
 723   // --- Energy calibration factors ------------------------------------
 724   Float_t valFromOCDB[6], calibEne[6];
 725   for(Int_t ij=0; ij<6; ij++){
 726     valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
 727     if(ij<4){
 728       if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
 729       else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
 730     }
 731     else calibEne[ij] = valFromOCDB[ij];
 732   }
 733
 734   // ******     Equalization of detector responses
 735   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 736   for(Int_t gi=0; gi<10; gi++){
 737      equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
 738      equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
 739      equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
 740      equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
 741   }
 742
 743   // ******     Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
 744   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 745   for(Int_t gi=0; gi<5; gi++){
 746        calibSumZN1[0] += equalTowZN1[gi];
 747        calibSumZP1[0] += equalTowZP1[gi];
 748        calibSumZN2[0] += equalTowZN2[gi];
 749        calibSumZP2[0] += equalTowZP2[gi];
 750        //
 751        calibSumZN1[1] += equalTowZN1[gi+5];
 752        calibSumZP1[1] += equalTowZP1[gi+5];
 753        calibSumZN2[1] += equalTowZN2[gi+5];
 754        calibSumZP2[1] += equalTowZP2[gi+5];
 755   }
 756   // High gain chain
 757   calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
 758   calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
 759   calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
 760   calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
 761   // Low gain chain
 762   calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
 763   calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
 764   calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
 765   calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
 766   //
 767   Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
 768   calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
 769   calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
 770   calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
 771   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
 772   for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
 773
 774   // ******     Energy calibration of detector responses
 775   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 776   for(Int_t gi=0; gi<5; gi++){
 777      // High gain chain
 778      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
 779      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
 780      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
 781      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
 782      // Low gain chain
 783      calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
 784      calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
 785      calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
 786      calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
 787   }
 788
 789   //  ******    Number of detected spectator nucleons
 790   Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
 791   if(fBeamEnergy!=0){
 792     nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
 793     nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
 794     nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
 795     nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
 796   }
 797   else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
 798   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
 799     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
 800     nDetSpecNRight, nDetSpecPRight);*/
 801
 802   if(fIsCalibrationMB == kFALSE){
 803     // ******   Reconstruction parameters ------------------
 804     // Ch. debug
 805     //fRecoParam->Print("");
 806     //
 807     TH2F *hZDCvsZEM  = fRecoParam->GethZDCvsZEM();
 808     TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
 809     TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
 810     TH1D *hNpartDist = fRecoParam->GethNpartDist();
 811     TH1D *hbDist = fRecoParam->GethbDist();
 812     Float_t ClkCenter = fRecoParam->GetClkCenter();
 813     //
 814     Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
 815     Double_t origin = xHighEdge*ClkCenter;
 816     // Ch. debug
 817     printf("\n\n  xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
 818     //
 819     // ====> Summed ZDC info (sideA+side C)
 820     TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
 821     Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
 822     Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
 823     line->SetParameter(0, y/(x-origin));
 824     line->SetParameter(1, -origin*y/(x-origin));
 825     // Ch. debug
 826     printf("  ***************** Summed ZDC info (sideA+side C) \n");
 827     printf("  E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f   ", x, y,y/(x-origin),-origin*y/(x-origin));
 828     //
 829     Double_t countPerc=0;
 830     Double_t xBinCenter=0, yBinCenter=0;
 831     for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
 832       for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
 833          xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
 834          yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
 835          //
 836          if(line->GetParameter(0)>0){
 837            if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
 838              countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
 839              // Ch. debug
 840              /*printf(" xBinCenter  %1.3f, yBinCenter %1.0f,  countPerc %1.0f\n",
 841                 xBinCenter, yBinCenter, countPerc);*/
 842            }
 843          }
 844          else{
 845            if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
 846              countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
 847              // Ch. debug
 848              /*printf(" xBinCenter  %1.3f, yBinCenter %1.0f,  countPerc %1.0f\n",
 849                 xBinCenter, yBinCenter, countPerc);*/
 850            }
 851          }
 852       }
 853     }
 854     //
 855     Double_t xSecPerc = 0.;
 856     if(hZDCvsZEM->GetEntries()!=0){
 857       xSecPerc = countPerc/hZDCvsZEM->GetEntries();
 858     }
 859     else{
 860       AliWarning("  Histogram hZDCvsZEM from OCDB has no entries!!!");
 861     }
 862     // Ch. debug
 863     //printf("  xSecPerc %1.4f  \n", xSecPerc);
 864
 865     // ====> side C
 866     TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
 867     Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
 868     lineC->SetParameter(0, yC/(x-origin));
 869     lineC->SetParameter(1, -origin*yC/(x-origin));
 870     // Ch. debug
 871     //printf("  ***************** Side C \n");
 872     //printf("  E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f   ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
 873     //
 874     Double_t countPercC=0;
 875     Double_t xBinCenterC=0, yBinCenterC=0;
 876     for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
 877       for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
 878          xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
 879          yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
 880          if(lineC->GetParameter(0)>0){
 881            if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
 882              countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
 883            }
 884          }
 885          else{
 886            if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
 887              countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
 888            }
 889          }
 890       }
 891     }
 892     //
 893     Double_t xSecPercC = 0.;
 894     if(hZDCCvsZEM->GetEntries()!=0){
 895       xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
 896     }
 897     else{
 898       AliWarning("  Histogram hZDCCvsZEM from OCDB has no entries!!!");
 899     }
 900     // Ch. debug
 901     //printf("  xSecPercC %1.4f  \n", xSecPercC);
 902
 903     // ====> side A
 904     TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
 905     Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
 906     lineA->SetParameter(0, yA/(x-origin));
 907     lineA->SetParameter(1, -origin*yA/(x-origin));
 908     //
 909     // Ch. debug
 910     //printf("  ***************** Side A \n");
 911     //printf("  E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f   ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
 912     //
 913     Double_t countPercA=0;
 914     Double_t xBinCenterA=0, yBinCenterA=0;
 915     for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
 916       for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
 917          xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
 918          yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
 919          if(lineA->GetParameter(0)>0){
 920            if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
 921              countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
 922            }
 923          }
 924          else{
 925            if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
 926              countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
 927            }
 928          }
 929       }
 930     }
 931     //
 932     Double_t xSecPercA = 0.;
 933     if(hZDCAvsZEM->GetEntries()!=0){
 934       xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
 935     }
 936     else{
 937       AliWarning("  Histogram hZDCAvsZEM from OCDB has no entries!!!");
 938     }
 939     // Ch. debug
 940     //printf("  xSecPercA %1.4f  \n", xSecPercA);
 941
 942     //  ******    Number of participants (from E_ZDC vs. E_ZEM correlation)
 943     Int_t nPart=0, nPartC=0, nPartA=0;
 944     Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
 945     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
 946       nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
 947       if((1.-nPartFrac) < xSecPerc){
 948         nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
 949         // Ch. debug
 950         //printf("  ***************** Summed ZDC info (sideA+side C) \n");
 951         //printf("  nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
 952         break;
 953       }
 954     }
 955     if(nPart<0) nPart=0;
 956     //
 957     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
 958       nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
 959       if((1.-nPartFracC) < xSecPercC){
 960         nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
 961         // Ch. debug
 962         //printf("  ***************** Side C \n");
 963         //printf("  nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
 964         break;
 965     }
 966     }
 967     if(nPartC<0) nPartC=0;
 968     //
 969     for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
 970       nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
 971       if((1.-nPartFracA) < xSecPercA){
 972         nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
 973         // Ch. debug
 974         //printf("  ***************** Side A \n");
 975         //printf("  nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
 976         break;
 977       }
 978     }
 979     if(nPartA<0) nPartA=0;
 980
 981     //  ******    Impact parameter (from E_ZDC vs. E_ZEM correlation)
 982     Float_t b=0, bC=0, bA=0;
 983     Double_t bFrac=0., bFracC=0., bFracA=0.;
 984     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
 985       bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
 986       if((1.-bFrac) < xSecPerc){
 987         b = hbDist->GetBinLowEdge(ibbin);
 988         break;
 989       }
 990     }
 991     //
 992     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
 993       bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
 994       if((1.-bFracC) < xSecPercC){
 995         bC = hbDist->GetBinLowEdge(ibbin);
 996         break;
 997       }
 998     }
 999     //
1000     for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1001       bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1002       if((1.-bFracA) < xSecPercA){
1003         bA = hbDist->GetBinLowEdge(ibbin);
1004         break;
1005       }
1006     }
1007
1008     //  ******  Number of spectator nucleons
1009     Int_t nGenSpec=0, nGenSpecC=0, nGenSpecA=0;
1010     //
1011     nGenSpec = 416 - nPart;
1012     nGenSpecC = 416 - nPartC;
1013     nGenSpecA = 416 - nPartA;
1014     if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1015     if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1016     if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1017
1018     //  Ch. debug
1019     /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1020         "  calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
1021         calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1022     printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n",
1023         nGenSpecLeft, nGenSpecRight);
1024     printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1025     */
1026
1027     // create the output tree
1028     AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1029                     calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1030                     calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1031                     nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1032                     nGenSpec, nGenSpecA, nGenSpecC,
1033                     nPart, nPartA, nPartC, b, bA, bC);
1034
1035     AliZDCReco* preco = &reco;
1036     const Int_t kBufferSize = 4000;
1037     clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
1038     // write the output tree
1039     clustersTree->Fill();
1040
1041     delete lineC;  delete lineA;
1042   } // ONLY IF fIsCalibrationMB==kFALSE
1043   else{
1044     // create the output tree
1045     AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1046                     calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1047                     calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1048                     nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1049                     0, 0, 0,
1050                     0, 0, 0, 0., 0., 0.);
1051
1052     AliZDCReco* preco = &reco;
1053     const Int_t kBufferSize = 4000;
1054     clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
1055     // write the output tree
1056     clustersTree->Fill();
1057   }
1058 }
1059
1060 //_____________________________________________________________________________
1061 void AliZDCReconstructor::BuildRecoParam(TH2F* hCorr, TH2F* hCorrC, TH2F* hCorrA,
1062                         Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1063 {
1064   // Calculate RecoParam object for Pb-Pb data
1065   hCorr->Fill(ZDCC+ZDCA, ZEM);
1066   hCorrC->Fill(ZDCC, ZEM);
1067   hCorrA->Fill(ZDCA, ZEM);
1068   //
1069   /*TH1D*       hNpartDist;
1070   TH1D* hbDist;
1071   Float_t clkCenter;*/
1072
1073 }
1074
1075 //_____________________________________________________________________________
1076 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1077 {
1078   // fill energies and number of participants to the ESD
1079
1080   if(fIsCalibrationMB==kTRUE) WritePbPbRecoParamInOCDB();
1081
1082   AliZDCReco reco;
1083   AliZDCReco* preco = &reco;
1084   clustersTree->SetBranchAddress("ZDC", &preco);
1085
1086   clustersTree->GetEntry(0);
1087   //
1088   AliESDZDC * esdzdc = esd->GetESDZDC();
1089   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1090   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1091   for(Int_t i=0; i<5; i++){
1092      tZN1Ene[i] = reco.GetZN1HREnTow(i);
1093      tZN2Ene[i] = reco.GetZN2HREnTow(i);
1094      tZP1Ene[i] = reco.GetZP1HREnTow(i);
1095      tZP2Ene[i] = reco.GetZP2HREnTow(i);
1096      //
1097      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1098      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1099      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1100      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1101   }
1102   //
1103   esdzdc->SetZN1TowerEnergy(tZN1Ene);
1104   esdzdc->SetZN2TowerEnergy(tZN2Ene);
1105   esdzdc->SetZP1TowerEnergy(tZP1Ene);
1106   esdzdc->SetZP2TowerEnergy(tZP2Ene);
1107   //
1108   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1109   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1110   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1111   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1112   //
1113   Int_t nPart  = reco.GetNParticipants();
1114   Int_t nPartA = reco.GetNPartSideA();
1115   Int_t nPartC = reco.GetNPartSideC();
1116   Double_t b  = reco.GetImpParameter();
1117   Double_t bA = reco.GetImpParSideA();
1118   Double_t bC = reco.GetImpParSideC();
1119  //
1120   esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
1121               reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1122               nPart, nPartA, nPartC, b, bA, bC, fRecoFlag);
1123
1124 }
1125
1126 //_____________________________________________________________________________
1127 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1128 {
1129   // Setting the storage
1130
1131   Bool_t deleteManager = kFALSE;
1132
1133   AliCDBManager *manager = AliCDBManager::Instance();
1134   AliCDBStorage *defstorage = manager->GetDefaultStorage();
1135
1136   if(!defstorage || !(defstorage->Contains("ZDC"))){
1137      AliWarning("No default storage set or default storage doesn't contain ZDC!");
1138      manager->SetDefaultStorage(uri);
1139      deleteManager = kTRUE;
1140   }
1141
1142   AliCDBStorage *storage = manager->GetDefaultStorage();
1143
1144   if(deleteManager){
1145     AliCDBManager::Instance()->UnsetDefaultStorage();
1146     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
1147   }
1148
1149   return storage;
1150 }
1151
1152 //_____________________________________________________________________________
1153 AliZDCPedestals* AliZDCReconstructor::GetPedData() const
1154 {
1155
1156   // Getting pedestal calibration object for ZDC set
1157
1158   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1159   if(!entry) AliFatal("No calibration data loaded!");
1160
1161   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
1162   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1163
1164   return calibdata;
1165 }
1166
1167 //_____________________________________________________________________________
1168 AliZDCEnCalib* AliZDCReconstructor::GetEnCalibData() const
1169 {
1170
1171   // Getting energy and equalization calibration object for ZDC set
1172
1173   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnCalib");
1174   if(!entry) AliFatal("No calibration data loaded!");
1175
1176   AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1177   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1178
1179   return calibdata;
1180 }
1181
1182 //_____________________________________________________________________________
1183 AliZDCTowerCalib* AliZDCReconstructor::GetTowCalibData() const
1184 {
1185
1186   // Getting energy and equalization calibration object for ZDC set
1187
1188   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowCalib");
1189   if(!entry) AliFatal("No calibration data loaded!");
1190
1191   AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1192   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
1193
1194   return calibdata;
1195 }
1196
1197 //_____________________________________________________________________________
1198 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1199 {
1200
1201   // Getting reconstruction parameters from OCDB
1202
1203   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1204   if(!entry) AliFatal("No RecoParam data found in OCDB!");
1205
1206   AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1207   if(!param)  AliFatal("No RecoParam object in OCDB entry!");
1208
1209   return param;
1210
1211 }
1212
1213 //_____________________________________________________________________________
1214 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1215 {
1216
1217   // Getting reconstruction parameters from OCDB
1218
1219   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1220   if(!entry) AliFatal("No RecoParam data found in OCDB!");
1221
1222   AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1223   if(!param)  AliFatal("No RecoParam object in OCDB entry!");
1224
1225   return param;
1226
1227 }
1228
1229 //_____________________________________________________________________________
1230 void AliZDCReconstructor::WritePbPbRecoParamInOCDB() const
1231 {
1232
1233   // Writing Pb-Pb reconstruction parameters from OCDB
1234
1235   AliCDBManager *man = AliCDBManager::Instance();
1236   AliCDBMetaData *md= new AliCDBMetaData();
1237   md->SetResponsible("Chiara Oppedisano");
1238   md->SetComment("ZDC Pb-Pb reconstruction parameters");
1239   md->SetObjectClassName("AliZDCRecoParamPbPb");
1240   AliCDBId id("ZDC/Calib/RecoParamPbPb",fNRun,AliCDBRunRange::Infinity());
1241   man->Put(fRecoParam, id, md);
1242
1243 }
1244