]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ZDC/AliZDCReconstructor.cxx
AOD particles with cluster index 0 were not checked for their shower shape, avoid...
[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx
index 279004930bd2a8344f8117934e38bdb8bf0d21f4..b95c7651e195c613bd55dec6ddd22a8982f4c951 100644 (file)
 #include "AliZDCReconstructor.h"
 #include "AliZDCPedestals.h"
 #include "AliZDCEnCalib.h"
+#include "AliZDCSaturationCalib.h"
 #include "AliZDCTowerCalib.h"
 #include "AliZDCMBCalib.h"
+#include "AliZDCTDCCalib.h"
 #include "AliZDCRecoParam.h"
 #include "AliZDCRecoParampp.h"
 #include "AliZDCRecoParamPbPb.h"
@@ -58,7 +60,9 @@ AliZDCMBCalib *AliZDCReconstructor::fgMBCalibData=0;  //calibration parameters f
 AliZDCReconstructor:: AliZDCReconstructor() :
   fPedData(GetPedestalData()),
   fEnCalibData(GetEnergyCalibData()),
+  fSatCalibData(GetSaturationCalibData()),
   fTowCalibData(GetTowerCalibData()),
+  fTDCCalibData(GetTDCCalibData()),
   fRecoMode(0),
   fBeamEnergy(0.),
   fNRun(0),
@@ -66,8 +70,7 @@ AliZDCReconstructor:: AliZDCReconstructor() :
   fPedSubMode(0),
   fSignalThreshold(7),
   fMeanPhase(0),
-  fESDZDC(NULL)
-{
+  fESDZDC(NULL){
   // **** Default constructor
 }
 
@@ -79,6 +82,7 @@ AliZDCReconstructor::~AliZDCReconstructor()
 //   if(fgRecoParam)    delete fgRecoParam;
    if(fPedData)      delete fPedData;    
    if(fEnCalibData)  delete fEnCalibData;
+   if(fSatCalibData)  delete fSatCalibData;
    if(fTowCalibData) delete fTowCalibData;
    if(fgMBCalibData) delete fgMBCalibData;
    if(fESDZDC)       delete fESDZDC;
@@ -106,12 +110,19 @@ void AliZDCReconstructor::Init()
   }*/
     
   fBeamEnergy = GetRunInfo()->GetBeamEnergy();
-  if(fBeamEnergy<0.01) AliWarning(" Beam energy value missing -> E_beam = 0");
-
+  if(fBeamEnergy<0.01){
+     AliWarning(" Beam energy value missing -> setting it to 1380 GeV ");
+     fBeamEnergy = 1380.;
+  }
+  
   if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
      ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)){
     fRecoMode=1;
   }
+  else if(((beamType.CompareTo("p-A"))==0) || ((beamType.CompareTo("A-p"))==0)
+     ||((beamType.CompareTo("P-A"))==0) || ((beamType.CompareTo("A-P"))==0)){
+    fRecoMode=1;
+  }
   else if((beamType.CompareTo("A-A")) == 0 || (beamType.CompareTo("AA")) == 0){
     fRecoMode=2;
     if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
@@ -122,14 +133,15 @@ void AliZDCReconstructor::Init()
 
   AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase"); 
   if (!entry) AliFatal("LHC clock-phase shift is not found in OCDB !");
-  AliLHCClockPhase *phaseLHC = (AliLHCClockPhase*)entry->GetObject();
-  // 4/2/2011 According to A. Di Mauro BEAM 2 measurement is more reliable 
-  // than BEAM1 and therefore laso than the average of the 2
-  fMeanPhase = phaseLHC->GetMeanPhaseB2();
-    
+  else{
+    AliLHCClockPhase *phaseLHC = (AliLHCClockPhase*)entry->GetObject();
+    // 4/2/2011 According to A. Di Mauro BEAM1 measurement is more reliable 
+    // than BEAM2 and therefore also than the average of the 2
+    fMeanPhase = phaseLHC->GetMeanPhaseB1();
+  }  
   if(fIsCalibrationMB==kFALSE)  
-    printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f + %1.0f GeV *****\n\n",
-       beamType.Data(), fBeamEnergy, fBeamEnergy);
+    AliInfo(Form("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f + %1.0f GeV *****\n\n",
+       beamType.Data(), fBeamEnergy, fBeamEnergy));
   
   // if EMD calibration run NO ENERGY CALIBRATION should be performed
   // pp-like reconstruction must be performed (E cailb. coeff. = 1)
