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