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