@@ -138,6 +150,8 @@ void AliZDCReconstructor::Init()
     fBeamEnergy = 1380.;
   }
   
+  AliInfo(Form("\n   ZDC reconstruction mode %d (1 -> p-p, 2-> A-A)\n\n",fRecoMode));
+  
   fESDZDC = new AliESDZDC();
 
 }
@@ -157,6 +171,10 @@ void AliZDCReconstructor::Init(TString beamType, Float_t beamEnergy)
      ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)){
     fRecoMode=1;
   }
+  else if(((beamType.CompareTo("p-A"))==0) || ((beamType.CompareTo("A-p"))==0)
+     ||((beamType.CompareTo("P-A"))==0) || ((beamType.CompareTo("A-P"))==0)){
+    fRecoMode=1;
+  }
   else if((beamType.CompareTo("A-A")) == 0 || (beamType.CompareTo("AA")) == 0){
     fRecoMode=2;
     if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
@@ -165,13 +183,14 @@ void AliZDCReconstructor::Init(TString beamType, Float_t beamEnergy)
 
   AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase"); 
   if (!entry) AliFatal("LHC clock-phase shift is not found in OCDB !");
-  AliLHCClockPhase *phaseLHC = (AliLHCClockPhase*)entry->GetObject();
-  fMeanPhase = phaseLHC->GetMeanPhase();
-  
+  else{
+    AliLHCClockPhase *phaseLHC = (AliLHCClockPhase*)entry->GetObject();
+    fMeanPhase = phaseLHC->GetMeanPhase();
+  }
   fESDZDC = new AliESDZDC();
   
-  printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f + %1.0f GeV *****\n\n",
-       beamType.Data(), fBeamEnergy, fBeamEnergy);
+  AliInfo(Form("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f + %1.0f GeV *****\n\n",
+       beamType.Data(), fBeamEnergy, fBeamEnergy));
   
 }
 
@@ -201,18 +220,15 @@ void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) co
 
   // loop over digits
   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10]; 
-  Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2]; 
-  for(Int_t i=0; i<10; i++){
-     tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
-     if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
-  }  
+  for(Int_t i=0; i<10; i++) tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;  
+  Float_t dZEM1Corr[2]={0,0}, dZEM2Corr[2]={0,0}, sPMRef1[2]={0,0}, sPMRef2[2]={0,0}; 
   
   Int_t digNentries = digitsTree->GetEntries();
   Float_t ootDigi[kNch]; Int_t i=0;
   // -- Reading out-of-time signals (last kNch entries) for current event
   if(fPedSubMode==1){
     for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
-       if(i<=kNch) ootDigi[i] = digitsTree->GetEntry(iDigit);
+       if(i<=kNch) ootDigi[i-1] = digitsTree->GetEntry(iDigit);
        else AliWarning(" Can't read more out of time values: index>kNch !!!\n");
        i++;
     }
@@ -248,56 +264,38 @@ void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) co
     if(det == 1){ // *** ZNC
        tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
        tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
-       if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
-       if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
     }
     else if(det == 2){ // *** ZP1
        tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
        tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
-       if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
-       if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
     }
     else if(det == 3){
        if(quad == 1){      // *** ZEM1  
          dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg); 
          dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg); 
-         if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
-         if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
        }
        else if(quad == 2){  // *** ZEM2
          dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg); 
          dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg); 
