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 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24
  25 #include <TF1.h>
  26
  27 #include "AliRunLoader.h"
  28 #include "AliRunInfo.h"
  29 #include "AliRawReader.h"
  30 #include "AliESDEvent.h"
  31 #include "AliESDZDC.h"
  32 #include "AliZDCDigit.h"
  33 #include "AliZDCRawStream.h"
  34 #include "AliZDCReco.h"
  35 #include "AliZDCReconstructor.h"
  36 #include "AliZDCPedestals.h"
  37 #include "AliZDCCalib.h"
  38 #include "AliZDCRecoParam.h"
  39 #include "AliZDCRecoParampp.h"
  40 #include "AliZDCRecoParamPbPb.h"
  41
  42
  43 ClassImp(AliZDCReconstructor)
  44 AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0;  //reconstruction parameters
  45
  46 //_____________________________________________________________________________
  47 AliZDCReconstructor:: AliZDCReconstructor() :
  48   fPedData(GetPedData()),
  49   fECalibData(GetECalibData())
  50 {
  51   // **** Default constructor
  52
  53 }
  54
  55
  56 //_____________________________________________________________________________
  57 AliZDCReconstructor::~AliZDCReconstructor()
  58 {
  59 // destructor
  60    if(fRecoParam)  delete fRecoParam;
  61    if(fPedData)    delete fPedData;
  62    if(fECalibData) delete fECalibData;
  63 }
  64
  65 //_____________________________________________________________________________
  66 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
  67 {
  68   // *** Local ZDC reconstruction for digits
  69   // Works on the current event
  70
  71   // Retrieving calibration data
  72   Float_t meanPed[48];
  73   for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
  74
  75   // get digits
  76   AliZDCDigit digit;
  77   AliZDCDigit* pdigit = &digit;
  78   digitsTree->SetBranchAddress("ZDC", &pdigit);
  79   //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
  80
  81   // loop over digits
  82   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
  83   Float_t dZEM1Corr[2], dZEM2Corr[2], PMRef1[2], PMRef2[2];
  84   for(Int_t i=0; i<10; i++){
  85      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
  86      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = PMRef1[i] = PMRef2[i] = 0.;
  87   }
  88   //
  89   for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
  90    digitsTree->GetEntry(iDigit);
  91    if (!pdigit) continue;
  92    //
  93    Int_t det = digit.GetSector(0);
  94    Int_t quad = digit.GetSector(1);
  95    Int_t pedindex = -1, kNch = 24;
  96    //printf("\n\t Digit #%d det %d quad %d", iDigit, det, quad);
  97    //
  98    if(quad != 5){ // ZDC (not reference PTMs!)
  99     if(det == 1){ // *** ZNC
 100        pedindex = quad;
 101        tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 102        if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
 103        tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 104        if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
 105        //printf("\t pedindex %d tZN1Corr[%d] = %1.0f tZN1Corr[%d] = %1.0f",
 106        //       pedindex, quad, tZN1Corr[quad], quad+5, tZN1Corr[quad+5]);
 107     }
 108     else if(det == 2){ // *** ZP1
 109        pedindex = quad+5;
 110        tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 111        if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
 112        tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 113        if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
 114        //printf("\t pedindex %d tZP1Corr[%d] = %1.0f tZP1Corr[%d] = %1.0f",
 115        //       pedindex, quad, tZP1Corr[quad], quad+5, tZP1Corr[quad+5]);
 116     }
 117     else if(det == 3){
 118        pedindex = quad+9;
 119        if(quad == 1){       // *** ZEM1
 120          dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 121          if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
 122          dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 123          if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
 124          //printf("\t pedindex %d tZEM1Corr[%d] = %1.0f tZEM1Corr[%d] = %1.0f",
 125          //     pedindex, quad, tZEM1Corr[quad], quad+1, tZEM1Corr[quad+1]);
 126        }
 127        else if(quad == 2){  // *** ZEM2
 128          dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 129          if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
 130          dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 131          if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
 132          //printf("\t pedindex %d tZEM2Corr[%d] = %1.0f tZEM2Corr[%d] = %1.0f",
 133          //     pedindex, quad, tZEM2Corr[quad], quad+1, tZEM2Corr[quad+1]);
 134        }
 135     }
 136     else if(det == 4){  // *** ZN2
 137        pedindex = quad+12;
 138        tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 139        if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
 140        tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 141        if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
 142        //printf("\t pedindex %d tZN2Corr[%d] = %1.0f tZN2Corr[%d] = %1.0f\n",
 143        //       pedindex, quad, tZN2Corr[quad], quad+5, tZN2Corr[quad+5]);
 144     }
 145     else if(det == 5){  // *** ZP2
 146        pedindex = quad+17;
 147        tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 148        if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
 149        tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 150        if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
 151        //printf("\t pedindex %d tZP2Corr[%d] = %1.0f tZP2Corr[%d] = %1.0f\n",
 152        //       pedindex, quad, tZP2Corr[quad], quad+5, tZP2Corr[quad+5]);
 153     }
 154    }
 155    else{ // Reference PMs
 156      pedindex = (det-1)/3+22;
 157      if(det == 1){
 158        PMRef1[0] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 159        if(PMRef1[0]<0.) PMRef1[0] = 0.;
 160        PMRef1[1] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 161        if(PMRef2[1]<0.) PMRef1[1] = 0.;
 162      }
 163      else if(det == 4){
 164        PMRef2[0] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 165        if(PMRef2[0]<0.) PMRef2[0] = 0.;
 166        PMRef2[1] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 167        if(PMRef2[1]<0.) PMRef2[1] = 0.;
 168      }
 169    }
 170   }
 171
 172   // Reconstruct the event according to beam type from AliInfo
 173   AliRunInfo *runInfo;
 174   if((strcmp(runInfo->GetBeamType(),"p-p") == 0)){
 175     printf("\n\t ZDC reconstructing p-p event\n");
 176     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 177         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
 178   }
 179   else if((strcmp(runInfo->GetBeamType(),"A-A") == 0)){
 180     printf("\n\t ZDC reconstructing A-A event\n");
 181     ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 182         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
 183   }
 184   else printf("\n\t ZDC does not reconstruct event! Beam type != p-p || A-A\n");
 185 }
 186
 187 //_____________________________________________________________________________
 188 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
 189 {
 190   // *** ZDC raw data reconstruction
 191   // Works on the current event
 192
 193   // Retrieving calibration data
 194   Float_t meanPed[48];
 195   for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
 196
 197   rawReader->Reset();
 198
 199   // loop over raw data
 200   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
 201   Float_t dZEM1Corr[2], dZEM2Corr[2], PMRef1[2], PMRef2[2];
 202   for(Int_t i=0; i<10; i++){
 203      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
 204      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = PMRef1[i] = PMRef2[i] = 0.;
 205   }
 206   //
 207   AliZDCRawStream rawData(rawReader);
 208   Int_t kNch = 24;
 209   while (rawData.Next()) {
 210     if(rawData.IsADCDataWord()){
 211      Int_t det = rawData.GetSector(0);
 212      Int_t quad = rawData.GetSector(1);
 213      Int_t gain = rawData.GetADCGain();
 214      Int_t pedindex=0;
 215      //
 216      if(quad !=5){ // ZDCs (not reference PTMs)
 217       if(det == 1){
 218         pedindex = quad;
 219         if(gain == 0) tZN1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 220         else tZN1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 221       }
 222       else if(det == 2){
 223         pedindex = quad+5;
 224         if(gain == 0) tZP1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 225         else tZP1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 226       }
 227       else if(det == 3){
 228         pedindex = quad+9;
 229         if(quad==1){
 230           if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 231           else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 232         }
 233         else if(quad==2){
 234           if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 235           else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 236         }
 237       }
 238       else if(det == 4){
 239         pedindex = quad+12;
 240         if(gain == 0) tZN2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 241         else tZN2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 242       }
 243       else if(det == 5){
 244         pedindex = quad+17;
 245         if(gain == 0) tZP2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 246         else tZP2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 247       }
 248       //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
 249       //  det,quad,gain, pedindex, meanPed[pedindex]);
 250      }
 251      else{ // reference PM
 252        pedindex = (det-1)/3 + 22;
 253        if(det == 1){
 254          if(gain==0) PMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 255          else PMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 256        }
 257        else if(det ==4){
 258          if(gain==0) PMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 259          else PMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 260        }
 261      }
 262     }//IsADCDataWord
 263   }
 264
 265   // Reconstruct the event according to beam type from AliInfo
 266   AliRunInfo *runInfo;
 267   if((strcmp(runInfo->GetBeamType(),"p-p") == 0)){
 268     printf("\n\t ZDC reconstructing p-p event\n");
 269     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 270         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
 271   }
 272   else if((strcmp(runInfo->GetBeamType(),"A-A") == 0)){
 273     printf("\n\t ZDC reconstructing A-A event\n");
 274     ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 275         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
 276   }
 277   else printf("\n\t ZDC does not reconstruct event! Beam type != p-p || A-A\n");
 278
 279 }
 280
 281 //_____________________________________________________________________________
 282 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* ZN1ADCCorr,
 283         Float_t* ZP1ADCCorr, Float_t* ZN2ADCCorr, Float_t* ZP2ADCCorr,
 284         Float_t* ZEM1ADCCorr, Float_t* ZEM2ADCCorr, Float_t* PMRef1, Float_t* PMRef2) const
 285 {
 286   // ***** Reconstruct one event
 287
 288   // *** RECONSTRUCTION FROM "REAL" DATA
 289   //
 290   // Retrieving calibration data
 291   // --- Equalization coefficients ---------------------------------------------
 292   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 293   for(Int_t ji=0; ji<5; ji++){
 294      equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
 295      equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
 296      equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
 297      equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
 298   }
 299   // --- Energy calibration factors ------------------------------------
 300   Float_t calibEne[4];
 301   // *********************************************************************
 302   // **** Until the beam type info isn't known @ reconstruction level ****
 303   // **** the energy calibration coefficient are manually set to 1    ****
 304   // **** as it will be in real life for pp data taking               ****
 305   // *********************************************************************
 306   //for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
 307   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = 1.;
 308
 309   // Equalization of detector responses
 310   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 311   for(Int_t gi=0; gi<5; gi++){
 312      equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
 313      equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
 314      equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
 315      equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
 316      equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
 317      equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
 318      equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
 319      equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
 320   }
 321
 322   // Energy calibration of detector responses
 323   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 324   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 325   for(Int_t gi=0; gi<10; gi++){
 326      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
 327      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
 328      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
 329      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
 330      //
 331      if(gi<5){
 332        calibSumZN1[0] += calibTowZN1[gi];
 333        calibSumZP1[0] += calibTowZP1[gi];
 334        calibSumZN2[0] += calibTowZN2[gi];
 335        calibSumZP2[0] += calibTowZP2[gi];
 336      }
 337      //
 338      else{
 339        calibSumZN1[1] += calibTowZN1[gi];
 340        calibSumZP1[1] += calibTowZP1[gi];
 341        calibSumZN2[1] += calibTowZN2[gi];
 342        calibSumZP2[1] += calibTowZP2[gi];
 343      }
 344   }
 345
 346   //
 347   // --- Reconstruction parameters ------------------
 348   if(!fRecoParam)  fRecoParam = (AliZDCRecoParampp*) AliZDCRecoParampp::GetppRecoParam();
 349
 350   //  ---      Number of detected spectator nucleons
 351   //  *** N.B. -> It works only in Pb-Pb
 352   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
 353   Float_t beamE = fRecoParam->GetBeamEnergy();
 354   nDetSpecNLeft = (Int_t) (calibSumZN1[0]/beamE);
 355   nDetSpecPLeft = (Int_t) (calibSumZP1[0]/beamE);
 356   nDetSpecNRight = (Int_t) (calibSumZN2[0]/beamE);
 357   nDetSpecPRight = (Int_t) (calibSumZP2[0]/beamE);
 358   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
 359     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
 360     nDetSpecNRight, nDetSpecPRight);*/
 361
 362   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
 363   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
 364   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
 365   Int_t nPartTotLeft=0, nPartTotRight=0;
 366   Double_t impPar=0.;
 367
 368   // create the output tree
 369   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 370                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
 371                   ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
 372                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 373                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
 374                   nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
 375
 376   AliZDCReco* preco = &reco;
 377   const Int_t kBufferSize = 4000;
 378   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
 379
 380   // write the output tree
 381   clustersTree->Fill();
 382 }
 383
 384 //_____________________________________________________________________________
 385 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree, Float_t* ZN1ADCCorr,
 386         Float_t* ZP1ADCCorr, Float_t* ZN2ADCCorr, Float_t* ZP2ADCCorr,
 387         Float_t* ZEM1ADCCorr, Float_t* ZEM2ADCCorr, Float_t* PMRef1, Float_t* PMRef2) const
 388 {
 389   // ***** Reconstruct one event
 390
 391   // *** RECONSTRUCTION FROM "REAL" DATA
 392   //
 393   // Retrieving calibration data
 394   // --- Equalization coefficients ---------------------------------------------
 395   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 396   for(Int_t ji=0; ji<5; ji++){
 397      equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
 398      equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
 399      equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
 400      equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
 401   }
 402   // --- Energy calibration factors ------------------------------------
 403   Float_t calibEne[4];
 404   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
 405
 406   // Equalization of detector responses
 407   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 408   for(Int_t gi=0; gi<5; gi++){
 409      equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
 410      equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
 411      equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
 412      equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
 413      equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
 414      equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
 415      equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
 416      equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
 417   }
 418
 419   // Energy calibration of detector responses
 420   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 421   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 422   for(Int_t gi=0; gi<10; gi++){
 423      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
 424      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
 425      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
 426      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
 427      //
 428      if(gi<5){
 429        calibSumZN1[0] += calibTowZN1[gi];
 430        calibSumZP1[0] += calibTowZP1[gi];
 431        calibSumZN2[0] += calibTowZN2[gi];
 432        calibSumZP2[0] += calibTowZP2[gi];
 433      }
 434      //
 435      else{
 436        calibSumZN1[1] += calibTowZN1[gi];
 437        calibSumZP1[1] += calibTowZP1[gi];
 438        calibSumZN2[1] += calibTowZN2[gi];
 439        calibSumZP2[1] += calibTowZP2[gi];
 440      }
 441   }
 442
 443   //
 444   // --- Reconstruction parameters ------------------
 445   if(!fRecoParam)  fRecoParam = (AliZDCRecoParamPbPb*) AliZDCRecoParamPbPb::GetPbPbRecoParam();
 446   //
 447   Float_t endPointZEM = fRecoParam->GetZEMEndValue();
 448   Float_t cutFractionZEM = fRecoParam->GetZEMCutFraction();
 449   Float_t dZEMSup = fRecoParam->GetDZEMSup();
 450   Float_t dZEMInf = fRecoParam->GetDZEMInf();
 451   //
 452   Float_t cutValueZEM = endPointZEM*cutFractionZEM;
 453   Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
 454   Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
 455   //
 456   Float_t maxValEZN1  = fRecoParam->GetEZN1MaxValue();
 457   Float_t maxValEZP1  = fRecoParam->GetEZP1MaxValue();
 458   Float_t maxValEZDC1 = fRecoParam->GetEZDC1MaxValue();
 459   Float_t maxValEZN2  = fRecoParam->GetEZN2MaxValue();
 460   Float_t maxValEZP2  = fRecoParam->GetEZP2MaxValue();
 461   Float_t maxValEZDC2 = fRecoParam->GetEZDC2MaxValue();
 462   //
 463   //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
 464   //   " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
 465
 466   //  ---      Number of detected spectator nucleons
 467   //  *** N.B. -> It works only in Pb-Pb
 468   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
 469   Float_t beamE = fRecoParam->GetBeamEnergy();
 470   nDetSpecNLeft = (Int_t) (calibSumZN1[0]/beamE);
 471   nDetSpecPLeft = (Int_t) (calibSumZP1[0]/beamE);
 472   nDetSpecNRight = (Int_t) (calibSumZN2[0]/beamE);
 473   nDetSpecPRight = (Int_t) (calibSumZP2[0]/beamE);
 474   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
 475     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
 476     nDetSpecNRight, nDetSpecPRight);*/
 477
 478   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
 479   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
 480   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
 481   Double_t impPar=0.;
 482   //
 483   Float_t corrADCZEMHG = ZEM1ADCCorr[0] + ZEM2ADCCorr[0];
 484   //
 485   if(corrADCZEMHG > supValueZEM){
 486     nGenSpecNLeft  = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN1[0]));
 487     nGenSpecPLeft  = (Int_t) ((fRecoParam->GetfZPCen())->Eval(calibSumZP1[0]));
 488     nGenSpecLeft   = (Int_t) ((fRecoParam->GetfZDCCen())->Eval(calibSumZN1[0]+calibSumZP1[0]));
 489     nGenSpecNRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN2[0]));
 490     nGenSpecPRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZP2[0]));
 491     nGenSpecRight  = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN2[0]+calibSumZP2[0]));
 492     impPar  = (fRecoParam->GetfbCen())->Eval(calibSumZN1[0]+calibSumZP1[0]);
 493   }
 494   else if(corrADCZEMHG < infValueZEM){
 495     nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZNPer())->Eval(calibSumZN1[0]));
 496     nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZPPer())->Eval(calibSumZP1[0]));
 497     nGenSpecLeft  = (Int_t) ((fRecoParam->GetfZDCPer())->Eval(calibSumZN1[0]+calibSumZP1[0]));
 498     impPar   = (fRecoParam->GetfbPer())->Eval(calibSumZN1[0]+calibSumZP1[0]);
 499   }
 500   else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
 501     nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 502     nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 503     nGenSpecLeft  = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 504     impPar   =  (fRecoParam->GetfZEMb())->Eval(corrADCZEMHG);
 505   }
 506   //
 507   if(calibSumZN1[0]/maxValEZN1>1.)  nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 508   if(calibSumZP1[0]/maxValEZP1>1.)  nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 509   if((calibSumZN1[0]+calibSumZP1[0]/maxValEZDC1)>1.){
 510      nGenSpecLeft = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 511      impPar = (fRecoParam->GetfZEMb())->Eval(corrADCZEMHG);
 512   }
 513   if(calibSumZN2[0]/maxValEZN2>1.)  nGenSpecNRight = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 514   if(calibSumZP2[0]/maxValEZP2>1.)  nGenSpecPRight = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 515   if((calibSumZN2[0]+calibSumZP2[0]/maxValEZDC2)>1.) nGenSpecRight = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 516   //
 517   if(nGenSpecNLeft>125)    nGenSpecNLeft=125;
 518   else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
 519   if(nGenSpecPLeft>82)     nGenSpecPLeft=82;
 520   else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
 521   if(nGenSpecLeft>207)     nGenSpecLeft=207;
 522   else if(nGenSpecLeft<0)  nGenSpecLeft=0;
 523
 524   //  ---      Number of generated participants (from HIJING parameterization)
 525   Int_t nPart, nPartTotLeft, nPartTotRight;
 526   nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
 527   nPartTotLeft = 207-nGenSpecLeft;
 528   nPartTotRight = 207-nGenSpecRight;
 529   if(nPart<0) nPart=0;
 530   if(nPartTotLeft<0) nPartTotLeft=0;
 531   if(nPartTotRight<0) nPartTotRight=0;
 532   //
 533   // *** DEBUG ***
 534   /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
 535       "  calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
 536       calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
 537   printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
 538       "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n",
 539       nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft,
 540       nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
 541   printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
 542   */
 543
 544   // create the output tree
 545   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 546                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
 547                   ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
 548                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 549                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
 550                   nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
 551
 552   AliZDCReco* preco = &reco;
 553   const Int_t kBufferSize = 4000;
 554   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
 555
 556   // write the output tree
 557   clustersTree->Fill();
 558 }
 559
 560 //_____________________________________________________________________________
 561 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
 562 {
 563   // fill energies and number of participants to the ESD
 564
 565   AliZDCReco reco;
 566   AliZDCReco* preco = &reco;
 567   clustersTree->SetBranchAddress("ZDC", &preco);
 568
 569   clustersTree->GetEntry(0);
 570   //
 571   AliESDZDC * esdzdc = esd->GetESDZDC();
 572   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
 573   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
 574   for(Int_t i=0; i<5; i++){
 575      tZN1Ene[i] = reco.GetZN1HREnTow(i);
 576      tZN2Ene[i] = reco.GetZN2HREnTow(i);
 577      tZP1Ene[i] = reco.GetZP1HREnTow(i);
 578      tZP2Ene[i] = reco.GetZP2HREnTow(i);
 579      //
 580      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
 581      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
 582      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
 583      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
 584   }
 585   esdzdc->SetZN1TowerEnergy(tZN1Ene);
 586   esdzdc->SetZN2TowerEnergy(tZN2Ene);
 587   esdzdc->SetZP1TowerEnergy(tZP1Ene);
 588   esdzdc->SetZP2TowerEnergy(tZP2Ene);
 589   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
 590   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
 591   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
 592   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
 593   //
 594   esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
 595               reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
 596               reco.GetNPartLeft(),reco.GetNPartRight());
 597   //
 598
 599 }
 600
 601 //_____________________________________________________________________________
 602 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
 603 {
 604   // Setting the storage
 605
 606   Bool_t deleteManager = kFALSE;
 607
 608   AliCDBManager *manager = AliCDBManager::Instance();
 609   AliCDBStorage *defstorage = manager->GetDefaultStorage();
 610
 611   if(!defstorage || !(defstorage->Contains("ZDC"))){
 612      AliWarning("No default storage set or default storage doesn't contain ZDC!");
 613      manager->SetDefaultStorage(uri);
 614      deleteManager = kTRUE;
 615   }
 616
 617   AliCDBStorage *storage = manager->GetDefaultStorage();
 618
 619   if(deleteManager){
 620     AliCDBManager::Instance()->UnsetDefaultStorage();
 621     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
 622   }
 623
 624   return storage;
 625 }
 626
 627 //_____________________________________________________________________________
 628 AliZDCPedestals* AliZDCReconstructor::GetPedData() const
 629 {
 630
 631   // Getting pedestal calibration object for ZDC set
 632
 633   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
 634   if(!entry) AliFatal("No calibration data loaded!");
 635
 636   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
 637   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 638
 639   return calibdata;
 640 }
 641
 642 //_____________________________________________________________________________
 643 AliZDCCalib* AliZDCReconstructor::GetECalibData() const
 644 {
 645
 646   // Getting energy and equalization calibration object for ZDC set
 647
 648   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EMDCalib");
 649   if(!entry) AliFatal("No calibration data loaded!");
 650
 651   AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*>  (entry->GetObject());
 652   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 653
 654   return calibdata;
 655 }
 656