]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ZDC/AliZDCReconstructor.cxx
Fiexs needed for Pb-Pb reconstruction
[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx
index 0789612ee4e6c83fee5297d8a037f8aab540fcc9..c8189eace99ca9cea06dc541645d8c3064f43541 100644 (file)
@@ -47,6 +47,7 @@
 #include "AliZDCRecoParampp.h"
 #include "AliZDCRecoParamPbPb.h"
 #include "AliRunInfo.h"
+#include "AliLHCClockPhase.h"
 
 
 ClassImp(AliZDCReconstructor)
@@ -63,7 +64,9 @@ AliZDCReconstructor:: AliZDCReconstructor() :
   fNRun(0),
   fIsCalibrationMB(kFALSE),
   fPedSubMode(0),
-  fSignalThreshold(7)
+  fSignalThreshold(7),
+  fMeanPhase(0),
+  fESDZDC(NULL)
 {
   // **** Default constructor
 }
@@ -77,15 +80,15 @@ AliZDCReconstructor::~AliZDCReconstructor()
    if(fPedData)      delete fPedData;    
    if(fEnCalibData)  delete fEnCalibData;
    if(fTowCalibData) delete fTowCalibData;
-   if(fgMBCalibData)  delete fgMBCalibData;
+   if(fgMBCalibData) delete fgMBCalibData;
+   if(fESDZDC)       delete fESDZDC;
 }
 
 //____________________________________________________________________________
 void AliZDCReconstructor::Init()
 {
-  // Setting reconstruction mode
-  // Getting beam type and beam energy from GRP calibration object
-  
+  // Setting reconstruction parameters
+    
   TString runType = GetRunInfo()->GetRunType();
   if((runType.CompareTo("CALIBRATION_MB")) == 0){
     fIsCalibrationMB = kTRUE;
@@ -101,6 +104,9 @@ void AliZDCReconstructor::Init()
     AliError("\t UNKNOWN beam type\n");
     return;
   }*/
+    
+  fBeamEnergy = GetRunInfo()->GetBeamEnergy();
+  if(fBeamEnergy<0.01) AliWarning(" Beam energy value missing -> E_beam = 0");
 
   if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
      ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)){
@@ -108,15 +114,65 @@ void AliZDCReconstructor::Init()
   }
   else if((beamType.CompareTo("A-A")) == 0 || (beamType.CompareTo("AA")) == 0){
     fRecoMode=2;
+    if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
+    if(fgRecoParam){
+      fgRecoParam->SetGlauberMCDist(fBeamEnergy);      
+    } 
   }
-    
-  fBeamEnergy = GetRunInfo()->GetBeamEnergy();
-  if(fBeamEnergy<0.01) AliWarning(" Beam energy value missing -> E_beam = 0");
+
+  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 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);
   
+  // if EMD calibration run NO ENERGY CALIBRATION should be performed
+  // pp-like reconstruction must be performed (E cailb. coeff. = 1)
+  if((runType.CompareTo("CALIBRATION_EMD")) == 0){
+    fRecoMode=1; 
+    fBeamEnergy = 1380.;
+  }
+  
+  fESDZDC = new AliESDZDC();
+
+}
+
+
+//____________________________________________________________________________
+void AliZDCReconstructor::Init(TString beamType, Float_t beamEnergy)
+{
+  // Setting reconstruction mode
+  // Needed to work in the HLT framework
+  
+  fIsCalibrationMB = kFALSE;
+     
+  fBeamEnergy = beamEnergy;
+  
+  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("A-A")) == 0 || (beamType.CompareTo("AA")) == 0){
+    fRecoMode=2;
+    if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
+    if( fgRecoParam ) fgRecoParam->SetGlauberMCDist(fBeamEnergy);      
+  }    
+
+  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();
+  
+  fESDZDC = new AliESDZDC();
+  
+  printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f + %1.0f GeV *****\n\n",
+       beamType.Data(), fBeamEnergy, fBeamEnergy);
+  
 }
 
 //_____________________________________________________________________________
@@ -257,10 +313,10 @@ void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) co
   }//digits loop
  
   UInt_t counts[32];
-  Int_t  tdc[32];
+  Int_t  tdc[32][4];
   for(Int_t jj=0; jj<32; jj++){
     counts[jj]=0;
-    tdc[jj]=0;
+    for(Int_t ii=0; ii<4; ii++) tdc[jj][ii]=0;
   }
   
   Int_t  evQualityBlock[4] = {1,0,0,0};
@@ -297,6 +353,7 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
   for(Int_t jj=0; jj<2*kNch; jj++){
      corrCoeff0[jj] =  fPedData->GetPedCorrCoeff0(jj);
      corrCoeff1[jj] =  fPedData->GetPedCorrCoeff1(jj);
+     //printf("  %d   %1.4f  %1.4f\n", jj,corrCoeff0[jj],corrCoeff1[jj]);
   }
 
   Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