-         if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
-         if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
        }
     }
     else if(det == 4){  // *** ZN2
        tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
        tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
-       if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
-       if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
    }
     else if(det == 5){  // *** ZP2 
        tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
        tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
-       if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
-       if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
     }
    }
    else{ // Reference PMs
      if(det == 1){
        sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
        sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
-       // Ch. debug
-       if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
-       if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
      }
      else if(det == 4){
        sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
        sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
-       // Ch. debug
-       if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
-       if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
      }
    }
 
@@ -374,11 +372,8 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
   }
   
   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10]; 
-  Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2]; 
-  for(Int_t i=0; i<10; i++){
-     tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
-     if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
-  }  
+  for(Int_t i=0; i<10; i++) tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
+  Float_t dZEM1Corr[2]={0,0}, dZEM2Corr[2]={0,0}, sPMRef1[2]={0,0}, sPMRef2[2]={0,0}; 
 
   Bool_t isScalerOn=kFALSE;
   Int_t jsc=0, itdc=0, iprevtdc=-1, ihittdc=0;
@@ -426,36 +421,15 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
        //  **** Pb-Pb data taking 2010 -> subtracting some ch. from correlation ****
        // Not interested in o.o.t. signals (ADC modules 2, 3)
        //if(adcMod == 2 || adcMod == 3) continue;
