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