New reconstruction algorithm
[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 // class for ZDC reconstruction                                              //
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24
25 #include <TF1.h>
26
27 #include "AliRunLoader.h"
28 #include "AliRawReader.h"
29 #include "AliESDEvent.h"
30 #include "AliZDCDigit.h"
31 #include "AliZDCRawStream.h"
32 #include "AliZDCReco.h"
33 #include "AliZDCReconstructor.h"
34 #include "AliZDCCalibData.h"
35
36
37 ClassImp(AliZDCReconstructor)
38
39
40 //_____________________________________________________________________________
41 AliZDCReconstructor:: AliZDCReconstructor() :
42
43   fZNCen(new TF1("fZNCen", 
44         "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.)),
45   fZNPer(new TF1("fZNPer",
46       "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.)),
47   fZPCen(new TF1("fZPCen",
48        "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.)),
49   fZPPer(new TF1("fZPPer",
50       "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.)),
51   fZDCCen(new TF1("fZDCCen",
52       "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.)),
53   fZDCPer(new TF1("fZDCPer",
54       "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.)),
55   fbCen(new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.)),
56   fbPer(new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.)),
57   fZEMn(new TF1("fZEMn","126.2-0.05399*x+0.000005679*x*x",0.,4000.)),
58   fZEMp(new TF1("fZEMp","82.49-0.03611*x+0.00000385*x*x",0.,4000.)),
59   fZEMsp(new TF1("fZEMsp","208.7-0.09006*x+0.000009526*x*x",0.,4000.)),
60   fZEMb(new TF1("fZEMb",
61         "16.06-0.01633*x+1.44e-5*x*x-6.778e-9*x*x*x+1.438e-12*x*x*x*x-1.112e-16*x*x*x*x*x",0.,4000.)),
62   //
63   fCalibData(GetCalibData())
64
65 {
66   // **** Default constructor
67
68 }
69
70
71 //_____________________________________________________________________________
72 AliZDCReconstructor::~AliZDCReconstructor()
73 {
74 // destructor
75
76   delete fZNCen;
77   delete fZNPer;
78   delete fZPCen;
79   delete fZPPer;
80   delete fZDCCen;
81   delete fZDCPer;
82   delete fbCen;
83   delete fbPer;
84   delete fZEMn;
85   delete fZEMp;
86   delete fZEMsp;
87   delete fZEMb;
88
89 }
90
91
92 //_____________________________________________________________________________
93 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
94 {
95   // *** Local ZDC reconstruction for digits
96   // Works on the current event
97     
98   // Retrieving calibration data  
99   Float_t meanPed[47];
100   for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
101
102   // get digits
103   AliZDCDigit digit;
104   AliZDCDigit* pdigit = &digit;
105   digitsTree->SetBranchAddress("ZDC", &pdigit);
106
107   // loop over digits
108     Float_t ZN1TowCorrHG[5], ZP1TowCorrHG[5], ZEMCorrHG=0., 
109             ZN2TowCorrHG[5], ZP2TowCorrHG[5];
110     Float_t ZN1TowCorrLG[5], ZP1TowCorrLG[5], ZEMCorrLG=0., 
111             ZN2TowCorrLG[5], ZP2TowCorrLG[5];
112     
113   for (Int_t iDigit = 0; iDigit < digitsTree->GetEntries(); iDigit++) {
114     digitsTree->GetEntry(iDigit);
115     if (!pdigit) continue;
116       
117     Int_t det = digit.GetSector(0);
118     Int_t quad = digit.GetSector(1);
119     Int_t pedindex;
120     //
121     if(det == 1){ // *** ZN1
122        pedindex = quad;
123        ZN1TowCorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
124        ZN1TowCorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
125     }
126     else if(det == 2){ // *** ZP1
127        pedindex = quad+10;
128        ZP1TowCorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
129        ZP1TowCorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
130     }
131     else if(det == 3){
132        if(quad == 1){       // *** ZEM1  
133          pedindex = quad+20;
134          ZEMCorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]); 
135          ZEMCorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+2]); 
136        }
137        else if(quad == 2){  // *** ZEM1
138          pedindex = quad+21;
139          ZEMCorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]); 
140          ZEMCorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+2]); 
141        }
142     }
143     else if(det == 4){  // *** ZN2
144        pedindex = quad+24;
145        ZN2TowCorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
146        ZN2TowCorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
147     }
148     else if(det == 5){  // *** ZP2 
149        pedindex = quad+34;
150        ZP2TowCorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
151        ZP2TowCorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
152     }
153   }
154
155   // reconstruct the event
156     ReconstructEvent(clustersTree, ZN1TowCorrHG, ZP1TowCorrHG, ZN2TowCorrHG, 
157         ZP2TowCorrHG, ZN1TowCorrLG, ZP1TowCorrLG, ZN2TowCorrLG, 
158         ZP2TowCorrLG, ZEMCorrHG);
159
160 }
161
162 //_____________________________________________________________________________
163 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
164 {
165   // *** ZDC raw data reconstruction
166   // Works on the current event
167   
168   // Retrieving calibration data  
169   Float_t meanPed[47];
170   for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
171
172   rawReader->Reset();
173
174   // loop over raw data rawDatas
175   Float_t ZN1TowCorrHG[5], ZP1TowCorrHG[5], ZEMCorrHG=0., 
176           ZN2TowCorrHG[5], ZP2TowCorrHG[5];
177   Float_t ZN1TowCorrLG[5], ZP1TowCorrLG[5], ZEMCorrLG=0., 
178           ZN2TowCorrLG[5], ZP2TowCorrLG[5];
179   //
180   AliZDCRawStream rawData(rawReader);
181   while (rawData.Next()) {
182     if(rawData.IsADCDataWord()){
183       Int_t det = rawData.GetSector(0);
184       Int_t quad = rawData.GetSector(1);
185       Int_t gain = rawData.GetADCGain();
186       Int_t pedindex;
187       //
188       if(det == 1){    
189         pedindex = quad;
190         if(gain == 0) ZN1TowCorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
191         else ZN1TowCorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]); 
192       }
193       else if(det == 2){ 
194         pedindex = quad+10;
195         if(gain == 0) ZP1TowCorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
196         else ZP1TowCorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]); 
197       }
198       else if(det == 3){ 
199         if(quad==1){     
200           pedindex = quad+20;
201           if(gain == 0) ZEMCorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
202           else ZEMCorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]); 
203         }
204         else if(quad==2){ 
205           pedindex = rawData.GetSector(1)+21;
206           if(gain == 0) ZEMCorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
207           else ZEMCorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]); 
208         }
209       }
210       else if(det == 4){       
211         pedindex = rawData.GetSector(1)+24;
212         if(gain == 0) ZN2TowCorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
213         else ZN2TowCorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]); 
214       }
215       else if(det == 5){
216         pedindex = rawData.GetSector(1)+34;
217         if(gain == 0) ZP2TowCorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
218         else ZP2TowCorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]); 
219       }
220     }
221   }
222     
223   // reconstruct the event
224     ReconstructEvent(clustersTree, ZN1TowCorrHG, ZP1TowCorrHG, ZN2TowCorrHG, 
225         ZP2TowCorrHG, ZN1TowCorrLG, ZP1TowCorrLG, ZN2TowCorrLG, 
226         ZP2TowCorrLG, ZEMCorrHG);
227
228 }
229
230 //_____________________________________________________________________________
231 void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree, 
232                 Float_t* ZN1ADCCorrHG, Float_t* ZP1ADCCorrHG, 
233                 Float_t* ZN2ADCCorrHG, Float_t* ZP2ADCCorrHG, 
234                 Float_t* ZN1ADCCorrLG, Float_t* ZP1ADCCorrLG, 
235                 Float_t* ZN2ADCCorrLG, Float_t* ZP2ADCCorrLG, 
236                 Float_t ZEMADCCorrHG) const
237 {
238   // ***** Reconstruct one event
239   
240   // *** RECONSTRUCTION FROM SIMULATED DATA
241   // It passes trhough the no. of phe which is known from simulations
242   //  ---      ADCchannel -> photoelectrons
243   // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
244   // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
245   //Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
246   //zn1phe  = ZN1Corr/convFactor;
247   //zp1phe  = ZP1Corr/convFactor;
248   //zemphe = ZEMCorr/convFactor;
249   //zn2phe  = ZN2Corr/convFactor;
250   //zp2phe  = ZP2Corr/convFactor;
251   ////if AliDebug(1,Form("\n    znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
252   //
253   ////  ---      Energy calibration
254   //// Conversion factors for hadronic ZDCs goes from phe yield to TRUE 
255   //// incident energy (conversion from GeV to TeV is included); while for EM 
256   //// calos conversion is from light yield to detected energy calculated by
257   //// GEANT NB -> ZN and ZP conversion factors are constant since incident
258   //// spectators have all the same energy, ZEM energy is obtained through a
259   //// fit over the whole range of incident particle energies 
260   //// (obtained with full HIJING simulations) 
261   //Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
262   //Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
263   //zn1energy  = zn1phe/zn1phexTeV;
264   //zp1energy  = zp1phe/zp1phexTeV;
265   //zdc1energy = zn1energy+zp1energy;
266   //zn2energy  = zn2phe/zn2phexTeV;
267   //zp2energy  = zp2phe/zp2phexTeV;
268   //zdc2energy = zn2energy+zp2energy;
269   //zemenergy = -4.81+0.3238*zemphe;
270   //if(zemenergy<0) zemenergy=0;
271   ////  if AliDebug(1,Form("    znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
272   ////                     "\n          zemenergy = %f TeV\n", znenergy, zpenergy, 
273   ////                     zdcenergy, zemenergy);
274   ////  if(zdcenergy==0)
275   ////    if AliDebug(1,Form("\n\n      ###     ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
276   ////                       " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy); 
277   
278   //
279   // *** RECONSTRUCTION FROM "REAL" DATA
280   //
281   // Retrieving calibration data
282   Float_t ZN1EqualCoeff[5], ZP1EqualCoeff[5], ZN2EqualCoeff[5], ZP2EqualCoeff[5];
283   for(Int_t ji=0; ji<5; ji++){
284      ZN1EqualCoeff[ji] = fCalibData->GetZN1EqualCoeff(ji);
285      ZP1EqualCoeff[ji] = fCalibData->GetZP1EqualCoeff(ji); 
286      ZN2EqualCoeff[ji] = fCalibData->GetZN2EqualCoeff(ji); 
287      ZP2EqualCoeff[ji] = fCalibData->GetZP2EqualCoeff(ji); 
288   }
289   //
290   Float_t CalibEne[4];
291   for(Int_t ij=0; ij<4; ij++) CalibEne[ij] = fCalibData->GetEnCalib(ij);
292   //
293   Float_t ZEMEndPoint = fCalibData->GetZEMEndValue();
294   Float_t ZEMCutFraction = fCalibData->GetZEMCutFraction();
295   Float_t DZEMSup = fCalibData->GetDZEMSup();
296   Float_t DZEMInf = fCalibData->GetDZEMInf();
297   //
298   Float_t ZEMCutValue = ZEMEndPoint*ZEMCutFraction;
299   Float_t ZEMSupValue = ZEMCutValue+(ZEMEndPoint*DZEMSup);
300   Float_t ZEMInfValue = ZEMCutValue-(ZEMEndPoint*DZEMInf);
301   //
302   Float_t EZN1MaxVal = fCalibData->GetEZN1MaxValue();
303   Float_t EZP1MaxVal = fCalibData->GetEZP1MaxValue();
304   Float_t EZDC1MaxVal = fCalibData->GetEZDC1MaxValue();
305   Float_t EZN2MaxVal = fCalibData->GetEZN1MaxValue();
306   Float_t EZP2MaxVal = fCalibData->GetEZP1MaxValue();
307   Float_t EZDC2MaxVal = fCalibData->GetEZDC1MaxValue();
308   
309   // Equalization of detector responses
310   Float_t ZN1EqualTowHG[5], ZN2EqualTowHG[5], ZP1EqualTowHG[5], ZP2EqualTowHG[5];
311   Float_t ZN1EqualTowLG[5], ZN2EqualTowLG[5], ZP1EqualTowLG[5], ZP2EqualTowLG[5];
312   for(Int_t gi=0; gi<5; gi++){
313      ZN1EqualTowHG[gi] = ZN1ADCCorrHG[gi]*ZN1EqualCoeff[gi];
314      ZP1EqualTowHG[gi] = ZP1ADCCorrHG[gi]*ZP1EqualCoeff[gi];
315      ZN2EqualTowHG[gi] = ZN2ADCCorrHG[gi]*ZN2EqualCoeff[gi];
316      ZP2EqualTowHG[gi] = ZP2ADCCorrHG[gi]*ZP2EqualCoeff[gi];
317      //
318      ZN1EqualTowLG[gi] = ZN1ADCCorrLG[gi]*ZN1EqualCoeff[gi];
319      ZP1EqualTowLG[gi] = ZP1ADCCorrLG[gi]*ZP1EqualCoeff[gi];
320      ZN2EqualTowLG[gi] = ZN2ADCCorrLG[gi]*ZN2EqualCoeff[gi];
321      ZP2EqualTowLG[gi] = ZP2ADCCorrLG[gi]*ZP2EqualCoeff[gi];
322   }
323   
324   // Energy calibration of detector responses
325   Float_t ZN1CalibTowHG[5], ZN2CalibTowHG[5], ZP1CalibTowHG[5], ZP2CalibTowHG[5];
326   Float_t ZN1CalibSumHG=0., ZN2CalibSumHG=0., ZP1CalibSumHG=0., ZP2CalibSumHG=0.;
327   Float_t ZN1CalibTowLG[5], ZN2CalibTowLG[5], ZP1CalibTowLG[5], ZP2CalibTowLG[5];
328   Float_t ZN1CalibSumLG=0., ZN2CalibSumLG=0., ZP1CalibSumLG=0., ZP2CalibSumLG=0.;
329   for(Int_t gi=0; gi<5; gi++){
330      ZN1CalibTowHG[gi] = ZN1EqualTowHG[gi]*CalibEne[0];
331      ZP1CalibTowHG[gi] = ZP1EqualTowHG[gi]*CalibEne[1];
332      ZN2CalibTowHG[gi] = ZN2EqualTowHG[gi]*CalibEne[2];
333      ZP2CalibTowHG[gi] = ZP2EqualTowHG[gi]*CalibEne[3];
334      ZN1CalibSumHG += ZN1CalibTowHG[gi];
335      ZP1CalibSumHG += ZP1CalibTowHG[gi];
336      ZN2CalibSumHG += ZN2CalibTowHG[gi];
337      ZP2CalibSumHG += ZP2CalibTowHG[gi];
338      //
339      ZN1CalibTowLG[gi] = ZN1EqualTowLG[gi]*CalibEne[0];
340      ZP1CalibTowLG[gi] = ZP1EqualTowLG[gi]*CalibEne[1];
341      ZN2CalibTowLG[gi] = ZN2EqualTowLG[gi]*CalibEne[2];
342      ZP2CalibTowLG[gi] = ZP2EqualTowLG[gi]*CalibEne[3];
343      ZN1CalibSumLG += ZN1CalibTowLG[gi];
344      ZP1CalibSumLG += ZP1CalibTowLG[gi];
345      ZN2CalibSumLG += ZN2CalibTowLG[gi];
346      ZP2CalibSumLG += ZP2CalibTowLG[gi];
347   }
348   
349   //  ---      Number of detected spectator nucleons
350   //  *** N.B. -> It works only in Pb-Pb
351   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
352   nDetSpecNLeft = (Int_t) (ZN1CalibSumHG/2.760);
353   nDetSpecPLeft = (Int_t) (ZP1CalibSumHG/2.760);
354   nDetSpecNRight = (Int_t) (ZN2CalibSumHG/2.760);
355   nDetSpecPRight = (Int_t) (ZP2CalibSumHG/2.760);
356
357   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
358   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
359   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
360   Double_t impPar=0.;
361   //
362   // *** RECONSTRUCTION FROM SIMULATED DATA
363   // Cut value for Ezem (GeV)
364   // ### Results from production  -> 0<b<18 fm (Apr 2002)
365   /*Float_t eZEMCut = 420.;
366   Float_t deltaEZEMSup = 690.; 
367   Float_t deltaEZEMInf = 270.; 
368   if(zemenergy > (eZEMCut+deltaEZEMSup)){
369     nGenSpecNLeft  = (Int_t) (fZNCen->Eval(ZN1CalibSum));
370     nGenSpecPLeft  = (Int_t) (fZPCen->Eval(ZP1CalibSum));
371     nGenSpecLeft   = (Int_t) (fZDCCen->Eval(ZN1CalibSum+ZP1CalibSum));
372     nGenSpecNRight = (Int_t) (fZNCen->Eval(ZN2CalibSum));
373     nGenSpecPRight = (Int_t) (fZNCen->Eval(ZP2CalibSum));
374     nGenSpecRight  = (Int_t) (fZNCen->Eval(ZN2CalibSum+ZP2CalibSum));
375     impPar  = fbCen->Eval(ZN1CalibSum+ZP1CalibSum);
376   }
377   else if(zemenergy < (eZEMCut-deltaEZEMInf)){
378     nGenSpecNLeft = (Int_t) (fZNPer->Eval(ZN1CalibSum)); 
379     nGenSpecPLeft = (Int_t) (fZPPer->Eval(ZP1CalibSum));
380     nGenSpecLeft  = (Int_t) (fZDCPer->Eval(ZN1CalibSum+ZP1CalibSum));
381     impPar   = fbPer->Eval(ZN1CalibSum+ZP1CalibSum);
382   }
383   else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
384     nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
385     nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
386     nGenSpecLeft  = (Int_t)(fZEMsp->Eval(zemenergy));
387     impPar   =  fZEMb->Eval(zemenergy);
388   }
389   // ### Results from production  -> 0<b<18 fm (Apr 2002)
390   if(ZN1CalibSum>162.)  nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
391   if(ZP1CalibSum>59.75)  nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
392   if(ZN1CalibSum+ZP1CalibSum>221.5) nGenSpecLeft  = (Int_t)(fZEMsp->Eval(zemenergy));
393   if(ZN1CalibSum+ZP1CalibSum>220.)  impPar    =  fZEMb->Eval(zemenergy);
394   */
395   //
396   //
397   // *** RECONSTRUCTION FROM REAL DATA
398   //
399   if(ZEMADCCorrHG > ZEMSupValue){
400     nGenSpecNLeft  = (Int_t) (fZNCen->Eval(ZN1CalibSumHG));
401     nGenSpecPLeft  = (Int_t) (fZPCen->Eval(ZP1CalibSumHG));
402     nGenSpecLeft   = (Int_t) (fZDCCen->Eval(ZN1CalibSumHG+ZP1CalibSumHG));
403     nGenSpecNRight = (Int_t) (fZNCen->Eval(ZN2CalibSumHG));
404     nGenSpecPRight = (Int_t) (fZNCen->Eval(ZP2CalibSumHG));
405     nGenSpecRight  = (Int_t) (fZNCen->Eval(ZN2CalibSumHG+ZP2CalibSumHG));
406     impPar  = fbCen->Eval(ZN1CalibSumHG+ZP1CalibSumHG);
407   }
408   else if(ZEMADCCorrHG < ZEMInfValue){
409     nGenSpecNLeft = (Int_t) (fZNPer->Eval(ZN1CalibSumHG)); 
410     nGenSpecPLeft = (Int_t) (fZPPer->Eval(ZP1CalibSumHG));
411     nGenSpecLeft  = (Int_t) (fZDCPer->Eval(ZN1CalibSumHG+ZP1CalibSumHG));
412     impPar   = fbPer->Eval(ZN1CalibSumHG+ZP1CalibSumHG);
413   }
414   else if(ZEMADCCorrHG >= ZEMInfValue && ZEMADCCorrHG <= ZEMSupValue){
415     nGenSpecNLeft = (Int_t) (fZEMn->Eval(ZEMADCCorrHG));
416     nGenSpecPLeft = (Int_t) (fZEMp->Eval(ZEMADCCorrHG));
417     nGenSpecLeft  = (Int_t)(fZEMsp->Eval(ZEMADCCorrHG));
418     impPar   =  fZEMb->Eval(ZEMADCCorrHG);
419   }
420   // 
421   if(ZN1CalibSumHG/EZN1MaxVal>1.)  nGenSpecNLeft = (Int_t) (fZEMn->Eval(ZEMADCCorrHG));
422   if(ZP1CalibSumHG/EZP1MaxVal>1.)  nGenSpecPLeft = (Int_t) (fZEMp->Eval(ZEMADCCorrHG));
423   if((ZN1CalibSumHG+ZP1CalibSumHG/EZDC1MaxVal)>1.){
424      nGenSpecLeft = (Int_t)(fZEMsp->Eval(ZEMADCCorrHG));
425      impPar = fZEMb->Eval(ZEMADCCorrHG);
426   }
427   if(ZN2CalibSumHG/EZN2MaxVal>1.)  nGenSpecNRight = (Int_t) (fZEMn->Eval(ZEMADCCorrHG));
428   if(ZP2CalibSumHG/EZP2MaxVal>1.)  nGenSpecPRight = (Int_t) (fZEMp->Eval(ZEMADCCorrHG));
429   if((ZN2CalibSumHG+ZP2CalibSumHG/EZDC2MaxVal)>1.) nGenSpecRight = (Int_t)(fZEMsp->Eval(ZEMADCCorrHG));
430   //
431   if(nGenSpecNLeft>125)    nGenSpecNLeft=125;
432   else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
433   if(nGenSpecPLeft>82)     nGenSpecPLeft=82;
434   else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
435   if(nGenSpecLeft>207)     nGenSpecLeft=207;
436   else if(nGenSpecLeft<0)  nGenSpecLeft=0;
437   
438   //  ---      Number of generated participants (from HIJING parameterization)
439   Int_t nPart, nPartTotLeft, nPartTotRight;
440   nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
441   nPartTotLeft = 207-nGenSpecLeft;
442   nPartTotRight = 207-nGenSpecRight;
443
444   // create the output tree
445   AliZDCReco reco(ZN1CalibSumHG, ZP1CalibSumHG, ZN2CalibSumHG, ZP2CalibSumHG, 
446                   ZN1CalibTowLG, ZN2CalibTowLG, ZP1CalibTowLG, ZP2CalibTowLG, 
447                   ZEMADCCorrHG, 
448                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight, 
449                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight, 
450                   nGenSpecPRight, nGenSpecRight,
451                   nPartTotLeft, nPartTotRight, impPar);
452                   
453   AliZDCReco* preco = &reco;
454   const Int_t kBufferSize = 4000;
455   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
456
457   // write the output tree
458   clustersTree->Fill();
459 }
460
461 //_____________________________________________________________________________
462 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
463 {
464   // fill energies and number of participants to the ESD
465
466   AliZDCReco reco;
467   AliZDCReco* preco = &reco;
468   clustersTree->SetBranchAddress("ZDC", &preco);
469
470   clustersTree->GetEntry(0);
471   esd->SetZDC(reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEMsignal(),
472               reco.GetZN2Energy(), reco.GetZP2Energy(), 
473               reco.GetNPartLeft());
474 }
475
476 //_____________________________________________________________________________
477 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri) 
478 {
479   // Setting the storage
480
481   Bool_t deleteManager = kFALSE;
482   
483   AliCDBManager *manager = AliCDBManager::Instance();
484   AliCDBStorage *defstorage = manager->GetDefaultStorage();
485   
486   if(!defstorage || !(defstorage->Contains("ZDC"))){ 
487      AliWarning("No default storage set or default storage doesn't contain ZDC!");
488      manager->SetDefaultStorage(uri);
489      deleteManager = kTRUE;
490   }
491  
492   AliCDBStorage *storage = manager->GetDefaultStorage();
493
494   if(deleteManager){
495     AliCDBManager::Instance()->UnsetDefaultStorage();
496     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
497   }
498
499   return storage; 
500 }
501
502 //_____________________________________________________________________________
503 AliZDCCalibData* AliZDCReconstructor::GetCalibData() const
504 {
505
506   // Getting calibration object for ZDC set
507
508   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Data");
509   if(!entry) AliFatal("No calibration data loaded!");  
510
511   AliZDCCalibData *calibdata = dynamic_cast<AliZDCCalibData*>  (entry->GetObject());
512   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
513
514   return calibdata;
515 }