-       if(((det==1 && quad==0) || (det==2 && quad==2) || (det==3)) && (quad != 5)){
-         if(det == 1){
-           if(adcMod==0 || adcMod==1){
-            if(gain==0) adcZN1[quad] = rawData.GetADCValue();
-             else adcZN1lg[quad] = rawData.GetADCValue();
-          }
-          else if(adcMod==2 || adcMod==3){
-            if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
-             else adcZN1ootlg[quad] = rawData.GetADCValue();
-          }
+       //  **** Pb-Pb data taking 2011 -> subtracting only ZEM from correlation ****
+       if(det==3){
+        if(adcMod==0 || adcMod==1){
+          if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
+           else adcZEMlg[quad-1] = rawData.GetADCValue();
         }
-        else if(det == 2){
-          if(adcMod==0 || adcMod==1){
-            if(gain==0) adcZP2[quad] = rawData.GetADCValue();
-             else adcZP2lg[quad] = rawData.GetADCValue();
-          }
-          else if(adcMod==2 || adcMod==3){
-            if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
-             else adcZP2ootlg[quad] = rawData.GetADCValue();
-          }
-        }
-        else if(det == 3){
-          if(adcMod==0 || adcMod==1){
-            if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
-             else adcZEMlg[quad-1] = rawData.GetADCValue();
-          }
-          else if(adcMod==2 || adcMod==3){ 
-            if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
-             else adcZEMootlg[quad-1] = rawData.GetADCValue();
-          }
+        else if(adcMod==2 || adcMod==3){ 
+          if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
+           else adcZEMootlg[quad-1] = rawData.GetADCValue();
         }
        }
        // When oot values are read the ADC modules 2, 3 can be skipped!!!
@@ -463,12 +437,12 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
        
        // *************************************************************************
        if(quad != 5){ // ZDCs (not reference PTMs)
-        if(det==1 && quad!=0){    
+        if(det==1){    
           pedindex = quad;
           if(gain == 0) tZN1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
           else tZN1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
         }
-        else if(det==2 && quad!=2){ 
+        else if(det==2){ 
           pedindex = quad+5;
           if(gain == 0) tZP1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
           else tZP1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
@@ -499,7 +473,7 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
          pedindex = (det-1)/3 + 22;
          if(det == 1){
            if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
-        else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
+          else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
          }
          else if(det == 4){
            if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
@@ -582,7 +556,7 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
    }// ADC DATA
    // ***************************** Reading Scaler
    else if(rawData.GetADCModule()==kScalerGeo){
-     if(rawData.IsScalerWord()==kTRUE && rawData.IsScEventGood()==kTRUE){
+     if(rawData.IsScalerWord()==kTRUE){
        isScalerOn = kTRUE;
        scalerData[jsc] = rawData.GetTriggerCount();
        // Ch. debug
@@ -597,9 +571,9 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
        if(itdc==iprevtdc) ihittdc++;
        else ihittdc=0;
        iprevtdc=itdc;
-       tdcData[itdc][ihittdc] = rawData.GetZDCTDCDatum();
+       if(ihittdc<4) tdcData[itdc][ihittdc] = rawData.GetZDCTDCDatum();
        // Ch. debug
-       //printf("   Reconstructed TDC[%d, %d] %d  ",itdc, ihittdc, tdcData[itdc][ihittdc]);
+       //if(ihittdc==0) printf("   TDC%d %d  ",itdc, tdcData[itdc][ihittdc]);
    }// ZDC TDC DATA
    // ***************************** Reading PU
    else if(rawData.GetADCModule()==kPUGeo){
@@ -639,16 +613,16 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
     sPMRef2[0] = pmRef[0]   - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
     sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
   }
-  else{
-    //  **** Pb-Pb data taking 2010 -> subtracting some ch. from correlation ****
-    tZN1Corr[0] = adcZN1[0] - (corrCoeff1[0]*adcZN1oot[0]+corrCoeff0[0]);
-    tZN1Corr[5] = adcZN1lg[0] - (corrCoeff1[kNch]*adcZN1ootlg[0]+corrCoeff0[kNch]);
+  if(fPedSubMode==0 && fRecoMode==2){
+    //  **** Pb-Pb data taking 2011 -> subtracting some ch. from correlation ****
+    //tZN1Corr[0] = adcZN1[0] - (corrCoeff1[0]*adcZN1oot[0]+corrCoeff0[0]);
+    //tZN1Corr[5] = adcZN1lg[0] - (corrCoeff1[kNch]*adcZN1ootlg[0]+corrCoeff0[kNch]);
     // Ch. debug
     //printf(" adcZN1 %d  adcZN1oot %d tZN1Corr %1.2f \n", adcZN1[0],adcZN1oot[0],tZN1Corr[0]);
     //printf(" adcZN1lg %d  adcZN1ootlg %d tZN1Corrlg %1.2f \n", adcZN1lg[0],adcZN1ootlg[0],tZN1Corr[5]);
     //
-    tZP1Corr[2] = adcZP1[2] - (corrCoeff1[2+5]*adcZP1oot[2]+corrCoeff0[2+5]);
-    tZP1Corr[2+5] = adcZP1lg[2] - (corrCoeff1[2+5+kNch]*adcZP1ootlg[2]+corrCoeff0[2+5+kNch]);
+    //tZP1Corr[2] = adcZP1[2] - (corrCoeff1[2+5]*adcZP1oot[2]+corrCoeff0[2+5]);
+    //tZP1Corr[2+5] = adcZP1lg[2] - (corrCoeff1[2+5+kNch]*adcZP1ootlg[2]+corrCoeff0[2+5+kNch]);
     //
     dZEM1Corr[0] = adcZEM[0]   - (corrCoeff1[10]*adcZEMoot[0]+corrCoeff0[10]);
     dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[10+kNch]*adcZEMootlg[0]+corrCoeff0[10+kNch]);
@@ -656,6 +630,15 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
     dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[11+kNch]*adcZEMootlg[1]+corrCoeff0[11+kNch]);
     // *************************************************************************
   }
+  /*else if(fPedSubMode==0 && fRecoMode==1){
+    //  **** p-p data taking 2011 -> temporary patch to overcome DA problem ****
+    //
+    dZEM1Corr[0] = adcZEM[0]   - meanPed[10];
+    dZEM1Corr[1] = adcZEMlg[0] - meanPed[10+kNch];
+    dZEM2Corr[0] = adcZEM[1]   - meanPed[11];
+    dZEM2Corr[1] = adcZEMlg[1] - meanPed[11+kNch];
+        // *************************************************************************
+  }*/
     
   if(fRecoMode==1) // p-p data
     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr, 
@@ -749,10 +732,14 @@ void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree,
      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji); 
   }
   // --- Energy calibration factors ------------------------------------
-  Float_t calibEne[6];
+  Float_t calibEne[6], calibSatZNA[4], calibSatZNC[4];
   // **** Energy calibration coefficient set to 1 
   // **** (no trivial way to calibrate in p-p runs)
   for(Int_t ij=0; ij<6; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
+  for(Int_t ij=0; ij<4; ij++){
+    calibSatZNA[ij] = fSatCalibData->GetZNASatCalib(ij);
+    calibSatZNC[ij] = fSatCalibData->GetZNCSatCalib(ij);
+  }
   
   // ******    Equalization of detector responses
   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
@@ -782,6 +769,24 @@ void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree,
        equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
   printf(" ----------------------------------------\n");*/
   
+  //  *** p-A RUN 2013 -> new calibration object
+  //      to take into account saturation in ZN PMC
+  //   -> 5th order pol. fun. to be applied BEFORE en. calibration 
+  equalTowZN1[0] = equalTowZN1[0] + calibSatZNC[0]*equalTowZN1[0]*equalTowZN1[0] +
+       calibSatZNC[1]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0] +
+       calibSatZNC[2]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0] +
+       calibSatZNC[3]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0];
+  equalTowZN2[0] = equalTowZN2[0] + calibSatZNA[0]*equalTowZN2[0]*equalTowZN2[0] +
+       calibSatZNA[1]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0] +
+       calibSatZNA[2]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0] +
+       calibSatZNA[3]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0];
+
+  // Ch. debug
+  /*printf("\n ------------- SATURATION CORRECTION -------------\n");
+  printf(" ZNC PMC %1.2f\n", equalTowZN1[0]);
+  printf(" ZNA PMC %1.2f\n", equalTowZN2[0]);
+  printf(" ----------------------------------------\n");*/
+  
   // ******    Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
   for(Int_t gi=0; gi<5; gi++){
@@ -821,12 +826,11 @@ void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree,
      calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
   }
   //
-  Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
+  Float_t calibZEM1[]={0,0}, calibZEM2[]={0,0};
   calibZEM1[0] = corrADCZEM1[0]*calibEne[4];
   calibZEM1[1] = corrADCZEM1[1]*calibEne[4];
   calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
-  for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
   // Ch. debug
   /*printf("\n ------------- CALIBRATION -------------\n");
   printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
@@ -848,6 +852,7 @@ void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree,
   Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
   Double_t impPar=0., impPar1=0., impPar2=0.;
   
+  Bool_t energyFlag = kFALSE;
   // create the output tree
   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2, 
                   calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2, 
@@ -856,7 +861,7 @@ void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree,
                   nGenSpec, nGenSpecLeft, nGenSpecRight, 
                   nPart, nPartTotLeft, nPartTotRight, 
                   impPar, impPar1, impPar2,
-                  recoFlag, isScalerOn, scaler, tdcData);
+                  recoFlag, energyFlag, isScalerOn, scaler, tdcData);
                  
   const Int_t kBufferSize = 4000;
   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
@@ -944,10 +949,14 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
      equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji); 
   }
   // --- Energy calibration factors ------------------------------------
-  Float_t calibEne[6];
-  // The energy calibration object already takes into account of E_beam 
-  // -> the value from the OCDB can be directly used (Jul 2010)
+  Float_t calibEne[6], calibSatZNA[4], calibSatZNC[4];
+  // **** Energy calibration coefficient set to 1 
+  // **** (no trivial way to calibrate in p-p runs)
   for(Int_t ij=0; ij<6; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
+  for(Int_t ij=0; ij<4; ij++){
+    calibSatZNA[ij] = fSatCalibData->GetZNASatCalib(ij);
+    calibSatZNC[ij] = fSatCalibData->GetZNCSatCalib(ij);
+  }
   
   // ******    Equalization of detector responses
   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
@@ -979,6 +988,18 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
   printf(" ----------------------------------------\n");
 */
   
+  //  *** p-A RUN 2013 -> new calibration object
+  //      to take into account saturation in ZN PMC
+  //   -> 5th order pol. fun. to be applied BEFORE en. calibration 
+  equalTowZN1[0] = equalTowZN1[0] + calibSatZNC[0]*equalTowZN1[0]*equalTowZN1[0] +
+       calibSatZNC[1]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0] +
+       calibSatZNC[2]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0] +
+       calibSatZNC[3]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0]*equalTowZN1[0];
+  equalTowZN2[0] = equalTowZN2[0] + calibSatZNA[0]*equalTowZN2[0]*equalTowZN2[0] +
+       calibSatZNA[1]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0] +
+       calibSatZNA[2]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0] +
+       calibSatZNA[3]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0]*equalTowZN2[0];
+  
   // ******    Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
   for(Int_t gi=0; gi<5; gi++){
@@ -1006,12 +1027,11 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
   calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
   calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
   //
-  Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
-  calibZEM1[0] = corrADCZEM1[0]*calibEne[4]*8.;
+  Float_t calibZEM1[]={0,0}, calibZEM2[]={0,0};
+  calibZEM1[0] = corrADCZEM1[0]*calibEne[4];
   calibZEM1[1] = corrADCZEM1[1]*calibEne[4];
-  calibZEM2[0] = corrADCZEM2[0]*calibEne[5]*8.;
+  calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
   calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
-  for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
     
   // ******    Energy calibration of detector responses
   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
@@ -1059,20 +1079,20 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
   Double_t b=0., bA=0., bC=0.;
   
   if(fIsCalibrationMB == kFALSE){
-    // ******  Reconstruction parameters ------------------ 
-    if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
-    if(!fgRecoParam){  
-      AliError("  RecoParam object not retrieved correctly: not reconstructing event!!!");
-      return;
-    }
-    TH1D* hNpartDist = fgRecoParam->GethNpartDist();
-    TH1D*   hbDist = fgRecoParam->GethbDist();    
-    Float_t  fClkCenter = fgRecoParam->GetClkCenter();
-    if(!hNpartDist || !hbDist){  
-       AliError("Something wrong in Glauber MC histos got from AliZDCREcoParamPbPb: NO EVENT RECO FOR ZDC DATA!!!\n\n");
-       return;
-    }
-     
+   // ******   Reconstruction parameters ------------------ 
+   if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
+   if(!fgRecoParam){  
+     AliError("  RecoParam object not retrieved correctly: not reconstructing ZDC event!!!");
+     return;
+   }
+   TH1D* hNpartDist = fgRecoParam->GethNpartDist();
+   TH1D*   hbDist = fgRecoParam->GethbDist();   
+   Float_t  fClkCenter = fgRecoParam->GetClkCenter();
+   if(!hNpartDist || !hbDist){  
+      AliError("Something wrong in Glauber MC histos got from AliZDCREcoParamPbPb: NO EVENT RECO FOR ZDC DATA!!!\n\n");
+      //return;
+   }
+   else{  
     if(!fgMBCalibData) fgMBCalibData = const_cast<AliZDCMBCalib*>(GetMBCalibData()); 
     TH2F *hZDCvsZEM  = fgMBCalibData->GethZDCvsZEM();
     TH2F *hZDCCvsZEM = fgMBCalibData->GethZDCCvsZEM();
@@ -1280,16 +1300,17 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
     
     delete line; 
     delete lineC;  delete lineA;
-
+   }
   } // ONLY IF fIsCalibrationMB==kFALSE
   
+  Bool_t energyFlag = kTRUE;  
   AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2, 
                  calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2, 
                  calibZEM1, calibZEM2, sPMRef1, sPMRef2,
                  nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight, 
                  nGenSpec, nGenSpecA, nGenSpecC, 
                  nPart, nPartA, nPartC, b, bA, bC,
-                 recoFlag, isScalerOn, scaler, tdcData);
+                 recoFlag, energyFlag, isScalerOn, scaler, tdcData);
                    
   const Int_t kBufferSize = 4000;
   clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
