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 #include <TMap.h>
  27
  28 #include "AliRunLoader.h"
  29 #include "AliRawReader.h"
  30 #include "AliGRPObject.h"
  31 #include "AliESDEvent.h"
  32 #include "AliESDZDC.h"
  33 #include "AliZDCDigit.h"
  34 #include "AliZDCRawStream.h"
  35 #include "AliZDCReco.h"
  36 #include "AliZDCReconstructor.h"
  37 #include "AliZDCPedestals.h"
  38 #include "AliZDCCalib.h"
  39 #include "AliZDCRecoParam.h"
  40 #include "AliZDCRecoParampp.h"
  41 #include "AliZDCRecoParamPbPb.h"
  42
  43
  44 ClassImp(AliZDCReconstructor)
  45 AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0;  //reconstruction parameters
  46
  47 //_____________________________________________________________________________
  48 AliZDCReconstructor:: AliZDCReconstructor() :
  49   fPedData(GetPedData()),
  50   fECalibData(GetECalibData()),
  51   fRecoMode(0),
  52   fBeamEnergy(0.)
  53 {
  54   // **** Default constructor
  55   SetRecoMode();
  56
  57 }
  58
  59
  60 //_____________________________________________________________________________
  61 AliZDCReconstructor::~AliZDCReconstructor()
  62 {
  63 // destructor
  64    if(fRecoParam)  delete fRecoParam;
  65    if(fPedData)    delete fPedData;
  66    if(fECalibData) delete fECalibData;
  67 }
  68
  69 //____________________________________________________________________________
  70 void AliZDCReconstructor::SetRecoMode()
  71 {
  72   // Setting reconstruction mode
  73
  74   // Initialization of the GRP entry
  75   AliGRPObject* grpData;
  76   AliCDBEntry*  entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
  77   if(entry){
  78     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
  79     if(m){
  80        m->Print();
  81        grpData = new AliGRPObject();
  82        grpData->ReadValuesFromMap(m);
  83     }
  84     else{
  85        grpData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
  86        entry->SetOwner(0);
  87     }
  88     AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
  89   }
  90   if(!grpData) AliError("No GRP entry found in OCDB!");
  91
  92   TString beamType = grpData->GetBeamType();
  93   if(beamType==AliGRPObject::GetInvalidString()){
  94     AliError("GRP/GRP/Data entry:  missing value for the beam energy !");
  95     AliError("\t ZDC does not reconstruct event 4 UNKNOWN beam type\n");
  96     return;
  97   }
  98   //
  99   if((beamType.CompareTo("p-p")) == 0){
 100     fRecoMode=0;
 101     fRecoParam = (AliZDCRecoParampp*) AliZDCRecoParampp::GetppRecoParam();
 102   }
 103   else if((beamType.CompareTo("A-A")) == 0){
 104     fRecoMode=1;
 105     fRecoParam = (AliZDCRecoParamPbPb*) AliZDCRecoParamPbPb::GetPbPbRecoParam();
 106   }
 107
 108   fBeamEnergy = grpData->GetBeamEnergy();
 109   if(fBeamEnergy==AliGRPObject::GetInvalidFloat()) {
 110     AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
 111     fBeamEnergy = 0.;
 112   }
 113
 114   printf("\n ***** ZDC reconstruction initialized for %s @ %1.3f GeV\n\n",beamType.Data(), fBeamEnergy);
 115 }
 116
 117 //_____________________________________________________________________________
 118 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
 119 {
 120   // *** Local ZDC reconstruction for digits
 121   // Works on the current event
 122
 123   // Retrieving calibration data
 124   Float_t meanPed[48];
 125   for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
 126
 127   // get digits
 128   AliZDCDigit digit;
 129   AliZDCDigit* pdigit = &digit;
 130   digitsTree->SetBranchAddress("ZDC", &pdigit);
 131   //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
 132
 133   // loop over digits
 134   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
 135   Float_t dZEM1Corr[2], dZEM2Corr[2], PMRef1[2], PMRef2[2];
 136   for(Int_t i=0; i<10; i++){
 137      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
 138      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = PMRef1[i] = PMRef2[i] = 0.;
 139   }
 140   //
 141   for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
 142    digitsTree->GetEntry(iDigit);
 143    if (!pdigit) continue;
 144    //
 145    Int_t det = digit.GetSector(0);
 146    Int_t quad = digit.GetSector(1);
 147    Int_t pedindex = -1, kNch = 24;
 148    //printf("\n\t Digit #%d det %d quad %d", iDigit, det, quad);
 149    //
 150    if(quad != 5){ // ZDC (not reference PTMs!)
 151     if(det == 1){ // *** ZNC
 152        pedindex = quad;
 153        tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 154        if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
 155        tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 156        if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
 157        //printf("\t pedindex %d tZN1Corr[%d] = %1.0f tZN1Corr[%d] = %1.0f",
 158        //       pedindex, quad, tZN1Corr[quad], quad+5, tZN1Corr[quad+5]);
 159     }
 160     else if(det == 2){ // *** ZP1
 161        pedindex = quad+5;
 162        tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 163        if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
 164        tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 165        if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
 166        //printf("\t pedindex %d tZP1Corr[%d] = %1.0f tZP1Corr[%d] = %1.0f",
 167        //       pedindex, quad, tZP1Corr[quad], quad+5, tZP1Corr[quad+5]);
 168     }
 169     else if(det == 3){
 170        pedindex = quad+9;
 171        if(quad == 1){       // *** ZEM1
 172          dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 173          if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
 174          dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 175          if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
 176          //printf("\t pedindex %d tZEM1Corr[%d] = %1.0f tZEM1Corr[%d] = %1.0f",
 177          //     pedindex, quad, tZEM1Corr[quad], quad+1, tZEM1Corr[quad+1]);
 178        }
 179        else if(quad == 2){  // *** ZEM2
 180          dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 181          if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
 182          dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 183          if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
 184          //printf("\t pedindex %d tZEM2Corr[%d] = %1.0f tZEM2Corr[%d] = %1.0f",
 185          //     pedindex, quad, tZEM2Corr[quad], quad+1, tZEM2Corr[quad+1]);
 186        }
 187     }
 188     else if(det == 4){  // *** ZN2
 189        pedindex = quad+12;
 190        tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 191        if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
 192        tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 193        if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
 194        //printf("\t pedindex %d tZN2Corr[%d] = %1.0f tZN2Corr[%d] = %1.0f\n",
 195        //       pedindex, quad, tZN2Corr[quad], quad+5, tZN2Corr[quad+5]);
 196     }
 197     else if(det == 5){  // *** ZP2
 198        pedindex = quad+17;
 199        tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 200        if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
 201        tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 202        if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
 203        //printf("\t pedindex %d tZP2Corr[%d] = %1.0f tZP2Corr[%d] = %1.0f\n",
 204        //       pedindex, quad, tZP2Corr[quad], quad+5, tZP2Corr[quad+5]);
 205     }
 206    }
 207    else{ // Reference PMs
 208      pedindex = (det-1)/3+22;
 209      if(det == 1){
 210        PMRef1[0] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 211        if(PMRef1[0]<0.) PMRef1[0] = 0.;
 212        PMRef1[1] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 213        if(PMRef2[1]<0.) PMRef1[1] = 0.;
 214      }
 215      else if(det == 4){
 216        PMRef2[0] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
 217        if(PMRef2[0]<0.) PMRef2[0] = 0.;
 218        PMRef2[1] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
 219        if(PMRef2[1]<0.) PMRef2[1] = 0.;
 220      }
 221    }
 222   }
 223
 224   // reconstruct the event
 225   if(fRecoMode==0)
 226     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 227         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
 228   else if(fRecoMode==1)
 229     ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 230         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
 231
 232 }
 233
 234 //_____________________________________________________________________________
 235 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
 236 {
 237   // *** ZDC raw data reconstruction
 238   // Works on the current event
 239
 240   // Retrieving calibration data
 241   Float_t meanPed[48];
 242   for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
 243
 244   rawReader->Reset();
 245
 246   // loop over raw data
 247   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
 248   Float_t dZEM1Corr[2], dZEM2Corr[2], PMRef1[2], PMRef2[2];
 249   for(Int_t i=0; i<10; i++){
 250      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
 251      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = PMRef1[i] = PMRef2[i] = 0.;
 252   }
 253   //
 254   AliZDCRawStream rawData(rawReader);
 255   Int_t kNch = 24;
 256   while (rawData.Next()) {
 257     if(rawData.IsADCDataWord()){
 258      Int_t det = rawData.GetSector(0);
 259      Int_t quad = rawData.GetSector(1);
 260      Int_t gain = rawData.GetADCGain();
 261      Int_t pedindex=0;
 262      //
 263      if(quad !=5){ // ZDCs (not reference PTMs)
 264       if(det == 1){
 265         pedindex = quad;
 266         if(gain == 0) tZN1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 267         else tZN1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 268       }
 269       else if(det == 2){
 270         pedindex = quad+5;
 271         if(gain == 0) tZP1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 272         else tZP1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 273       }
 274       else if(det == 3){
 275         pedindex = quad+9;
 276         if(quad==1){
 277           if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 278           else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 279         }
 280         else if(quad==2){
 281           if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 282           else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 283         }
 284       }
 285       else if(det == 4){
 286         pedindex = quad+12;
 287         if(gain == 0) tZN2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 288         else tZN2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 289       }
 290       else if(det == 5){
 291         pedindex = quad+17;
 292         if(gain == 0) tZP2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 293         else tZP2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
 294       }
 295       //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
 296       //  det,quad,gain, pedindex, meanPed[pedindex]);
 297      }
 298      else{ // reference PM
 299        pedindex = (det-1)/3 + 22;
 300        if(det == 1){
 301          if(gain==0) PMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 302          else PMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 303        }
 304        else if(det ==4){
 305          if(gain==0) PMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 306          else PMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
 307        }
 308      }
 309     }//IsADCDataWord
 310   }
 311
 312   // reconstruct the event
 313   if(fRecoMode==0)
 314     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 315         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
 316   else if(fRecoMode==1)
 317     ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
 318         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
 319
 320 }
 321
 322 //_____________________________________________________________________________
 323 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* ZN1ADCCorr,
 324         Float_t* ZP1ADCCorr, Float_t* ZN2ADCCorr, Float_t* ZP2ADCCorr,
 325         Float_t* ZEM1ADCCorr, Float_t* ZEM2ADCCorr, Float_t* PMRef1, Float_t* PMRef2) const
 326 {
 327   // ***** Reconstruct one event
 328
 329   // *** RECONSTRUCTION FROM "REAL" DATA
 330   //
 331   // Retrieving calibration data
 332   // --- Equalization coefficients ---------------------------------------------
 333   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 334   for(Int_t ji=0; ji<5; ji++){
 335      equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
 336      equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
 337      equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
 338      equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
 339   }
 340   // --- Energy calibration factors ------------------------------------
 341   Float_t calibEne[4];
 342   // *********************************************************************
 343   // **** Until the beam type info isn't known @ reconstruction level ****
 344   // **** the energy calibration coefficient are manually set to 1    ****
 345   // **** as it will be in real life for pp data taking               ****
 346   // *********************************************************************
 347   //for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
 348   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = 1.;
 349
 350   // Equalization of detector responses
 351   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 352   for(Int_t gi=0; gi<5; gi++){
 353      equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
 354      equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
 355      equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
 356      equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
 357      equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
 358      equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
 359      equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
 360      equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
 361   }
 362
 363   // Energy calibration of detector responses
 364   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 365   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 366   for(Int_t gi=0; gi<10; gi++){
 367      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
 368      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
 369      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
 370      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
 371      //
 372      if(gi<5){
 373        calibSumZN1[0] += calibTowZN1[gi];
 374        calibSumZP1[0] += calibTowZP1[gi];
 375        calibSumZN2[0] += calibTowZN2[gi];
 376        calibSumZP2[0] += calibTowZP2[gi];
 377      }
 378      //
 379      else{
 380        calibSumZN1[1] += calibTowZN1[gi];
 381        calibSumZP1[1] += calibTowZP1[gi];
 382        calibSumZN2[1] += calibTowZN2[gi];
 383        calibSumZP2[1] += calibTowZP2[gi];
 384      }
 385   }
 386
 387   //  ---      Number of detected spectator nucleons
 388   //  *** N.B. -> It works only in Pb-Pb
 389   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
 390   if(fBeamEnergy!=0){
 391    nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
 392    nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
 393    nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
 394    nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
 395   }
 396   else AliWarning(" ATTENTION -> fBeamEnergy = 0\n");
 397   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
 398     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
 399     nDetSpecNRight, nDetSpecPRight);*/
 400
 401   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
 402   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
 403   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
 404   Int_t nPartTotLeft=0, nPartTotRight=0;
 405   Double_t impPar=0.;
 406
 407   // create the output tree
 408   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 409                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
 410                   ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
 411                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 412                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
 413                   nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
 414
 415   AliZDCReco* preco = &reco;
 416   const Int_t kBufferSize = 4000;
 417   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
 418
 419   // write the output tree
 420   clustersTree->Fill();
 421 }
 422
 423 //_____________________________________________________________________________
 424 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree, Float_t* ZN1ADCCorr,
 425         Float_t* ZP1ADCCorr, Float_t* ZN2ADCCorr, Float_t* ZP2ADCCorr,
 426         Float_t* ZEM1ADCCorr, Float_t* ZEM2ADCCorr, Float_t* PMRef1, Float_t* PMRef2) const
 427 {
 428   // ***** Reconstruct one event
 429
 430   // *** RECONSTRUCTION FROM "REAL" DATA
 431   //
 432   // Retrieving calibration data
 433   // --- Equalization coefficients ---------------------------------------------
 434   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
 435   for(Int_t ji=0; ji<5; ji++){
 436      equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
 437      equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji);
 438      equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji);
 439      equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji);
 440   }
 441   // --- Energy calibration factors ------------------------------------
 442   Float_t calibEne[4];
 443   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
 444
 445   // Equalization of detector responses
 446   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
 447   for(Int_t gi=0; gi<5; gi++){
 448      equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
 449      equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
 450      equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
 451      equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
 452      equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
 453      equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
 454      equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
 455      equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
 456   }
 457
 458   // Energy calibration of detector responses
 459   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
 460   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
 461   for(Int_t gi=0; gi<10; gi++){
 462      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
 463      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
 464      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
 465      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
 466      //
 467      if(gi<5){
 468        calibSumZN1[0] += calibTowZN1[gi];
 469        calibSumZP1[0] += calibTowZP1[gi];
 470        calibSumZN2[0] += calibTowZN2[gi];
 471        calibSumZP2[0] += calibTowZP2[gi];
 472      }
 473      //
 474      else{
 475        calibSumZN1[1] += calibTowZN1[gi];
 476        calibSumZP1[1] += calibTowZP1[gi];
 477        calibSumZN2[1] += calibTowZN2[gi];
 478        calibSumZP2[1] += calibTowZP2[gi];
 479      }
 480   }
 481
 482   //
 483   // --- Reconstruction parameters ------------------
 484   if(!fRecoParam)  fRecoParam = (AliZDCRecoParamPbPb*) AliZDCRecoParamPbPb::GetPbPbRecoParam();
 485   //
 486   Float_t endPointZEM = fRecoParam->GetZEMEndValue();
 487   Float_t cutFractionZEM = fRecoParam->GetZEMCutFraction();
 488   Float_t dZEMSup = fRecoParam->GetDZEMSup();
 489   Float_t dZEMInf = fRecoParam->GetDZEMInf();
 490   //
 491   Float_t cutValueZEM = endPointZEM*cutFractionZEM;
 492   Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
 493   Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
 494   //
 495   Float_t maxValEZN1  = fRecoParam->GetEZN1MaxValue();
 496   Float_t maxValEZP1  = fRecoParam->GetEZP1MaxValue();
 497   Float_t maxValEZDC1 = fRecoParam->GetEZDC1MaxValue();
 498   Float_t maxValEZN2  = fRecoParam->GetEZN2MaxValue();
 499   Float_t maxValEZP2  = fRecoParam->GetEZP2MaxValue();
 500   Float_t maxValEZDC2 = fRecoParam->GetEZDC2MaxValue();
 501   //
 502   //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
 503   //   " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
 504
 505   //  ---      Number of detected spectator nucleons
 506   //  *** N.B. -> It works only in Pb-Pb
 507   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
 508   if(fBeamEnergy!=0){
 509     nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
 510     nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
 511     nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
 512     nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
 513   }
 514   else AliWarning(" ATTENTION -> fBeamEnergy = 0\n");
 515   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
 516     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
 517     nDetSpecNRight, nDetSpecPRight);*/
 518
 519   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
 520   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
 521   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
 522   Double_t impPar=0.;
 523   //
 524   Float_t corrADCZEMHG = ZEM1ADCCorr[0] + ZEM2ADCCorr[0];
 525   //
 526   if(corrADCZEMHG > supValueZEM){
 527     nGenSpecNLeft  = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN1[0]));
 528     nGenSpecPLeft  = (Int_t) ((fRecoParam->GetfZPCen())->Eval(calibSumZP1[0]));
 529     nGenSpecLeft   = (Int_t) ((fRecoParam->GetfZDCCen())->Eval(calibSumZN1[0]+calibSumZP1[0]));
 530     nGenSpecNRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN2[0]));
 531     nGenSpecPRight = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZP2[0]));
 532     nGenSpecRight  = (Int_t) ((fRecoParam->GetfZNCen())->Eval(calibSumZN2[0]+calibSumZP2[0]));
 533     impPar  = (fRecoParam->GetfbCen())->Eval(calibSumZN1[0]+calibSumZP1[0]);
 534   }
 535   else if(corrADCZEMHG < infValueZEM){
 536     nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZNPer())->Eval(calibSumZN1[0]));
 537     nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZPPer())->Eval(calibSumZP1[0]));
 538     nGenSpecLeft  = (Int_t) ((fRecoParam->GetfZDCPer())->Eval(calibSumZN1[0]+calibSumZP1[0]));
 539     impPar   = (fRecoParam->GetfbPer())->Eval(calibSumZN1[0]+calibSumZP1[0]);
 540   }
 541   else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
 542     nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 543     nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 544     nGenSpecLeft  = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 545     impPar   =  (fRecoParam->GetfZEMb())->Eval(corrADCZEMHG);
 546   }
 547   //
 548   if(calibSumZN1[0]/maxValEZN1>1.)  nGenSpecNLeft = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 549   if(calibSumZP1[0]/maxValEZP1>1.)  nGenSpecPLeft = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 550   if((calibSumZN1[0]+calibSumZP1[0]/maxValEZDC1)>1.){
 551      nGenSpecLeft = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 552      impPar = (fRecoParam->GetfZEMb())->Eval(corrADCZEMHG);
 553   }
 554   if(calibSumZN2[0]/maxValEZN2>1.)  nGenSpecNRight = (Int_t) ((fRecoParam->GetfZEMn())->Eval(corrADCZEMHG));
 555   if(calibSumZP2[0]/maxValEZP2>1.)  nGenSpecPRight = (Int_t) ((fRecoParam->GetfZEMp())->Eval(corrADCZEMHG));
 556   if((calibSumZN2[0]+calibSumZP2[0]/maxValEZDC2)>1.) nGenSpecRight = (Int_t)((fRecoParam->GetfZEMsp())->Eval(corrADCZEMHG));
 557   //
 558   if(nGenSpecNLeft>125)    nGenSpecNLeft=125;
 559   else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
 560   if(nGenSpecPLeft>82)     nGenSpecPLeft=82;
 561   else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
 562   if(nGenSpecLeft>207)     nGenSpecLeft=207;
 563   else if(nGenSpecLeft<0)  nGenSpecLeft=0;
 564
 565   //  ---      Number of generated participants (from HIJING parameterization)
 566   Int_t nPart, nPartTotLeft, nPartTotRight;
 567   nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
 568   nPartTotLeft = 207-nGenSpecLeft;
 569   nPartTotRight = 207-nGenSpecRight;
 570   if(nPart<0) nPart=0;
 571   if(nPartTotLeft<0) nPartTotLeft=0;
 572   if(nPartTotRight<0) nPartTotRight=0;
 573   //
 574   // *** DEBUG ***
 575   /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
 576       "  calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
 577       calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
 578   printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
 579       "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n",
 580       nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft,
 581       nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
 582   printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
 583   */
 584
 585   // create the output tree
 586   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
 587                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
 588                   ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
 589                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
 590                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight,
 591                   nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
 592
 593   AliZDCReco* preco = &reco;
 594   const Int_t kBufferSize = 4000;
 595   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
 596
 597   // write the output tree
 598   clustersTree->Fill();
 599 }
 600
 601 //_____________________________________________________________________________
 602 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
 603 {
 604   // fill energies and number of participants to the ESD
 605
 606   AliZDCReco reco;
 607   AliZDCReco* preco = &reco;
 608   clustersTree->SetBranchAddress("ZDC", &preco);
 609
 610   clustersTree->GetEntry(0);
 611   //
 612   AliESDZDC * esdzdc = esd->GetESDZDC();
 613   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
 614   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
 615   for(Int_t i=0; i<5; i++){
 616      tZN1Ene[i] = reco.GetZN1HREnTow(i);
 617      tZN2Ene[i] = reco.GetZN2HREnTow(i);
 618      tZP1Ene[i] = reco.GetZP1HREnTow(i);
 619      tZP2Ene[i] = reco.GetZP2HREnTow(i);
 620      //
 621      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
 622      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
 623      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
 624      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
 625   }
 626   esdzdc->SetZN1TowerEnergy(tZN1Ene);
 627   esdzdc->SetZN2TowerEnergy(tZN2Ene);
 628   esdzdc->SetZP1TowerEnergy(tZP1Ene);
 629   esdzdc->SetZP2TowerEnergy(tZP2Ene);
 630   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
 631   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
 632   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
 633   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
 634   //
 635   esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
 636               reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
 637               reco.GetNPartLeft(), reco.GetNPartRight());
 638   //
 639
 640 }
 641
 642 //_____________________________________________________________________________
 643 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
 644 {
 645   // Setting the storage
 646
 647   Bool_t deleteManager = kFALSE;
 648
 649   AliCDBManager *manager = AliCDBManager::Instance();
 650   AliCDBStorage *defstorage = manager->GetDefaultStorage();
 651
 652   if(!defstorage || !(defstorage->Contains("ZDC"))){
 653      AliWarning("No default storage set or default storage doesn't contain ZDC!");
 654      manager->SetDefaultStorage(uri);
 655      deleteManager = kTRUE;
 656   }
 657
 658   AliCDBStorage *storage = manager->GetDefaultStorage();
 659
 660   if(deleteManager){
 661     AliCDBManager::Instance()->UnsetDefaultStorage();
 662     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
 663   }
 664
 665   return storage;
 666 }
 667
 668 //_____________________________________________________________________________
 669 AliZDCPedestals* AliZDCReconstructor::GetPedData() const
 670 {
 671
 672   // Getting pedestal calibration object for ZDC set
 673
 674   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
 675   if(!entry) AliFatal("No calibration data loaded!");
 676
 677   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
 678   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 679
 680   return calibdata;
 681 }
 682
 683 //_____________________________________________________________________________
 684 AliZDCCalib* AliZDCReconstructor::GetECalibData() const
 685 {
 686
 687   // Getting energy and equalization calibration object for ZDC set
 688
 689   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EMDCalib");
 690   if(!entry) AliFatal("No calibration data loaded!");
 691
 692   AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*>  (entry->GetObject());
 693   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 694
 695   return calibdata;
 696 }
 697