@@ -324,25 +381,25 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
   }  
 
   Bool_t isScalerOn=kFALSE;
-  Int_t jsc=0, itdc=0;
+  Int_t jsc=0, itdc=0, iprevtdc=-1, ihittdc=0;
   UInt_t scalerData[32];
-  Int_t tdcData[32];   
+  Int_t tdcData[32][4];        
   for(Int_t k=0; k<32; k++){
     scalerData[k]=0;
-    tdcData[k]=0;
+    for(Int_t i=0; i<4; i++) tdcData[k][i]=0;
   }
   
+  
   Int_t  evQualityBlock[4] = {1,0,0,0};
   Int_t  triggerBlock[4] = {0,0,0,0};
   Int_t  chBlock[3] = {0,0,0};
   UInt_t puBits=0;
 
-  //fNRun = (Int_t) rawReader->GetRunNumber();
   Int_t kFirstADCGeo=0, kLastADCGeo=3, kScalerGeo=8, kZDCTDCGeo=4, kPUGeo=29;
   //Int_t kTrigScales=30, kTrigHistory=31;
 
   // loop over raw data
-  rawReader->Reset();
+  //rawReader->Reset();
   AliZDCRawStream rawData(rawReader);
   while(rawData.Next()){
    
@@ -366,21 +423,47 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
       //
       // Mean pedestal value subtraction -------------------------------------------------------
       if(fPedSubMode == 0){
+       //  **** 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==3))){
+         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();
+          }
+        }
+        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();
+          }
+        }
+       }
+       // When oot values are read the ADC modules 2, 3 can be skipped!!!
        if(adcMod == 2 || adcMod == 3) continue;
-       //
+       
+       // *************************************************************************
        if(quad != 5){ // ZDCs (not reference PTMs)
-        if(det == 1){    
+        if(det==1 && quad!=0){    
           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){ 
+        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]); 
         }
-        else if(det == 3){ 
+        /*else if(det == 3){ 
           pedindex = quad+9;
           if(quad==1){     
             if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
@@ -390,7 +473,7 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
             if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
             else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
           }
-        }
+        }*/
         else if(det == 4){      
           pedindex = quad+12;
           if(gain == 0) tZN2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
@@ -500,11 +583,13 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
    }// VME SCALER DATA
    // ***************************** Reading ZDC TDC
    else if(rawData.GetADCModule()==kZDCTDCGeo && rawData.IsZDCTDCDatum()==kTRUE){
-       tdcData[itdc] = rawData.GetZDCTDCDatum();
+       itdc = rawData.GetChannel(); 
+       if(itdc==iprevtdc) ihittdc++;
+       else ihittdc=0;
+       iprevtdc=itdc;
+       tdcData[itdc][ihittdc] = rawData.GetZDCTDCDatum();
        // Ch. debug
-       //printf("   Reconstructed VME Scaler: %d %d  ",jsc,scalerData[jsc]);
-       //
-       itdc++;
+       //printf("   Reconstructed TDC[%d, %d] %d  ",itdc, ihittdc, tdcData[itdc][ihittdc]);
    }// ZDC TDC DATA
    // ***************************** Reading PU
    else if(rawData.GetADCModule()==kPUGeo){
@@ -533,49 +618,34 @@ void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTr
        //
        tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
        tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
-       // 0---------0 Ch. debug 0---------0
-/*       printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
-       printf("\tCorrCoeff0\tCorrCoeff1\n");
-       printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
-       printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
-       printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
-       printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
-       printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
-       printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
-       printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
-       printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
-       
-       printf("ZN1 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
-                       adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
-       printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
-                       adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
-       //
-       printf("ZP1 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
-                       adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
-       printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
-                       adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
-       //
-       printf("ZN2 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
-                       adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
-       printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
-                       adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
-       //
-       printf("ZP2 ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
-                       adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
-       printf(" lg ->  rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
-                       adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
-*/
     }
-    dZEM1Corr[0] = adcZEM[0]   - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
-    dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
-    dZEM2Corr[0] = adcZEM[1]   - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
-    dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
+    dZEM1Corr[0] = adcZEM[0]   - (corrCoeff1[10]*adcZEMoot[0]+corrCoeff0[10]);
+    dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[10+kNch]*adcZEMootlg[0]+corrCoeff0[10+kNch]);
+    dZEM2Corr[0] = adcZEM[1]   - (corrCoeff1[11]*adcZEMoot[1]+corrCoeff0[11]);
+    dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[11+kNch]*adcZEMootlg[1]+corrCoeff0[11+kNch]);
     //
     sPMRef1[0] = pmRef[0]   - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
     sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
     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]);