@@ -1305,6 +1326,13 @@ void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd)
 {
   // fill energies and number of participants to the ESD
 
+  // Retrieving TDC calibration data  
+  // Parameters for TDC centering around zero
+  int const knTDC = 6;
+  Float_t tdcOffset[knTDC];
+  for(Int_t jj=0; jj<knTDC; jj++) tdcOffset[jj] = fTDCCalibData->GetMeanTDC(jj);
+  //fTDCCalibData->Print("");
+
   AliZDCReco reco;
   AliZDCReco* preco = &reco;
   clustersTree->SetBranchAddress("ZDC", &preco);
@@ -1355,24 +1383,45 @@ void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd)
     fESDZDC->SetZDCScaler(counts);
   }    
   
-  // Writing TDC data into ZDC ESDs
-  Int_t tdcValues[32][4]; 
-  Float_t tdcCorrected[32][4];
+  Int_t tdcValues[32][4] = {{0,}}; 
+  Float_t tdcCorrected[32][4] = {{9999.,}};
   for(Int_t jk=0; jk<32; jk++){
     for(Int_t lk=0; lk<4; lk++){
       tdcValues[jk][lk] = reco.GetZDCTDCData(jk, lk);
+      //
+      if(jk==8 && TMath::Abs(tdcValues[jk][lk])>1e-09)      fESDZDC->SetZEM1TDChit(kTRUE);
+      else if(jk==9 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZEM2TDChit(kTRUE);
+      else if(jk==10 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZNCTDChit(kTRUE);
+      else if(jk==11 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZPCTDChit(kTRUE);
+      else if(jk==12 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZNATDChit(kTRUE);
+      else if(jk==13 && TMath::Abs(tdcValues[jk][lk])>1e-09) fESDZDC->SetZPATDChit(kTRUE);
+      //Ch debug
+      //if((jk>=8 && jk<=13 && lk==0) || jk==15) printf(" *** ZDC: tdc%d =  %d = %f ns \n",jk,tdcValues[jk][lk],0.025*tdcValues[jk][lk]);
     }
   }
+  
+  // Writing TDC data into ZDC ESDs
   // 4/2/2011 -> Subtracting L0 (tdcValues[15]) instead of ADC gate 
   // we try to keep the TDC oscillations as low as possible!
   for(Int_t jk=0; jk<32; jk++){
     for(Int_t lk=0; lk<4; lk++){
-      if(tdcValues[jk][lk]!=0.) tdcCorrected[jk][lk] = 0.025*(tdcValues[jk][lk]-tdcValues[15][0])+fMeanPhase;
+      if(tdcValues[jk][lk]!=0.){
+        // Feb2013 _-> TDC correct entry is there ONLY IF tdc has a hit!
+        if(TMath::Abs(tdcValues[jk][lk])>1e-09){
+          tdcCorrected[jk][lk] = 0.025*(tdcValues[jk][lk]-tdcValues[15][0])+fMeanPhase;
+           // Sep 2011: TDC ch. from 8 to 13 centered around 0 using OCDB 
+          if(jk>=8 && jk<=13) tdcCorrected[jk][lk] =  tdcCorrected[jk][lk] - tdcOffset[jk-8];
+          //Ch. debug
+          //if(jk>=8 && jk<=13) printf(" *** tdcOffset%d %f  tdcCorr%d %f \n",jk,tdcOffset[jk-8],tdcCorrected[jk][lk]);
+        }
+      }
     }
   }
+
   fESDZDC->SetZDCTDCData(tdcValues);
   fESDZDC->SetZDCTDCCorrected(tdcCorrected);
-  fESDZDC->AliESDZDC::SetBit(AliESDZDC::kCorrectedTDCFilled, kTRUE);
+  fESDZDC->AliESDZDC::SetBit(AliESDZDC::kCorrectedTDCFilled, reco.GetEnergyFlag());
+  fESDZDC->AliESDZDC::SetBit(AliESDZDC::kEnergyCalibratedSignal, kTRUE);
   
   if(esd) esd->SetZDCData(fESDZDC);
 }
@@ -1408,14 +1457,16 @@ AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
 {
 
   // Getting pedestal calibration object for ZDC set
-
+  AliZDCPedestals *calibdata = 0x0;
   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
   if(!entry) AliFatal("No calibration data loaded!");
-  entry->SetOwner(kFALSE);
+  else{
+    entry->SetOwner(kFALSE);
 
-  AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
-  if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+    calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
+    if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 
+  }
   return calibdata;
 }
 
@@ -1424,14 +1475,32 @@ AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
 {
 
   // Getting energy and equalization calibration object for ZDC set
-
+  AliZDCEnCalib *calibdata = 0x0;
   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
   if(!entry) AliFatal("No calibration data loaded!");  
-  entry->SetOwner(kFALSE);
+  else{
+    entry->SetOwner(kFALSE);
 
-  AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
-  if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+    calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
+    if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+  }
+  return calibdata;
+}
 
+//_____________________________________________________________________________
+AliZDCSaturationCalib* AliZDCReconstructor::GetSaturationCalibData() const
+{
+
+  // Getting energy and equalization calibration object for ZDC set
+  AliZDCSaturationCalib *calibdata = 0x0;
+  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/SaturationCalib");
+  if(!entry) AliFatal("No calibration data loaded!");  
+  else{
+    entry->SetOwner(kFALSE);
+
+    calibdata = dynamic_cast<AliZDCSaturationCalib*> (entry->GetObject());
+    if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+  }
   return calibdata;
 }
 
@@ -1440,14 +1509,15 @@ AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
 {
 
   // Getting energy and equalization calibration object for ZDC set
-
+  AliZDCTowerCalib *calibdata = 0x0;
   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
   if(!entry) AliFatal("No calibration data loaded!");  
-  entry->SetOwner(kFALSE);
-
-  AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
-  if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+  else{
+    entry->SetOwner(kFALSE);
 
+    calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
+    if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+  }
   return calibdata;
 }
 
@@ -1456,13 +1526,32 @@ AliZDCMBCalib* AliZDCReconstructor::GetMBCalibData() const
 {
 
   // Getting energy and equalization calibration object for ZDC set
-
+  AliZDCMBCalib *calibdata = 0x0;
   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/MBCalib");
   if(!entry) AliFatal("No calibration data loaded!");  
-  entry->SetOwner(kFALSE);
+  else{
+    entry->SetOwner(kFALSE);
+
+    calibdata = dynamic_cast<AliZDCMBCalib*> (entry->GetObject());
+    if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+  }
+  return calibdata;
+}
+
+//_____________________________________________________________________________
+AliZDCTDCCalib* AliZDCReconstructor::GetTDCCalibData() const
+{
+
+  // Getting TDC object for ZDC 
+  AliZDCTDCCalib *calibdata = 0x0;
+  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TDCCalib");
+  if(!entry) AliFatal("No calibration data loaded!");  
+  else{
+    entry->SetOwner(kFALSE);
 
-  AliZDCMBCalib *calibdata = dynamic_cast<AliZDCMBCalib*> (entry->GetObject());
-  if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+    calibdata = dynamic_cast<AliZDCTDCCalib*> (entry->GetObject());
+    if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 
+  }
   return calibdata;
 }