Reconstruction classes updated (all cabled channels included in reco object)
[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<10; gi++){
339      equalTowZN1[gi] = ZN1ADCCorr[gi]*equalCoeffZN1[gi];
340      equalTowZP1[gi] = ZP1ADCCorr[gi]*equalCoeffZP1[gi];
341      equalTowZN2[gi] = ZN2ADCCorr[gi]*equalCoeffZN2[gi];
342      equalTowZP2[gi] = ZP2ADCCorr[gi]*equalCoeffZP2[gi];
343   }
344   
345   // Energy calibration of detector responses
346   Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
347   Float_t calibSumZN1[2], calibSumZN2[2], calibSumZP1[2], calibSumZP2[2];
348   for(Int_t gi=0; gi<10; gi++){
349      calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
350      calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
351      calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
352      calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
353      //
354      if(gi<5){
355        calibSumZN1[0] += calibTowZN1[gi];
356        calibSumZP1[0] += calibTowZP1[gi];
357        calibSumZN2[0] += calibTowZN2[gi];
358        calibSumZP2[0] += calibTowZP2[gi];
359      }
360      //
361      else{
362        calibSumZN1[1] += calibTowZN1[gi];
363        calibSumZP1[1] += calibTowZP1[gi];
364        calibSumZN2[1] += calibTowZN2[gi];
365        calibSumZP2[1] += calibTowZP2[gi];
366      }
367   }
368   
369   //  ---      Number of detected spectator nucleons
370   //  *** N.B. -> It works only in Pb-Pb
371   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
372   nDetSpecNLeft = (Int_t) (calibSumZN1[0]/2.760);
373   nDetSpecPLeft = (Int_t) (calibSumZP1[0]/2.760);
374   nDetSpecNRight = (Int_t) (calibSumZN2[0]/2.760);
375   nDetSpecPRight = (Int_t) (calibSumZP2[0]/2.760);
376   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
377     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft, 
378     nDetSpecNRight, nDetSpecPRight);*/
379
380   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
381   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
382   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
383   Double_t impPar=0.;
384   //
385   //
386   Float_t corrADCZEMHG = ZEM1ADCCorr[0] + ZEM2ADCCorr[0];
387   //
388   if(corrADCZEMHG > supValueZEM){
389     nGenSpecNLeft  = (Int_t) (fZNCen->Eval(calibSumZN1[0]));
390     nGenSpecPLeft  = (Int_t) (fZPCen->Eval(calibSumZP1[0]));
391     nGenSpecLeft   = (Int_t) (fZDCCen->Eval(calibSumZN1[0]+calibSumZP1[0]));
392     nGenSpecNRight = (Int_t) (fZNCen->Eval(calibSumZN2[0]));
393     nGenSpecPRight = (Int_t) (fZNCen->Eval(calibSumZP2[0]));
394     nGenSpecRight  = (Int_t) (fZNCen->Eval(calibSumZN2[0]+calibSumZP2[0]));
395     impPar  = fbCen->Eval(calibSumZN1[0]+calibSumZP1[0]);
396   }
397   else if(corrADCZEMHG < infValueZEM){
398     nGenSpecNLeft = (Int_t) (fZNPer->Eval(calibSumZN1[0])); 
399     nGenSpecPLeft = (Int_t) (fZPPer->Eval(calibSumZP1[0]));
400     nGenSpecLeft  = (Int_t) (fZDCPer->Eval(calibSumZN1[0]+calibSumZP1[0]));
401     impPar   = fbPer->Eval(calibSumZN1[0]+calibSumZP1[0]);
402   }
403   else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
404     nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
405     nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
406     nGenSpecLeft  = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
407     impPar   =  fZEMb->Eval(corrADCZEMHG);
408   }
409   // 
410   if(calibSumZN1[0]/maxValEZN1>1.)  nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
411   if(calibSumZP1[0]/maxValEZP1>1.)  nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
412   if((calibSumZN1[0]+calibSumZP1[0]/maxValEZDC1)>1.){
413      nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
414      impPar = fZEMb->Eval(corrADCZEMHG);
415   }
416   if(calibSumZN2[0]/maxValEZN2>1.)  nGenSpecNRight = (Int_t) (fZEMn->Eval(corrADCZEMHG));
417   if(calibSumZP2[0]/maxValEZP2>1.)  nGenSpecPRight = (Int_t) (fZEMp->Eval(corrADCZEMHG));
418   if((calibSumZN2[0]+calibSumZP2[0]/maxValEZDC2)>1.) nGenSpecRight = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
419   //
420   if(nGenSpecNLeft>125)    nGenSpecNLeft=125;
421   else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
422   if(nGenSpecPLeft>82)     nGenSpecPLeft=82;
423   else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
424   if(nGenSpecLeft>207)     nGenSpecLeft=207;
425   else if(nGenSpecLeft<0)  nGenSpecLeft=0;
426   
427   //  ---      Number of generated participants (from HIJING parameterization)
428   Int_t nPart, nPartTotLeft, nPartTotRight;
429   nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
430   nPartTotLeft = 207-nGenSpecLeft;
431   nPartTotRight = 207-nGenSpecRight;
432   if(nPart<0) nPart=0;
433   if(nPartTotLeft<0) nPartTotLeft=0;
434   if(nPartTotRight<0) nPartTotRight=0;
435   //
436   // *** DEBUG ***
437   /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
438       "  calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n", 
439       calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
440   printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
441       "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n", 
442       nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, 
443       nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
444   printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
445   */
446   
447   // create the output tree
448   AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2, 
449                   calibTowZN1, calibTowZN2, calibTowZP1, calibTowZP2, 
450                   ZEM1ADCCorr, ZEM2ADCCorr, PMRef1, PMRef2,
451                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight, 
452                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight, 
453                   nGenSpecPRight, nGenSpecRight, nPartTotLeft, nPartTotRight, impPar);
454                   
455   AliZDCReco* preco = &reco;
456   const Int_t kBufferSize = 4000;
457   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
458
459   // write the output tree
460   clustersTree->Fill();
461 }
462
463 //_____________________________________________________________________________
464 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
465 {
466   // fill energies and number of participants to the ESD
467
468   AliZDCReco reco;
469   AliZDCReco* preco = &reco;
470   clustersTree->SetBranchAddress("ZDC", &preco);
471
472   clustersTree->GetEntry(0);
473   //
474   AliESDZDC * esdzdc = esd->GetESDZDC();
475   Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
476   Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
477   for(Int_t i=0; i<5; i++){
478      tZN1Ene[i] = reco.GetZN1HREnTow(i);
479      tZN2Ene[i] = reco.GetZN2HREnTow(i);
480      tZP1Ene[i] = reco.GetZP1HREnTow(i);
481      tZP2Ene[i] = reco.GetZP2HREnTow(i);
482      //
483      tZN1EneLR[i] = reco.GetZN1LREnTow(i);
484      tZN2EneLR[i] = reco.GetZN2LREnTow(i);
485      tZP1EneLR[i] = reco.GetZP1LREnTow(i);
486      tZP2EneLR[i] = reco.GetZP2LREnTow(i);
487   }
488   esdzdc->SetZN1TowerEnergy(tZN1Ene);
489   esdzdc->SetZN2TowerEnergy(tZN2Ene);
490   esdzdc->SetZP1TowerEnergy(tZP1Ene);
491   esdzdc->SetZP2TowerEnergy(tZP2Ene);
492   esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
493   esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
494   esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
495   esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
496   // 
497   esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(), 
498               reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(), 
499               reco.GetNPartLeft());
500   //
501   
502 }
503
504 //_____________________________________________________________________________
505 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri) 
506 {
507   // Setting the storage
508
509   Bool_t deleteManager = kFALSE;
510   
511   AliCDBManager *manager = AliCDBManager::Instance();
512   AliCDBStorage *defstorage = manager->GetDefaultStorage();
513   
514   if(!defstorage || !(defstorage->Contains("ZDC"))){ 
515      AliWarning("No default storage set or default storage doesn't contain ZDC!");
516      manager->SetDefaultStorage(uri);
517      deleteManager = kTRUE;
518   }
519  
520   AliCDBStorage *storage = manager->GetDefaultStorage();
521
522   if(deleteManager){
523     AliCDBManager::Instance()->UnsetDefaultStorage();
524     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
525   }
526
527   return storage; 
528 }
529
530 //_____________________________________________________________________________
531 AliZDCPedestals* AliZDCReconstructor::GetPedData() const
532 {
533
534   // Getting pedestal calibration object for ZDC set
535
536   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
537   if(!entry) AliFatal("No calibration data loaded!");  
538
539   AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
540   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
541
542   return calibdata;
543 }
544
545 //_____________________________________________________________________________
546 AliZDCCalib* AliZDCReconstructor::GetECalibData() const
547 {
548
549   // Getting energy and equalization calibration object for ZDC set
550
551   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Calib");
552   if(!entry) AliFatal("No calibration data loaded!");  
553
554   AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*>  (entry->GetObject());
555   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
556
557   return calibdata;
558 }
559
560 //_____________________________________________________________________________
561 AliZDCRecParam* AliZDCReconstructor::GetRecParams() const
562 {
563
564   // Getting energy and equalization calibration object for ZDC set
565
566   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecParam");
567   if(!entry) AliFatal("No calibration data loaded!");  
568
569   AliZDCRecParam *calibdata = dynamic_cast<AliZDCRecParam*>  (entry->GetObject());
570   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
571
572   return calibdata;
573 }