+    // 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]);
+    //
+    dZEM1Corr[0] = adcZEM[0]   - (corrCoeff1[10]*adcZEMoot[0]+corrCoeff0[10]);
+    dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[10+kNch]*adcZEMootlg[0]+corrCoeff0[10+kNch]);
+    dZEM2Corr[0] = adcZEM[1]   - (corrCoeff1[11]*adcZEMoot[1]+corrCoeff0[11]);
+    dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[11+kNch]*adcZEMootlg[1]+corrCoeff0[11+kNch]);
+    // *************************************************************************
+  }
     
   if(fRecoMode==1) // p-p data
     ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr, 
@@ -595,7 +665,7 @@ void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree,
        const Float_t* const corrADCZN2, const Float_t* const corrADCZP2,
        const Float_t* const corrADCZEM1, const Float_t* const corrADCZEM2,
        Float_t* sPMRef1, Float_t* sPMRef2, Bool_t isScalerOn, UInt_t* scaler, 
-       Int_t* tdcData, const Int_t* const evQualityBlock, 
+       Int_t tdcData[32][4], const Int_t* const evQualityBlock, 
        const Int_t* const triggerBlock, const Int_t* const chBlock, UInt_t puBits) const
 {
   // ****************** Reconstruct one event ******************
@@ -768,6 +838,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, 
@@ -776,12 +847,13 @@ 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);
   // write the output tree
   clustersTree->Fill();
+  delete reco;
 }
 
 //_____________________________________________________________________________
@@ -790,7 +862,7 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
        const Float_t* const corrADCZN2, const Float_t* const corrADCZP2,
        const Float_t* const corrADCZEM1, const Float_t* const corrADCZEM2,
        Float_t* sPMRef1, Float_t* sPMRef2, Bool_t isScalerOn, UInt_t* scaler, 
-       Int_t* tdcData, const Int_t* const evQualityBlock, 
+       Int_t tdcData[32][4], const Int_t* const evQualityBlock, 
        const Int_t* const triggerBlock, const Int_t* const chBlock, UInt_t puBits) const
 {
   // ****************** Reconstruct one event ******************
@@ -936,15 +1008,15 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
   for(Int_t gi=0; gi<5; gi++){
      // High gain chain
-     calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]*8.;
-     calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]*8.;
-     calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]*8.;
-     calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]*8.;
+     calibTowZN1[gi] = equalTowZN1[gi]*2*calibEne[0]*8.;
+     calibTowZP1[gi] = equalTowZP1[gi]*2*calibEne[1]*8.;
+     calibTowZN2[gi] = equalTowZN2[gi]*2*calibEne[2]*8.;
+     calibTowZP2[gi] = equalTowZP2[gi]*2*calibEne[3]*8.;
      // Low gain chain
-     calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
-     calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
-     calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
-     calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
+     calibTowZN1[gi+5] = equalTowZN1[gi+5]*2*calibEne[0];
+     calibTowZP1[gi+5] = equalTowZP1[gi+5]*2*calibEne[1];
+     calibTowZN2[gi+5] = equalTowZN2[gi+5]*2*calibEne[2];
+     calibTowZP2[gi+5] = equalTowZP2[gi+5]*2*calibEne[3];
   }
 
   // Ch. debug
@@ -969,9 +1041,9 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
     nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
   }
   else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
-  printf("\n\t AliZDCReconstructor -> fBeamEnergy %1.0f: nDetSpecNsideA %d, nDetSpecPsideA %d,"
+  /*printf("\n\t AliZDCReconstructor -> fBeamEnergy %1.0f: nDetSpecNsideA %d, nDetSpecPsideA %d,"
     " nDetSpecNsideC %d, nDetSpecPsideC %d\n",fBeamEnergy,nDetSpecNLeft, nDetSpecPLeft, 
-    nDetSpecNRight, nDetSpecPRight);
+    nDetSpecNRight, nDetSpecPRight);*/
   
   Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
   Int_t nPart=0, nPartA=0, nPartC=0;
@@ -979,20 +1051,26 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
   
   if(fIsCalibrationMB == kFALSE){
     // ******  Reconstruction parameters ------------------ 
-    if(!fgMBCalibData) fgMBCalibData = const_cast<AliZDCMBCalib*>(GetMBCalibData()); 
-    if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam()); 
-    fgRecoParam->SetGlauberMCDist(fBeamEnergy);
+    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;
+    }
      
