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 "AliZDCRecParam.h"
  38
  39
  40 ClassImp(AliZDCReconstructor)
  41
  42
  43 //_____________________________________________________________________________
  44 AliZDCReconstructor:: AliZDCReconstructor() :
  45
  46   fZNCen(new TF1("fZNCen",
  47         "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.)),
  48   fZNPer(new TF1("fZNPer",
  49       "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.)),
  50   fZPCen(new TF1("fZPCen",
  51        "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.)),
  52   fZPPer(new TF1("fZPPer",
  53       "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.)),
  54   fZDCCen(new TF1("fZDCCen",
  55       "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.)),
  56   fZDCPer(new TF1("fZDCPer",
  57       "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.)),
  58   fbCen(new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.)),
  59   fbPer(new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.)),
  60   //
  61   fZEMn(new TF1("fZEMn","121.7-0.1934*x+0.00007565*x*x",0.,1200.)),
  62   fZEMp(new TF1("fZEMp","80.05-0.1315*x+0.00005327*x*x",0.,1200.)),
  63   fZEMsp(new TF1("fZEMsp","201.7-0.325*x+0.0001292*x*x",0.,1200.)),
  64   fZEMb(new TF1("fZEMb",
  65         "13.83-0.02851*x+5.101e-5*x*x-7.305e-8*x*x*x+5.101e-11*x*x*x*x-1.25e-14*x*x*x*x*x",0.,1200.)),
  66   //
  67   fPedData(GetPedData()),
  68   fECalibData(GetECalibData()),
  69   fRecParam(GetRecParams())
  70 {
  71   // **** Default constructor
  72
  73 }
  74
  75
  76 //_____________________________________________________________________________
  77 AliZDCReconstructor::~AliZDCReconstructor()
  78 {
  79 // destructor
  80
  81   delete fZNCen;
  82   delete fZNPer;
  83   delete fZPCen;
  84   delete fZPPer;
  85   delete fZDCCen;
  86   delete fZDCPer;
  87   delete fbCen;
  88   delete fbPer;
  89   delete fZEMn;
  90   delete fZEMp;
  91   delete fZEMsp;
  92   delete fZEMb;
  93
  94 }
  95
  96
  97 //_____________________________________________________________________________
  98 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
  99 {
 100   // *** Local ZDC reconstruction for digits
 101   // Works on the current event
 102
 103   // Retrieving calibration data
 104   Float_t meanPed[48];
 105   for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
 106
 107   // get digits
 108   AliZDCDigit digit;
 109   AliZDCDigit* pdigit = &digit;
 110   digitsTree->SetBranchAddress("ZDC", &pdigit);
 111
 112   // loop over digits
 113   Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
 114   Float_t dZEM1CorrHG=0., dZEM2CorrHG=0.;
 115   Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
 116   Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
 117   Float_t dZEM1CorrLG=0., dZEM2CorrLG=0.;
 118   Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
 119
 120   //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
 121   for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
 122    digitsTree->GetEntry(iDigit);
 123    if (!pdigit) continue;
 124    //pdigit->Print("");
 125    //
 126    Int_t det = digit.GetSector(0);
 127    Int_t quad = digit.GetSector(1);
 128    Int_t pedindex = -1;
 129    //printf("\n\t Digit #%d det %d quad %d", iDigit, det, quad);
 130    //
 131    if(quad != 5){ // ZDC (not reference PTMs!)
 132     if(det == 1){ // *** ZN1
 133        pedindex = quad;
 134        tZN1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 135        if(tZN1CorrHG[quad]<0.) tZN1CorrHG[quad] = 0.;
 136        tZN1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
 137        if(tZN1CorrLG[quad]<0.) tZN1CorrLG[quad] = 0.;
 138        //printf("\t pedindex %d tZN1CorrHG[%d] = %1.0f tZN1CorrLG[%d] = %1.0f",
 139        //       pedindex, quad, tZN1CorrHG[quad], quad, tZN1CorrLG[quad]);
 140     }
 141     else if(det == 2){ // *** ZP1
 142        pedindex = quad+5;
 143        tZP1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 144        if(tZP1CorrLG[quad]<0.) tZP1CorrLG[quad] = 0.;
 145        tZP1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
 146        if(tZP1CorrHG[quad]<0.) tZP1CorrHG[quad] = 0.;
 147        //printf("\t pedindex %d tZP1CorrHG[%d] = %1.0f tZP1CorrLG[%d] = %1.0f",
 148        //       pedindex, quad, tZP1CorrHG[quad], quad, tZP1CorrLG[quad]);
 149     }
 150     else if(det == 3){
 151        if(quad == 1){       // *** ZEM1
 152          pedindex = quad+9;
 153          dZEM1CorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 154          if(dZEM1CorrHG<0.) dZEM1CorrHG = 0.;
 155          dZEM1CorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
 156          if(dZEM1CorrLG<0.) dZEM1CorrLG = 0.;
 157          //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n",
 158          //     pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEM1CorrHG);
 159          //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n",
 160          //     pedindex+24, digit.GetADCValue(1), meanPed[pedindex+24], dZEM1CorrLG);
 161        }
 162        else if(quad == 2){  // *** ZEM2
 163          pedindex = quad+9;
 164          dZEM2CorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 165          if(dZEM2CorrHG<0.) dZEM2CorrHG = 0.;
 166          dZEM2CorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
 167          if(dZEM2CorrLG<0.) dZEM2CorrLG = 0.;
 168          //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n",
 169          //     pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEM2CorrHG);
 170          //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n",
 171          //     pedindex+2, digit.GetADCValue(1),meanPed[pedindex+2], dZEM2CorrLG);
 172        }
 173     }
 174     else if(det == 4){  // *** ZN2
 175        pedindex = quad+12;
 176        tZN2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 177        if(tZN2CorrHG[quad]<0.) tZN2CorrHG[quad] = 0.;
 178        tZN2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
 179        if(tZN2CorrLG[quad]<0.) tZN2CorrLG[quad] = 0.;
 180        //printf("\t pedindex %d tZN2CorrHG[%d] = %1.0f tZN2CorrLG[%d] = %1.0f\n",
 181        //       pedindex, quad, tZN2CorrHG[quad], quad, tZN2CorrLG[quad]);
 182     }
 183     else if(det == 5){  // *** ZP2
 184        pedindex = quad+17;
 185        tZP2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 186        if(tZP2CorrHG[quad]<0.) tZP2CorrHG[quad] = 0.;
 187        tZP2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+24]);
 188        if(tZP2CorrLG[quad]<0.) tZP2CorrLG[quad] = 0.;
 189        //printf("\t pedindex %d tZP2CorrHG[%d] = %1.0f tZP2CorrLG[%d] = %1.0f\n",
 190        //       pedindex, quad, tZP2CorrHG[quad], quad, tZP2CorrLG[quad]);
 191     }
 192    }
 193   }
 194
 195   // reconstruct the event
 196     ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG,
 197         tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG,
 198         tZP2CorrLG, dZEM1CorrHG, dZEM2CorrHG);
 199
 200 }
 201
 202 //_____________________________________________________________________________
 203 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
 204 {
 205   // *** ZDC raw data reconstruction
 206   // Works on the current event
 207
 208   // Retrieving calibration data
 209   Float_t meanPed[48];
 210   for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
 211
 212   rawReader->Reset();
 213
 214   // loop over raw data rawDatas
 215   Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
 216   Float_t dZEM1CorrHG=0., dZEM2CorrHG=0.;
 217   Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
 218   Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
 219   Float_t dZEM1CorrLG=0., dZEM2CorrLG=0.;
 220   Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
 221   //
 222   AliZDCRawStream rawData(rawReader);
 223   while (rawData.Next()) {
 224     if(rawData.IsADCDataWord()){
 225      Int_t det = rawData.GetSector(0);
 226      Int_t quad = rawData.GetSector(1);
 227      Int_t gain = rawData.GetADCGain();
 228      Int_t pedindex=0;
 229      //
 230      if(quad !=5){ // ZDCs (not reference PTMs)
 231       if(det == 1){
 232         pedindex = quad;
 233         if(gain == 0) tZN1CorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 234         else tZN1CorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
 235       }
 236       else if(det == 2){
 237         pedindex = quad+5;
 238         if(gain == 0) tZP1CorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 239         else tZP1CorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
 240       }
 241       else if(det == 3){
 242         if(quad==1){
 243           pedindex = quad+9;
 244           if(gain == 0) dZEM1CorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 245           else dZEM1CorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
 246         }
 247         else if(quad==2){
 248           pedindex = quad+9;
 249           if(gain == 0) dZEM2CorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 250           else dZEM2CorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
 251         }
 252       }
 253       else if(det == 4){
 254         pedindex = quad+12;
 255         if(gain == 0) tZN2CorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 256         else tZN2CorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
 257       }
 258       else if(det == 5){
 259         pedindex = quad+17;
 260         if(gain == 0) tZP2CorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 261         else tZP2CorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+24]);
 262       }
 263       printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
 264         det,quad,gain, pedindex, meanPed[pedindex]);
 265      }
 266     }//IsADCDataWord
 267   }
 268
 269   // reconstruct the event
 270     ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG,
 271         tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG,
 272         tZP2CorrLG, dZEM1CorrHG, dZEM2CorrHG);
 273
 274 }
 275
 276 //_____________________________________________________________________________
 277 void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree,
 278                 Float_t* ZN1ADCCorrHG, Float_t* ZP1ADCCorrHG,
 279                 Float_t* ZN2ADCCorrHG, Float_t* ZP2ADCCorrHG,
 280                 Float_t* ZN1ADCCorrLG, Float_t* ZP1ADCCorrLG,
 281                 Float_t* ZN2ADCCorrLG, Float_t* ZP2ADCCorrLG,
 282                 Float_t corrADCZEM1HG, Float_t corrADCZEM2HG) const
 283 {
 284   // ***** Reconstruct one event
 285
 286   // *** RECONSTRUCTION FROM SIMULATED DATA
 287   // It passes trhough the no. of phe which is known from simulations
 288   //  ---      ADCchannel -> photoelectrons
 289   // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
 290   // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
 291   //Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
 292   //zn1phe  = ZN1Corr/convFactor;
 293   //zp1phe  = ZP1Corr/convFactor;
 294   //zemphe = ZEMCorr/convFactor;
 295   //zn2phe  = ZN2Corr/convFactor;
 296   //zp2phe  = ZP2Corr/convFactor;
 297   ////if AliDebug(1,Form("\n    znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
 298   //
 299   ////  ---      Energy calibration
 300   //// Conversion factors for hadronic ZDCs goes from phe yield to TRUE
 301   //// incident energy (conversion from GeV to TeV is included); while for EM
 302   //// calos conversion is from light yield to detected energy calculated by
 303   //// GEANT NB -> ZN and ZP conversion factors are constant since incident
 304   //// spectators have all the same energy, ZEM energy is obtained through a
 305   //// fit over the whole range of incident particle energies
 306   //// (obtained with full HIJING simulations)
 307   //Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
 308   //Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
 309   //zn1energy  = zn1phe/zn1phexTeV;
 310   //zp1energy  = zp1phe/zp1phexTeV;
 311   //zdc1energy = zn1energy+zp1energy;
 312   //zn2energy  = zn2phe/zn2phexTeV;
 313   //zp2energy  = zp2phe/zp2phexTeV;
 314   //zdc2energy = zn2energy+zp2energy;
 315   //zemenergy = -4.81+0.3238*zemphe;
 316   //if(zemenergy<0) zemenergy=0;
 317   ////  if AliDebug(1,Form("    znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
 318   ////                     "\n          zemenergy = %f TeV\n", znenergy, zpenergy,
 319   ////                     zdcenergy, zemenergy);
 320   ////  if(zdcenergy==0)
 321   ////    if AliDebug(1,Form("\n\n      ###     ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
 322   ////                       " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy);
 323
 324   //
 325   // *** RECONSTRUCTION FROM "REAL" DATA
 326   //
 327   // Retrieving calibration data
 328   // --- Equalization coefficients ---------------------------------------------
 329   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 330   for(Int_t ji=0; ji<5; ji++){
 331      equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
 332      equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
 333      equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
 334      equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
 335   }
 336   // --- Energy calibration factors ------------------------------------
 337   Float_t calibEne[4];
 338   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
 339   //
 340   // --- Reconstruction parameters ------------------
 341   Float_t endPointZEM = fRecParam->GetZEMEndValue();
 342   Float_t cutFractionZEM = fRecParam->GetZEMCutFraction();
 343   Float_t dZEMSup = fRecParam->GetDZEMSup();
 344   Float_t dZEMInf = fRecParam->GetDZEMInf();
 345   //
 346   Float_t cutValueZEM = endPointZEM*cutFractionZEM;
 347   Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
 348   Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
 349   //
 350   Float_t maxValEZN1 = fRecParam->GetEZN1MaxValue();
 351   Float_t maxValEZP1 = fRecParam->GetEZP1MaxValue();
 352   Float_t maxValEZDC1 = fRecParam->GetEZDC1MaxValue();
 353   Float_t maxValEZN2 = fRecParam->GetEZN2MaxValue();
 354   Float_t maxValEZP2 = fRecParam->GetEZP2MaxValue();
 355   Float_t maxValEZDC2 = fRecParam->GetEZDC2MaxValue();
 356   //
 357   //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
 358   //   " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
 359
 360   // Equalization of detector responses
 361   Float_t equalTowZN1HG[5], equalTowZN2HG[5], equalTowZP1HG[5], equalTowZP2HG[5];
 362   Float_t equalTowZN1LG[5], equalTowZN2LG[5], equalTowZP1LG[5], equalTowZP2LG[5];
 363   for(Int_t gi=0; gi<5; gi++){
 364      equalTowZN1HG[gi] = ZN1ADCCorrHG[gi]*equalCoeffZN1[gi];
 365      equalTowZP1HG[gi] = ZP1ADCCorrHG[gi]*equalCoeffZP1[gi];
 366      equalTowZN2HG[gi] = ZN2ADCCorrHG[gi]*equalCoeffZN2[gi];
 367      equalTowZP2HG[gi] = ZP2ADCCorrHG[gi]*equalCoeffZP2[gi];
 368      //
 369      equalTowZN1LG[gi] = ZN1ADCCorrLG[gi]*equalCoeffZN1[gi];
 370      equalTowZP1LG[gi] = ZP1ADCCorrLG[gi]*equalCoeffZP1[gi];
 371      equalTowZN2LG[gi] = ZN2ADCCorrLG[gi]*equalCoeffZN2[gi];
 372      equalTowZP2LG[gi] = ZP2ADCCorrLG[gi]*equalCoeffZP2[gi];
 373   }
 374
 375   // Energy calibration of detector responses
 376   Float_t calibTowZN1HG[5], calibTowZN2HG[5], calibTowZP1HG[5], calibTowZP2HG[5];
 377   Float_t calibSumZN1HG=0., calibSumZN2HG=0., calibSumZP1HG=0., calibSumZP2HG=0.;
 378   Float_t calibTowZN1LG[5], calibTowZN2LG[5], calibTowZP1LG[5], calibTowZP2LG[5];
 379   Float_t calibSumZN1LG=0., calibSumZN2LG=0., calibSumZ12LG=0., calibSumZP2LG=0.;
 380   for(Int_t gi=0; gi<5; gi++){
 381      calibTowZN1HG[gi] = equalTowZN1HG[gi]*calibEne[0];
 382      calibTowZP1HG[gi] = equalTowZP1HG[gi]*calibEne[1];
 383      calibTowZN2HG[gi] = equalTowZN2HG[gi]*calibEne[2];
 384      calibTowZP2HG[gi] = equalTowZP2HG[gi]*calibEne[3];
 385      calibSumZN1HG += calibTowZN1HG[gi];
 386      calibSumZP1HG += calibTowZP1HG[gi];
 387      calibSumZN2HG += calibTowZN2HG[gi];
 388      calibSumZP2HG += calibTowZP2HG[gi];
 389      //
 390      calibTowZN1LG[gi] = equalTowZN1LG[gi]*calibEne[0];
 391      calibTowZP1LG[gi] = equalTowZP1LG[gi]*calibEne[1];
 392      calibTowZN2LG[gi] = equalTowZN2LG[gi]*calibEne[2];
 393      calibTowZP2LG[gi] = equalTowZP2LG[gi]*calibEne[3];
 394      calibSumZN1LG += calibTowZN1LG[gi];
 395      calibSumZ12LG += calibTowZP1LG[gi];
 396      calibSumZN2LG += calibTowZN2LG[gi];
 397      calibSumZP2LG += calibTowZP2LG[gi];
 398   }
 399
 400   //  ---      Number of detected spectator nucleons
 401   //  *** N.B. -> It works only in Pb-Pb
 402   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
 403   nDetSpecNLeft = (Int_t) (calibSumZN1HG/2.760);
 404   nDetSpecPLeft = (Int_t) (calibSumZP1HG/2.760);
 405   nDetSpecNRight = (Int_t) (calibSumZN2HG/2.760);
 406   nDetSpecPRight = (Int_t) (calibSumZP2HG/2.760);
 407   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
 408     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
 409     nDetSpecNRight, nDetSpecPRight);*/
 410
 411   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
 412   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
 413   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
 414   Double_t impPar=0.;
 415   //
 416   // *** RECONSTRUCTION FROM SIMULATED DATA
 417   // Cut value for Ezem (GeV)
 418   // ### Results from production  -> 0<b<18 fm (Apr 2002)
 419   /*Float_t eZEMCut = 420.;
 420   Float_t deltaEZEMSup = 690.;
 421   Float_t deltaEZEMInf = 270.;
 422   if(zemenergy > (eZEMCut+deltaEZEMSup)){
 423     nGenSpecNLeft  = (Int_t) (fZNCen->Eval(ZN1CalibSum));
 424     nGenSpecPLeft  = (Int_t) (fZPCen->Eval(ZP1CalibSum));
 425     nGenSpecLeft   = (Int_t) (fZDCCen->Eval(ZN1CalibSum+ZP1CalibSum));
 426     nGenSpecNRight = (Int_t) (fZNCen->Eval(ZN2CalibSum));
 427     nGenSpecPRight = (Int_t) (fZNCen->Eval(ZP2CalibSum));
 428     nGenSpecRight  = (Int_t) (fZNCen->Eval(ZN2CalibSum+ZP2CalibSum));
 429     impPar  = fbCen->Eval(ZN1CalibSum+ZP1CalibSum);
 430   }
 431   else if(zemenergy < (eZEMCut-deltaEZEMInf)){
 432     nGenSpecNLeft = (Int_t) (fZNPer->Eval(ZN1CalibSum));
 433     nGenSpecPLeft = (Int_t) (fZPPer->Eval(ZP1CalibSum));
 434     nGenSpecLeft  = (Int_t) (fZDCPer->Eval(ZN1CalibSum+ZP1CalibSum));
 435     impPar   = fbPer->Eval(ZN1CalibSum+ZP1CalibSum);
 436   }
 437   else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
 438     nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
 439     nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
 440     nGenSpecLeft  = (Int_t)(fZEMsp->Eval(zemenergy));
 441     impPar   =  fZEMb->Eval(zemenergy);
 442   }
 443   // ### Results from production  -> 0<b<18 fm (Apr 2002)
 444   if(ZN1CalibSum>162.)  nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
 445   if(ZP1CalibSum>59.75)  nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
 446   if(ZN1CalibSum+ZP1CalibSum>221.5) nGenSpecLeft  = (Int_t)(fZEMsp->Eval(zemenergy));
 447   if(ZN1CalibSum+ZP1CalibSum>220.)  impPar    =  fZEMb->Eval(zemenergy);
 448   */
 449   //
 450   //
 451   // *** RECONSTRUCTION FROM REAL DATA
 452   //
 453   Float_t corrADCZEMHG = corrADCZEM1HG + corrADCZEM2HG;
 454   //
 455   if(corrADCZEMHG > supValueZEM){
 456     nGenSpecNLeft  = (Int_t) (fZNCen->Eval(calibSumZN1HG));
 457     nGenSpecPLeft  = (Int_t) (fZPCen->Eval(calibSumZP1HG));
 458     nGenSpecLeft   = (Int_t) (fZDCCen->Eval(calibSumZN1HG+calibSumZP1HG));
 459     nGenSpecNRight = (Int_t) (fZNCen->Eval(calibSumZN2HG));
 460     nGenSpecPRight = (Int_t) (fZNCen->Eval(calibSumZP2HG));
 461     nGenSpecRight  = (Int_t) (fZNCen->Eval(calibSumZN2HG+calibSumZP2HG));
 462     impPar  = fbCen->Eval(calibSumZN1HG+calibSumZP1HG);
 463   }
 464   else if(corrADCZEMHG < infValueZEM){
 465     nGenSpecNLeft = (Int_t) (fZNPer->Eval(calibSumZN1HG));
 466     nGenSpecPLeft = (Int_t) (fZPPer->Eval(calibSumZP1HG));
 467     nGenSpecLeft  = (Int_t) (fZDCPer->Eval(calibSumZN1HG+calibSumZP1HG));
 468     impPar   = fbPer->Eval(calibSumZN1HG+calibSumZP1HG);
 469   }
 470   else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
 471     nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
 472     nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
 473     nGenSpecLeft  = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
 474     impPar   =  fZEMb->Eval(corrADCZEMHG);
 475   }
 476   //
 477   if(calibSumZN1HG/maxValEZN1>1.)  nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
 478   if(calibSumZP1HG/maxValEZP1>1.)  nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
 479   if((calibSumZN1HG+calibSumZP1HG/maxValEZDC1)>1.){
 480      nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
 481      impPar = fZEMb->Eval(corrADCZEMHG);
 482   }
 483   if(calibSumZN2HG/maxValEZN2>1.)  nGenSpecNRight = (Int_t) (fZEMn->Eval(corrADCZEMHG));
 484   if(calibSumZP2HG/maxValEZP2>1.)  nGenSpecPRight = (Int_t) (fZEMp->Eval(corrADCZEMHG));
 485   if((calibSumZN2HG+calibSumZP2HG/maxValEZDC2)>1.) nGenSpecRight = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
 486   //
 487   if(nGenSpecNLeft>125)    nGenSpecNLeft=125;
 488   else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
 489   if(nGenSpecPLeft>82)     nGenSpecPLeft=82;
 490   else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
 491   if(nGenSpecLeft>207)     nGenSpecLeft=207;
 492   else if(nGenSpecLeft<0)  nGenSpecLeft=0;
 493
 494   //  ---      Number of generated participants (from HIJING parameterization)
 495   Int_t nPart, nPartTotLeft, nPartTotRight;
 496   nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
 497   nPartTotLeft = 207-nGenSpecLeft;
 498   nPartTotRight = 207-nGenSpecRight;
 499   if(nPart<0) nPart=0;
 500   if(nPartTotLeft<0) nPartTotLeft=0;
 501   if(nPartTotRight<0) nPartTotRight=0;
 502   //
 503   // *** DEBUG ***
 504   printf("\n\t AliZDCReconstructor -> calibSumZN1HG %1.0f, calibSumZP1HG %1.0f,"
 505       "  calibSumZN2HG %1.0f, calibSumZP2HG %1.0f, corrADCZEMHG %1.0f\n",
 506       calibSumZN1HG,calibSumZP1HG,calibSumZN2HG,calibSumZP2HG,corrADCZEMHG);
 507   printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
 508       "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n",
 509       nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft,
 510       nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
 511   printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
 512
 513   // create the output tree
 514   AliZDCReco reco(calibSumZN1HG, calibSumZP1HG, calibSumZN2HG, calibSumZP2HG,
 515                   calibTowZN1LG, calibTowZN2LG, calibTowZP1LG, calibTowZP2LG,
 516                   calibTowZN1LG, calibTowZP1LG, calibTowZN2LG, calibTowZP2LG,
 517                   corrADCZEM1HG, corrADCZEM2HG,
 518                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 519                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
 520                   nGenSpecPRight, nGenSpecRight,
 521                   nPartTotLeft, nPartTotRight, impPar);
 522
 523   AliZDCReco* preco = &reco;
 524   const Int_t kBufferSize = 4000;
 525   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
 526
 527   // write the output tree
 528   clustersTree->Fill();
 529 }
 530
 531 //_____________________________________________________________________________
 532 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
 533 {
 534   // fill energies and number of participants to the ESD
 535
 536   AliZDCReco reco;
 537   AliZDCReco* preco = &reco;
 538   clustersTree->SetBranchAddress("ZDC", &preco);
 539
 540   clustersTree->GetEntry(0);
 541   //
 542   AliESDZDC * esdzdc = esd->GetESDZDC();
 543   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
 544   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
 545   for(Int_t i=0; i<5; i++){
 546      tZN1Ene[i] = reco.GetZN1EnTow(i);
 547      tZN2Ene[i] = reco.GetZN2EnTow(i);
 548      tZP1Ene[i] = reco.GetZP1EnTow(i);
 549      tZP2Ene[i] = reco.GetZP2EnTow(i);
 550      tZN1EneLR[i] = reco.GetZN1SigLowRes(i);
 551      tZN2EneLR[i] = reco.GetZP1SigLowRes(i);
 552      tZP1EneLR[i] = reco.GetZN2SigLowRes(i);
 553      tZP2EneLR[i] = reco.GetZP2SigLowRes(i);
 554   }
 555   esdzdc->SetZN1TowerEnergy(tZN1Ene);
 556   esdzdc->SetZN2TowerEnergy(tZN2Ene);
 557   esdzdc->SetZP1TowerEnergy(tZP1Ene);
 558   esdzdc->SetZP2TowerEnergy(tZP2Ene);
 559   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
 560   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
 561   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
 562   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
 563   //
 564   esd->SetZDC(reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEM1signal(),
 565               reco.GetZEM2signal(), reco.GetZN2Energy(), reco.GetZP2Energy(),
 566               reco.GetNPartLeft());
 567   //
 568
 569 }
 570
 571 //_____________________________________________________________________________
 572 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
 573 {
 574   // Setting the storage
 575
 576   Bool_t deleteManager = kFALSE;
 577
 578   AliCDBManager *manager = AliCDBManager::Instance();
 579   AliCDBStorage *defstorage = manager->GetDefaultStorage();
 580
 581   if(!defstorage || !(defstorage->Contains("ZDC"))){
 582      AliWarning("No default storage set or default storage doesn't contain ZDC!");
 583      manager->SetDefaultStorage(uri);
 584      deleteManager = kTRUE;
 585   }
 586
 587   AliCDBStorage *storage = manager->GetDefaultStorage();
 588
 589   if(deleteManager){
 590     AliCDBManager::Instance()->UnsetDefaultStorage();
 591     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
 592   }
 593
 594   return storage;
 595 }
 596
 597 //_____________________________________________________________________________
 598 AliZDCPedestals* AliZDCReconstructor::GetPedData() const
 599 {
 600
 601   // Getting pedestal calibration object for ZDC set
 602
 603   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
 604   if(!entry) AliFatal("No calibration data loaded!");
 605
 606   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
 607   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 608
 609   return calibdata;
 610 }
 611
 612 //_____________________________________________________________________________
 613 AliZDCCalib* AliZDCReconstructor::GetECalibData() const
 614 {
 615
 616   // Getting energy and equalization calibration object for ZDC set
 617
 618   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Calib");
 619   if(!entry) AliFatal("No calibration data loaded!");
 620
 621   AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*>  (entry->GetObject());
 622   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 623
 624   return calibdata;
 625 }
 626
 627 //_____________________________________________________________________________
 628 AliZDCRecParam* AliZDCReconstructor::GetRecParams() const
 629 {
 630
 631   // Getting energy and equalization calibration object for ZDC set
 632
 633   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecParam");
 634   if(!entry) AliFatal("No calibration data loaded!");
 635
 636   AliZDCRecParam *calibdata = dynamic_cast<AliZDCRecParam*>  (entry->GetObject());
 637   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 638
 639   return calibdata;
 640 }