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] = 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] = 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   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
 282
 283   // Equalization of detector responses
 284   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 285   for(Int_t gi=0; gi<5; gi++){
 286      equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
 287      equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
 288      equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
 289      equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
 290      equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
 291      equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
 292      equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
 293      equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
 294   }
 295
 296   // Energy calibration of detector responses
 297   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 298   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 299   for(Int_t gi=0; gi<10; gi++){
 300      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
 301      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
 302      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
 303      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
 304      //
 305      if(gi<5){
 306        calibSumZN1[0] += calibTowZN1[gi];
 307        calibSumZP1[0] += calibTowZP1[gi];
 308        calibSumZN2[0] += calibTowZN2[gi];
 309        calibSumZP2[0] += calibTowZP2[gi];
 310      }
 311      //
 312      else{
 313        calibSumZN1[1] += calibTowZN1[gi];
 314        calibSumZP1[1] += calibTowZP1[gi];
 315        calibSumZN2[1] += calibTowZN2[gi];
 316        calibSumZP2[1] += calibTowZP2[gi];
 317      }
 318   }
 319
 320   //
 321   // --- Reconstruction parameters ------------------
 322   if(!fRecoParam)  fRecoParam = (AliZDCRecoParampp*) AliZDCRecoParampp::GetppRecoParam();
 323
 324   //  ---      Number of detected spectator nucleons
 325   //  *** N.B. -> It works only in Pb-Pb
 326   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
 327   Float_t beamE = fRecoParam->GetBeamEnergy();
 328   nDetSpecNLeft = (Int_t) (calibSumZN1[0]/beamE);
 329   nDetSpecPLeft = (Int_t) (calibSumZP1[0]/beamE);
 330   nDetSpecNRight = (Int_t) (calibSumZN2[0]/beamE);
 331   nDetSpecPRight = (Int_t) (calibSumZP2[0]/beamE);
 332   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
 333     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
 334     nDetSpecNRight, nDetSpecPRight);*/
 335
 336   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
 337   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
 338   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
 339   Int_t nPartTotLeft=0, nPartTotRight=0;
 340   Double_t impPar=0.;
 341
 342   // create the output tree
 343   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 344                   calibTowZN1, calibTowZN2, calibTowZP1, calibTowZP2,
 345                   ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
 346                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 347                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
 348                   nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
 349
 350   AliZDCReco* preco = &reco;
 351   const Int_t kBufferSize = 4000;
 352   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
 353
 354   // write the output tree
 355   clustersTree->Fill();
 356 }
 357
 358 //_____________________________________________________________________________
 359 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree, Float_t* ZN1ADCCorr,
 360         Float_t* ZP1ADCCorr, Float_t* ZN2ADCCorr, Float_t* ZP2ADCCorr,
 361         Float_t* ZEM1ADCCorr, Float_t* ZEM2ADCCorr, Float_t* PMRef1, Float_t* PMRef2) const
 362 {
 363   // ***** Reconstruct one event
 364
 365   // *** RECONSTRUCTION FROM "REAL" DATA
 366   //
 367   // Retrieving calibration data
 368   // --- Equalization coefficients ---------------------------------------------
 369   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 370   for(Int_t ji=0; ji<5; ji++){
 371      equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
 372      equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
 373      equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
 374      equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
 375   }
 376   // --- Energy calibration factors ------------------------------------
 377   Float_t calibEne[4];
 378   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
 379
 380   // Equalization of detector responses
 381   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 382   for(Int_t gi=0; gi<5; gi++){
 383      equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
 384      equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
 385      equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
 386      equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
 387      equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
 388      equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
 389      equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
 390      equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
 391   }
 392
 393   // Energy calibration of detector responses
 394   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 395   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 396   for(Int_t gi=0; gi<10; gi++){
 397      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
 398      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
 399      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
 400      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
 401      //
 402      if(gi<5){
 403        calibSumZN1[0] += calibTowZN1[gi];
 404        calibSumZP1[0] += calibTowZP1[gi];
 405        calibSumZN2[0] += calibTowZN2[gi];
 406        calibSumZP2[0] += calibTowZP2[gi];
 407      }
 408      //
 409      else{
 410        calibSumZN1[1] += calibTowZN1[gi];
 411        calibSumZP1[1] += calibTowZP1[gi];
 412        calibSumZN2[1] += calibTowZN2[gi];
 413        calibSumZP2[1] += calibTowZP2[gi];
 414      }
 415   }
 416
 417   //
 418   // --- Reconstruction parameters ------------------
 419   if(!fRecoParam)  fRecoParam = (AliZDCRecoParamPbPb*) AliZDCRecoParamPbPb::GetPbPbRecoParam();
 420   //
 421   Float_t endPointZEM = fRecoParam->GetZEMEndValue();
 422   Float_t cutFractionZEM = fRecoParam->GetZEMCutFraction();
 423   Float_t dZEMSup = fRecoParam->GetDZEMSup();
 424   Float_t dZEMInf = fRecoParam->GetDZEMInf();
 425   //
 426   Float_t cutValueZEM = endPointZEM*cutFractionZEM;
 427   Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
 428   Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
 429   //
 430   Float_t maxValEZN1  = fRecoParam->GetEZN1MaxValue();
 431   Float_t maxValEZP1  = fRecoParam->GetEZP1MaxValue();
 432   Float_t maxValEZDC1 = fRecoParam->GetEZDC1MaxValue();
 433   Float_t maxValEZN2  = fRecoParam->GetEZN2MaxValue();
 434   Float_t maxValEZP2  = fRecoParam->GetEZP2MaxValue();
 435   Float_t maxValEZDC2 = fRecoParam->GetEZDC2MaxValue();
 436   //
 437   //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
 438   //   " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
 439
 440   //  ---      Number of detected spectator nucleons
 441   //  *** N.B. -> It works only in Pb-Pb
 442   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
 443   Float_t beamE = fRecoParam->GetBeamEnergy();
 444   nDetSpecNLeft = (Int_t) (calibSumZN1[0]/beamE);
 445   nDetSpecPLeft = (Int_t) (calibSumZP1[0]/beamE);
 446   nDetSpecNRight = (Int_t) (calibSumZN2[0]/beamE);
 447   nDetSpecPRight = (Int_t) (calibSumZP2[0]/beamE);
 448   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
 449     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
 450     nDetSpecNRight, nDetSpecPRight);*/
 451
 452   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
 453   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
 454   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
 455   Double_t impPar=0.;
 456   //
 457   Float_t corrADCZEMHG = ZEM1ADCCorr[0] + ZEM2ADCCorr[0];
 458   //
 459   if(corrADCZEMHG > supValueZEM){
 460     nGenSpecNLeft  = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN1[0]));
 461     nGenSpecPLeft  = (Int_t) ((fRecoParam->GetfZPCen())->Eval(calibSumZP1[0]));
 462     nGenSpecLeft   = (Int_t) ((fRecoParam->GetfZDCCen())->Eval(calibSumZN1[0]+calibSumZP1[0]));
 463     nGenSpecNRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN2[0]));
 464     nGenSpecPRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZP2[0]));
 465     nGenSpecRight  = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN2[0]+calibSumZP2[0]));
 466     impPar  = (fRecoParam->GetfbCen())->Eval(calibSumZN1[0]+calibSumZP1[0]);
 467   }
 468   else if(corrADCZEMHG < infValueZEM){
 469     nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZNPer())->Eval(calibSumZN1[0]));
 470     nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZPPer())->Eval(calibSumZP1[0]));
 471     nGenSpecLeft  = (Int_t) ((fRecoParam->GetfZDCPer())->Eval(calibSumZN1[0]+calibSumZP1[0]));
 472     impPar   = (fRecoParam->GetfbPer())->Eval(calibSumZN1[0]+calibSumZP1[0]);
 473   }
 474   else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
 475     nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 476     nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 477     nGenSpecLeft  = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 478     impPar   =  (fRecoParam->GetfZEMb())->Eval(corrADCZEMHG);
 479   }
 480   //
 481   if(calibSumZN1[0]/maxValEZN1>1.)  nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 482   if(calibSumZP1[0]/maxValEZP1>1.)  nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 483   if((calibSumZN1[0]+calibSumZP1[0]/maxValEZDC1)>1.){
 484      nGenSpecLeft = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 485      impPar = (fRecoParam->GetfZEMb())->Eval(corrADCZEMHG);
 486   }
 487   if(calibSumZN2[0]/maxValEZN2>1.)  nGenSpecNRight = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 488   if(calibSumZP2[0]/maxValEZP2>1.)  nGenSpecPRight = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 489   if((calibSumZN2[0]+calibSumZP2[0]/maxValEZDC2)>1.) nGenSpecRight = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 490   //
 491   if(nGenSpecNLeft>125)    nGenSpecNLeft=125;
 492   else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
 493   if(nGenSpecPLeft>82)     nGenSpecPLeft=82;
 494   else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
 495   if(nGenSpecLeft>207)     nGenSpecLeft=207;
 496   else if(nGenSpecLeft<0)  nGenSpecLeft=0;
 497
 498   //  ---      Number of generated participants (from HIJING parameterization)
 499   Int_t nPart, nPartTotLeft, nPartTotRight;
 500   nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
 501   nPartTotLeft = 207-nGenSpecLeft;
 502   nPartTotRight = 207-nGenSpecRight;
 503   if(nPart<0) nPart=0;
 504   if(nPartTotLeft<0) nPartTotLeft=0;
 505   if(nPartTotRight<0) nPartTotRight=0;
 506   //
 507   // *** DEBUG ***
 508   /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
 509       "  calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
 510       calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
 511   printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
 512       "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n",
 513       nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft,
 514       nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
 515   printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
 516   */
 517
 518   // create the output tree
 519   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 520                   calibTowZN1, calibTowZN2, calibTowZP1, calibTowZP2,
 521                   ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
 522                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 523                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
 524                   nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
 525
 526   AliZDCReco* preco = &reco;
 527   const Int_t kBufferSize = 4000;
 528   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
 529
 530   // write the output tree
 531   clustersTree->Fill();
 532 }
 533
 534 //_____________________________________________________________________________
 535 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
 536 {
 537   // fill energies and number of participants to the ESD
 538
 539   AliZDCReco reco;
 540   AliZDCReco* preco = &reco;
 541   clustersTree->SetBranchAddress("ZDC", &preco);
 542
 543   clustersTree->GetEntry(0);
 544   //
 545   AliESDZDC * esdzdc = esd->GetESDZDC();
 546   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
 547   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
 548   for(Int_t i=0; i<5; i++){
 549      tZN1Ene[i] = reco.GetZN1HREnTow(i);
 550      tZN2Ene[i] = reco.GetZN2HREnTow(i);
 551      tZP1Ene[i] = reco.GetZP1HREnTow(i);
 552      tZP2Ene[i] = reco.GetZP2HREnTow(i);
 553      //
 554      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
 555      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
 556      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
 557      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
 558   }
 559   esdzdc->SetZN1TowerEnergy(tZN1Ene);
 560   esdzdc->SetZN2TowerEnergy(tZN2Ene);
 561   esdzdc->SetZP1TowerEnergy(tZP1Ene);
 562   esdzdc->SetZP2TowerEnergy(tZP2Ene);
 563   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
 564   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
 565   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
 566   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
 567   //
 568   esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
 569               reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
 570               reco.GetNPartLeft());
 571   //
 572
 573 }
 574
 575 //_____________________________________________________________________________
 576 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
 577 {
 578   // Setting the storage
 579
 580   Bool_t deleteManager = kFALSE;
 581
 582   AliCDBManager *manager = AliCDBManager::Instance();
 583   AliCDBStorage *defstorage = manager->GetDefaultStorage();
 584
 585   if(!defstorage || !(defstorage->Contains("ZDC"))){
 586      AliWarning("No default storage set or default storage doesn't contain ZDC!");
 587      manager->SetDefaultStorage(uri);
 588      deleteManager = kTRUE;
 589   }
 590
 591   AliCDBStorage *storage = manager->GetDefaultStorage();
 592
 593   if(deleteManager){
 594     AliCDBManager::Instance()->UnsetDefaultStorage();
 595     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
 596   }
 597
 598   return storage;
 599 }
 600
 601 //_____________________________________________________________________________
 602 AliZDCPedestals* AliZDCReconstructor::GetPedData() const
 603 {
 604
 605   // Getting pedestal calibration object for ZDC set
 606
 607   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
 608   if(!entry) AliFatal("No calibration data loaded!");
 609
 610   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
 611   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 612
 613   return calibdata;
 614 }
 615
 616 //_____________________________________________________________________________
 617 AliZDCCalib* AliZDCReconstructor::GetECalibData() const
 618 {
 619
 620   // Getting energy and equalization calibration object for ZDC set
 621
 622   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Calib");
 623   if(!entry) AliFatal("No calibration data loaded!");
 624
 625   AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*>  (entry->GetObject());
 626   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 627
 628   return calibdata;
 629 }
 630