+    if(!fgMBCalibData) fgMBCalibData = const_cast<AliZDCMBCalib*>(GetMBCalibData()); 
     TH2F *hZDCvsZEM  = fgMBCalibData->GethZDCvsZEM();
     TH2F *hZDCCvsZEM = fgMBCalibData->GethZDCCvsZEM();
     TH2F *hZDCAvsZEM = fgMBCalibData->GethZDCAvsZEM();
     //
-    TH1D *hNpartDist = fgRecoParam->GethNpartDist();
-    TH1D *hbDist = fgRecoParam->GethbDist();    
-    Float_t clkCenter = fgRecoParam->GetClkCenter();
-    //
     Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
-    Double_t origin = xHighEdge*clkCenter;
+    Double_t origin = xHighEdge*fClkCenter;
     // Ch. debug
     //printf("\n\n  xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
     //
@@ -1190,24 +1268,27 @@ void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
     if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
     if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
     if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;    
-  
+    
+    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);
   //reco->Print("");
   // write the output tree
   clustersTree->Fill();
+  delete reco;
 }
 
 
@@ -1219,10 +1300,8 @@ void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd)
   AliZDCReco reco;
   AliZDCReco* preco = &reco;
   clustersTree->SetBranchAddress("ZDC", &preco);
-
   clustersTree->GetEntry(0);
   //
-  AliESDZDC * esdzdc = esd->GetESDZDC();
   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
   for(Int_t i=0; i<5; i++){
@@ -1237,15 +1316,15 @@ void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd)
      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
   }
   //
-  esdzdc->SetZN1TowerEnergy(tZN1Ene);
-  esdzdc->SetZN2TowerEnergy(tZN2Ene);
-  esdzdc->SetZP1TowerEnergy(tZP1Ene);
-  esdzdc->SetZP2TowerEnergy(tZP2Ene);
+  fESDZDC->SetZN1TowerEnergy(tZN1Ene);
+  fESDZDC->SetZN2TowerEnergy(tZN2Ene);
+  fESDZDC->SetZP1TowerEnergy(tZP1Ene);
+  fESDZDC->SetZP2TowerEnergy(tZP2Ene);
   //
-  esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
-  esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
-  esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
-  esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
+  fESDZDC->SetZN1TowerEnergyLR(tZN1EneLR);
+  fESDZDC->SetZN2TowerEnergyLR(tZN2EneLR);
+  fESDZDC->SetZP1TowerEnergyLR(tZP1EneLR);
+  fESDZDC->SetZP2TowerEnergyLR(tZP2EneLR);
   // 
   Int_t nPart  = reco.GetNParticipants();
   Int_t nPartA = reco.GetNPartSideA();
@@ -1254,23 +1333,41 @@ void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd)
   Double_t bA = reco.GetImpParSideA();
   Double_t bC = reco.GetImpParSideC();
   UInt_t recoFlag = reco.GetRecoFlag();
-  //
-  esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(), 
-             reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(), 
-             nPart, nPartA, nPartC, b, bA, bC, recoFlag);
+  
+  fESDZDC->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), 
+       reco.GetZEM1HRsignal(), reco.GetZEM2HRsignal(), 
+       reco.GetZN2HREnergy(), reco.GetZP2HREnergy(), 
+       nPart, nPartA, nPartC, b, bA, bC, recoFlag);
   
   // Writing ZDC scaler for cross section calculation
   // ONLY IF the scaler has been read during the event
   if(reco.IsScalerOn()==kTRUE){
     UInt_t counts[32];
     for(Int_t jk=0; jk<32; jk++) counts[jk] = reco.GetZDCScaler(jk);
-    esd->SetZDCScaler(counts);
-  }
+    fESDZDC->SetZDCScaler(counts);
+  }    
   
   // Writing TDC data into ZDC ESDs
-  Int_t tdcValues[32];
-  for(Int_t jk=0; jk<32; jk++) tdcValues[jk] = reco.GetZDCTDCData(jk);
-  esd->SetZDCTDC(tdcValues);
+  Int_t tdcValues[32][4]; 
+  Float_t tdcCorrected[32][4];
+  for(Int_t jk=0; jk<32; jk++){
+    for(Int_t lk=0; lk<4; lk++){
+      tdcValues[jk][lk] = reco.GetZDCTDCData(jk, lk);
+    }
+  }
+  // 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;
+    }
+  }
+  fESDZDC->SetZDCTDCData(tdcValues);
+  fESDZDC->SetZDCTDCCorrected(tdcCorrected);
+  fESDZDC->AliESDZDC::SetBit(AliESDZDC::kCorrectedTDCFilled, reco.GetEnergyFlag());
+  fESDZDC->AliESDZDC::SetBit(AliESDZDC::kEnergyCalibratedSignal, kTRUE);
+  
+  if(esd) esd->SetZDCData(fESDZDC);
 }
 
 //_____________________________________________________________________________