Bug fix (Chiara)
[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 "AliESDZDC.h"
31 #include "AliZDCDigit.h"
32 #include "AliZDCRawStream.h"
33 #include "AliZDCReco.h"
34 #include "AliZDCReconstructor.h"
35 #include "AliZDCPedestals.h"
36 #include "AliZDCCalib.h"
37 #include "AliZDCRecParam.h"
38
39
40 ClassImp(AliZDCReconstructor)
41
42
43 //_____________________________________________________________________________
44 AliZDCReconstructor:: AliZDCReconstructor() :
45
46   fZNCen(new TF1("fZNCen", 
47         "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.)),
48   fZNPer(new TF1("fZNPer",
49       "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.)),
50   fZPCen(new TF1("fZPCen",
51        "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.)),
52   fZPPer(new TF1("fZPPer",
53       "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.)),
54   fZDCCen(new TF1("fZDCCen",
55       "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.)),
56   fZDCPer(new TF1("fZDCPer",
57       "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.)),
58   fbCen(new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.)),
59   fbPer(new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.)),
60   //
61   fZEMn(new TF1("fZEMn","121.7-0.1934*x+0.00007565*x*x",0.,1200.)),
62   fZEMp(new TF1("fZEMp","80.05-0.1315*x+0.00005327*x*x",0.,1200.)),
63   fZEMsp(new TF1("fZEMsp","201.7-0.325*x+0.0001292*x*x",0.,1200.)),
64   fZEMb(new TF1("fZEMb",
65         "13.83-0.02851*x+5.101e-5*x*x-7.305e-8*x*x*x+5.101e-11*x*x*x*x-1.25e-14*x*x*x*x*x",0.,1200.)),
66   //
67   fPedData(GetPedData()),
68   fECalibData(GetECalibData()),
69   fRecParam(GetRecParams())
70 {
71   // **** Default constructor
72
73 }
74
75
76 //_____________________________________________________________________________
77 AliZDCReconstructor::~AliZDCReconstructor()
78 {
79 // destructor
80
81   delete fZNCen;
82   delete fZNPer;
83   delete fZPCen;
84   delete fZPPer;
85   delete fZDCCen;
86   delete fZDCPer;
87   delete fbCen;
88   delete fbPer;
89   delete fZEMn;
90   delete fZEMp;
91   delete fZEMsp;
92   delete fZEMb;
93
94 }
95
96
97 //_____________________________________________________________________________
98 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
99 {
100   // *** Local ZDC reconstruction for digits
101   // Works on the current event
102     
103   // Retrieving calibration data  
104   Float_t meanPed[48];
105   for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
106
107   // get digits
108   AliZDCDigit digit;
109   AliZDCDigit* pdigit = &digit;
110   digitsTree->SetBranchAddress("ZDC", &pdigit);
111   //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
112
113   // loop over digits
114   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10]; 
115   Float_t dZEM1Corr[2], dZEM2Corr[2], PMRef1[2], PMRef2[2]; 
116   for(Int_t i=0; i<10; i++){
117      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
118      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = PMRef1[i] = PMRef2[i] = 0.;
119   }  
120   //
121   for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
122    digitsTree->GetEntry(iDigit);
123    if (!pdigit) continue;
124    //  
125    Int_t det = digit.GetSector(0);
126    Int_t quad = digit.GetSector(1);
127    Int_t pedindex = -1, kNch = 24;
128    //printf("\n\t Digit #%d det %d quad %d", iDigit, det, quad);
129    //
130    if(quad != 5){ // ZDC (not reference PTMs!)
131     if(det == 1){ // *** ZNC
132        pedindex = quad;
133        tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
134        if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
135        tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
136        if(tZN1Corr[quad+5]<0.) tZN1Corr[quad] = 0.;
137        //printf("\t pedindex %d tZN1Corr[%d] = %1.0f tZN1Corr[%d] = %1.0f", 
138        //       pedindex, quad, tZN1Corr[quad], quad+5, tZN1Corr[quad+5]);
139     }
140     else if(det == 2){ // *** ZP1
141        pedindex = quad+5;
142        tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
143        if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
144        tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
145        if(tZP1Corr[quad+5]<0.) tZP1Corr[quad] = 0.;
146        //printf("\t pedindex %d tZP1Corr[%d] = %1.0f tZP1Corr[%d] = %1.0f", 
147        //       pedindex, quad, tZP1Corr[quad], quad+5, tZP1Corr[quad+5]);
148     }
149     else if(det == 3){
150        pedindex = quad+9;
151        if(quad == 1){       // *** ZEM1  
152          dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]); 
153          if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
154          dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]); 
155          if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
156          //printf("\t pedindex %d tZEM1Corr[%d] = %1.0f tZEM1Corr[%d] = %1.0f", 
157          //     pedindex, quad, tZEM1Corr[quad], quad+1, tZEM1Corr[quad+1]);
158        }
159        else if(quad == 2){  // *** ZEM2
160          dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]); 
161          if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
162          dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]); 
163          if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
164          //printf("\t pedindex %d tZEM2Corr[%d] = %1.0f tZEM2Corr[%d] = %1.0f", 
165          //     pedindex, quad, tZEM2Corr[quad], quad+1, tZEM2Corr[quad+1]);
166        }
167     }
168     else if(det == 4){  // *** ZN2
169        pedindex = quad+12;
170        tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
171        if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
172        tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
173        if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
174        //printf("\t pedindex %d tZN2Corr[%d] = %1.0f tZN2Corr[%d] = %1.0f\n", 
175        //       pedindex, quad, tZN2Corr[quad], quad+5, tZN2Corr[quad+5]);
176     }
177     else if(det == 5){  // *** ZP2 
178        pedindex = quad+17;
179        tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
180        if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
181        tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
182        if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
183        //printf("\t pedindex %d tZP2Corr[%d] = %1.0f tZP2Corr[%d] = %1.0f\n", 
184        //       pedindex, quad, tZP2Corr[quad], quad+5, tZP2Corr[quad+5]);
185     }
186    }
187    else{ // Reference PMs
188      pedindex = (det-1)/3+22;
189      if(det == 1){
190        PMRef1[0] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
191        if(PMRef1[0]<0.) PMRef1[0] = 0.;
192        PMRef1[1] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
193        if(PMRef2[1]<0.) PMRef1[1] = 0.;
194      }
195      else if(det == 4){
196        PMRef2[0] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
197        if(PMRef2[0]<0.) PMRef2[0] = 0.;
198        PMRef2[1] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+kNch]);
199        if(PMRef2[1]<0.) PMRef2[1] = 0.;
200      }
201    }
202   }
203
204   // reconstruct the event
205     ReconstructEvent(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr, 
206         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
207
208 }
209
210 //_____________________________________________________________________________
211 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
212 {
213   // *** ZDC raw data reconstruction
214   // Works on the current event
215   
216   // Retrieving calibration data  
217   Float_t meanPed[48];
218   for(Int_t jj=0; jj<48; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
219
220   rawReader->Reset();
221
222   // loop over raw data
223   Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10]; 
224   Float_t dZEM1Corr[2], dZEM2Corr[2], PMRef1[2], PMRef2[2]; 
225   for(Int_t i=0; i<10; i++){
226      tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
227      if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = PMRef1[i] = PMRef2[i] = 0.;
228   }  
229   //
230   AliZDCRawStream rawData(rawReader);
231   Int_t kNch = 24;
232   while (rawData.Next()) {
233     if(rawData.IsADCDataWord()){
234      Int_t det = rawData.GetSector(0);
235      Int_t quad = rawData.GetSector(1);
236      Int_t gain = rawData.GetADCGain();
237      Int_t pedindex=0;
238      //
239      if(quad !=5){ // ZDCs (not reference PTMs)
240       if(det == 1){    
241         pedindex = quad;
242         if(gain == 0) tZN1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
243         else tZN1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
244       }
245       else if(det == 2){ 
246         pedindex = quad+5;
247         if(gain == 0) tZP1Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
248         else tZP1Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
249       }
250       else if(det == 3){ 
251         pedindex = quad+9;
252         if(quad==1){     
253           if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
254           else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
255         }
256         else if(quad==2){ 
257           if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
258           else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
259         }
260       }
261       else if(det == 4){       
262         pedindex = quad+12;
263         if(gain == 0) tZN2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
264         else tZN2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
265       }
266       else if(det == 5){
267         pedindex = quad+17;
268         if(gain == 0) tZP2Corr[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
269         else tZP2Corr[quad+5]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]); 
270       }
271       //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n", 
272       //  det,quad,gain, pedindex, meanPed[pedindex]);
273      }
274      else{ // reference PM
275        pedindex = (det-1)/3 + 22;
276        if(det == 1){
277          if(gain==0) PMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
278          else PMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
279        }
280        else if(det ==4){
281          if(gain==0) PMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
282          else PMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
283        }
284      }
285     }//IsADCDataWord
286   }
287     
288   // reconstruct the event
289     ReconstructEvent(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr, 
290         dZEM1Corr, dZEM2Corr, PMRef1, PMRef2);
291
292 }
293
294 //_____________________________________________________________________________
295 void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree, Float_t* ZN1ADCCorr, 
296         Float_t* ZP1ADCCorr, Float_t* ZN2ADCCorr, Float_t* ZP2ADCCorr,
297         Float_t* ZEM1ADCCorr, Float_t* ZEM2ADCCorr, Float_t* PMRef1, Float_t* PMRef2) const
298 {
299   // ***** Reconstruct one event
300   
301   // *** RECONSTRUCTION FROM "REAL" DATA
302   //
303   // Retrieving calibration data
304   // --- Equalization coefficients ---------------------------------------------
305   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
306   for(Int_t ji=0; ji<5; ji++){
307      equalCoeffZN1[ji] = fECalibData->GetZN1EqualCoeff(ji);
308      equalCoeffZP1[ji] = fECalibData->GetZP1EqualCoeff(ji); 
309      equalCoeffZN2[ji] = fECalibData->GetZN2EqualCoeff(ji); 
310      equalCoeffZP2[ji] = fECalibData->GetZP2EqualCoeff(ji); 
311   }
312   // --- Energy calibration factors ------------------------------------
313   Float_t calibEne[4];
314   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fECalibData->GetEnCalib(ij);
315   //
316   // --- Reconstruction parameters ------------------
317   Float_t endPointZEM = fRecParam->GetZEMEndValue();
318   Float_t cutFractionZEM = fRecParam->GetZEMCutFraction();
319   Float_t dZEMSup = fRecParam->GetDZEMSup();
320   Float_t dZEMInf = fRecParam->GetDZEMInf();
321   //
322   Float_t cutValueZEM = endPointZEM*cutFractionZEM;
323   Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
324   Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
325   //
326   Float_t maxValEZN1  = fRecParam->GetEZN1MaxValue();
327   Float_t maxValEZP1  = fRecParam->GetEZP1MaxValue();
328   Float_t maxValEZDC1 = fRecParam->GetEZDC1MaxValue();
329   Float_t maxValEZN2  = fRecParam->GetEZN2MaxValue();
330   Float_t maxValEZP2  = fRecParam->GetEZP2MaxValue();
331   Float_t maxValEZDC2 = fRecParam->GetEZDC2MaxValue();
332   //
333   //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
334   //   " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
335   
336   // Equalization of detector responses
337   Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
338   for(Int_t gi=0; gi<5; gi++){
339      equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
340      equalTowZN1[gi+5] = ZN1ADCCorr[gi+5]*equalCoeffZN1[gi];
341      equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
342      equalTowZP1[gi+5] = ZP1ADCCorr[gi+5]*equalCoeffZP1[gi];
343      equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
344      equalTowZN2[gi+5] = ZN2ADCCorr[gi+5]*equalCoeffZN2[gi];
345      equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
346      equalTowZP2[gi+5] = ZP2ADCCorr[gi+5]*equalCoeffZP2[gi];
347   }
348   
349   // Energy calibration of detector responses
350   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
351   Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
352   for(Int_t gi=0; gi<10; gi++){
353      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
354      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
355      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
356      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
357      //
358      if(gi<5){
359        calibSumZN1[0] += calibTowZN1[gi];
360        calibSumZP1[0] += calibTowZP1[gi];
361        calibSumZN2[0] += calibTowZN2[gi];
362        calibSumZP2[0] += calibTowZP2[gi];
363      }
364      //
365      else{
366        calibSumZN1[1] += calibTowZN1[gi];
367        calibSumZP1[1] += calibTowZP1[gi];
368        calibSumZN2[1] += calibTowZN2[gi];
369        calibSumZP2[1] += calibTowZP2[gi];
370      }
371   }
372   
373   //  ---      Number of detected spectator nucleons
374   //  *** N.B. -> It works only in Pb-Pb
375   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
376   nDetSpecNLeft = (Int_t) (calibSumZN1[0]/2.760);
377   nDetSpecPLeft = (Int_t) (calibSumZP1[0]/2.760);
378   nDetSpecNRight = (Int_t) (calibSumZN2[0]/2.760);
379   nDetSpecPRight = (Int_t) (calibSumZP2[0]/2.760);
380   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
381     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft, 
382     nDetSpecNRight, nDetSpecPRight);*/
383
384   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
385   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
386   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
387   Double_t impPar=0.;
388   //
389   //
390   Float_t corrADCZEMHG = ZEM1ADCCorr[0] + ZEM2ADCCorr[0];
391   //
392   if(corrADCZEMHG > supValueZEM){
393     nGenSpecNLeft  = (Int_t) (fZNCen->Eval(calibSumZN1[0]));
394     nGenSpecPLeft  = (Int_t) (fZPCen->Eval(calibSumZP1[0]));
395     nGenSpecLeft   = (Int_t) (fZDCCen->Eval(calibSumZN1[0]+calibSumZP1[0]));
396     nGenSpecNRight = (Int_t) (fZNCen->Eval(calibSumZN2[0]));
397     nGenSpecPRight = (Int_t) (fZNCen->Eval(calibSumZP2[0]));
398     nGenSpecRight  = (Int_t) (fZNCen->Eval(calibSumZN2[0]+calibSumZP2[0]));
399     impPar  = fbCen->Eval(calibSumZN1[0]+calibSumZP1[0]);
400   }
401   else if(corrADCZEMHG < infValueZEM){
402     nGenSpecNLeft = (Int_t) (fZNPer->Eval(calibSumZN1[0])); 
403     nGenSpecPLeft = (Int_t) (fZPPer->Eval(calibSumZP1[0]));
404     nGenSpecLeft  = (Int_t) (fZDCPer->Eval(calibSumZN1[0]+calibSumZP1[0]));
405     impPar   = fbPer->Eval(calibSumZN1[0]+calibSumZP1[0]);
406   }
407   else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
408     nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
409     nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
410     nGenSpecLeft  = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
411     impPar   =  fZEMb->Eval(corrADCZEMHG);
412   }
413   // 
414   if(calibSumZN1[0]/maxValEZN1>1.)  nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
415   if(calibSumZP1[0]/maxValEZP1>1.)  nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
416   if((calibSumZN1[0]+calibSumZP1[0]/maxValEZDC1)>1.){
417      nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
418      impPar = fZEMb->Eval(corrADCZEMHG);
419   }
420   if(calibSumZN2[0]/maxValEZN2>1.)  nGenSpecNRight = (Int_t) (fZEMn->Eval(corrADCZEMHG));
421   if(calibSumZP2[0]/maxValEZP2>1.)  nGenSpecPRight = (Int_t) (fZEMp->Eval(corrADCZEMHG));
422   if((calibSumZN2[0]+calibSumZP2[0]/maxValEZDC2)>1.) nGenSpecRight = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
423   //
424   if(nGenSpecNLeft>125)    nGenSpecNLeft=125;
425   else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
426   if(nGenSpecPLeft>82)     nGenSpecPLeft=82;
427   else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
428   if(nGenSpecLeft>207)     nGenSpecLeft=207;
429   else if(nGenSpecLeft<0)  nGenSpecLeft=0;
430   
431   //  ---      Number of generated participants (from HIJING parameterization)
432   Int_t nPart, nPartTotLeft, nPartTotRight;
433   nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
434   nPartTotLeft = 207-nGenSpecLeft;
435   nPartTotRight = 207-nGenSpecRight;
436   if(nPart<0) nPart=0;
437   if(nPartTotLeft<0) nPartTotLeft=0;
438   if(nPartTotRight<0) nPartTotRight=0;
439   //
440   // *** DEBUG ***
441   /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
442       "  calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n", 
443       calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
444   printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
445       "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n", 
446       nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, 
447       nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
448   printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
449   */
450   
451   // create the output tree
452   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2, 
453                   calibTowZN1, calibTowZN2, calibTowZP1, calibTowZP2, 
454                   ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
455                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight, 
456                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight, 
457                   nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
458                   
459   AliZDCReco* preco = &reco;
460   const Int_t kBufferSize = 4000;
461   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
462
463   // write the output tree
464   clustersTree->Fill();
465 }
466
467 //_____________________________________________________________________________
468 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
469 {
470   // fill energies and number of participants to the ESD
471
472   AliZDCReco reco;
473   AliZDCReco* preco = &reco;
474   clustersTree->SetBranchAddress("ZDC", &preco);
475
476   clustersTree->GetEntry(0);
477   //
478   AliESDZDC * esdzdc = esd->GetESDZDC();
479   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
480   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
481   for(Int_t i=0; i<5; i++){
482      tZN1Ene[i] = reco.GetZN1HREnTow(i);
483      tZN2Ene[i] = reco.GetZN2HREnTow(i);
484      tZP1Ene[i] = reco.GetZP1HREnTow(i);
485      tZP2Ene[i] = reco.GetZP2HREnTow(i);
486      //
487      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
488      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
489      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
490      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
491   }
492   esdzdc->SetZN1TowerEnergy(tZN1Ene);
493   esdzdc->SetZN2TowerEnergy(tZN2Ene);
494   esdzdc->SetZP1TowerEnergy(tZP1Ene);
495   esdzdc->SetZP2TowerEnergy(tZP2Ene);
496   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
497   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
498   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
499   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
500   // 
501   esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(), 
502               reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(), 
503               reco.GetNPartLeft());
504   //
505   
506 }
507
508 //_____________________________________________________________________________
509 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri) 
510 {
511   // Setting the storage
512
513   Bool_t deleteManager = kFALSE;
514   
515   AliCDBManager *manager = AliCDBManager::Instance();
516   AliCDBStorage *defstorage = manager->GetDefaultStorage();
517   
518   if(!defstorage || !(defstorage->Contains("ZDC"))){ 
519      AliWarning("No default storage set or default storage doesn't contain ZDC!");
520      manager->SetDefaultStorage(uri);
521      deleteManager = kTRUE;
522   }
523  
524   AliCDBStorage *storage = manager->GetDefaultStorage();
525
526   if(deleteManager){
527     AliCDBManager::Instance()->UnsetDefaultStorage();
528     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
529   }
530
531   return storage; 
532 }
533
534 //_____________________________________________________________________________
535 AliZDCPedestals* AliZDCReconstructor::GetPedData() const
536 {
537
538   // Getting pedestal calibration object for ZDC set
539
540   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
541   if(!entry) AliFatal("No calibration data loaded!");  
542
543   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
544   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
545
546   return calibdata;
547 }
548
549 //_____________________________________________________________________________
550 AliZDCCalib* AliZDCReconstructor::GetECalibData() const
551 {
552
553   // Getting energy and equalization calibration object for ZDC set
554
555   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Calib");
556   if(!entry) AliFatal("No calibration data loaded!");  
557
558   AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*>  (entry->GetObject());
559   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
560
561   return calibdata;
562 }
563
564 //_____________________________________________________________________________
565 AliZDCRecParam* AliZDCReconstructor::GetRecParams() const
566 {
567
568   // Getting energy and equalization calibration object for ZDC set
569
570   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecParam");
571   if(!entry) AliFatal("No calibration data loaded!");  
572
573   AliZDCRecParam *calibdata = dynamic_cast<AliZDCRecParam*>  (entry->GetObject());
574   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
575
576   return calibdata;
577 }