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