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