ZEM signal not already splitted in ZDC ESD
[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   //
58   fZEMn(new TF1("fZEMn","121.7-0.1934*x+0.00007565*x*x",0.,1200.)),
59   fZEMp(new TF1("fZEMp","80.05-0.1315*x+0.00005327*x*x",0.,1200.)),
60   fZEMsp(new TF1("fZEMsp","201.7-0.325*x+0.0001292*x*x",0.,1200.)),
61   fZEMb(new TF1("fZEMb",
62         "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.)),
63   //
64   fCalibData(GetCalibData())
65
66 {
67   // **** Default constructor
68
69 }
70
71
72 //_____________________________________________________________________________
73 AliZDCReconstructor::~AliZDCReconstructor()
74 {
75 // destructor
76
77   delete fZNCen;
78   delete fZNPer;
79   delete fZPCen;
80   delete fZPPer;
81   delete fZDCCen;
82   delete fZDCPer;
83   delete fbCen;
84   delete fbPer;
85   delete fZEMn;
86   delete fZEMp;
87   delete fZEMsp;
88   delete fZEMb;
89
90 }
91
92
93 //_____________________________________________________________________________
94 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
95 {
96   // *** Local ZDC reconstruction for digits
97   // Works on the current event
98     
99   // Retrieving calibration data  
100   Float_t meanPed[47];
101   for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
102
103   // get digits
104   AliZDCDigit digit;
105   AliZDCDigit* pdigit = &digit;
106   digitsTree->SetBranchAddress("ZDC", &pdigit);
107
108   // loop over digits
109   Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.}; 
110   Float_t dZEM1CorrHG=0., dZEM2CorrHG=0.; 
111   Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
112   Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
113   Float_t dZEM1CorrLG=0., dZEM2CorrLG=0.; 
114   Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
115   
116   //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
117   for (Int_t iDigit = 0; iDigit < (digitsTree->GetEntries()/2); iDigit++) {
118     digitsTree->GetEntry(iDigit);
119     if (!pdigit) continue;
120     //pdigit->Print("");
121     //  
122     Int_t det = digit.GetSector(0);
123     Int_t quad = digit.GetSector(1);
124     Int_t pedindex = -1;
125     //printf("\n\t #%d det %d quad %d", iDigit, det, quad);
126     //
127     if(det == 1){ // *** ZN1
128        pedindex = quad;
129        tZN1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
130        if(tZN1CorrHG[quad]<0.) tZN1CorrHG[quad] = 0.;
131        tZN1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
132        if(tZN1CorrLG[quad]<0.) tZN1CorrLG[quad] = 0.;
133        //printf("\t pedindex %d tZN1CorrHG[%d] = %1.0f tZN1CorrLG[%d] = %1.0f", 
134        //       pedindex, quad, tZN1CorrHG[quad], quad, tZN1CorrLG[quad]);
135     }
136     else if(det == 2){ // *** ZP1
137        pedindex = quad+10;
138        tZP1CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
139        if(tZP1CorrLG[quad]<0.) tZP1CorrLG[quad] = 0.;
140        tZP1CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
141        if(tZP1CorrHG[quad]<0.) tZP1CorrHG[quad] = 0.;
142        //printf("\t pedindex %d tZP1CorrHG[%d] = %1.0f tZP1CorrLG[%d] = %1.0f", 
143        //       pedindex, quad, tZP1CorrHG[quad], quad, tZP1CorrLG[quad]);
144     }
145     else if(det == 3){
146        if(quad == 1){       // *** ZEM1  
147          pedindex = quad+19;
148          dZEM1CorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]); 
149          if(dZEM1CorrHG<0.) dZEM1CorrHG = 0.;
150          dZEM1CorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+2]); 
151          if(dZEM1CorrLG<0.) dZEM1CorrLG = 0.;
152          //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n", 
153          //     pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEM1CorrHG);
154          //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n", 
155          //     pedindex+2, digit.GetADCValue(1), meanPed[pedindex+2], dZEM1CorrLG);
156          ////printf("\t pedindex %d dZEMCorrHG = %1.0f dZEMCorrLG = %1.0f\n", pedindex, dZEM1CorrHG, dZEM1CorrLG);
157        }
158        else if(quad == 2){  // *** ZEM2
159          pedindex = quad+19;
160          dZEM2CorrHG += (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]); 
161          if(dZEM2CorrHG<0.) dZEM2CorrHG = 0.;
162          dZEM2CorrLG += (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+2]); 
163          if(dZEM2CorrLG<0.) dZEM2CorrLG = 0.;
164          //printf("\t pedindex %d ADC(0) = %d ped = %1.0f ADCCorr = %1.0f\n", 
165          //     pedindex, digit.GetADCValue(0), meanPed[pedindex], dZEM2CorrHG);
166          //printf("\t pedindex %d ADC(1) = %d ped = %1.0f ADCCorr = %1.0f\n", 
167          //     pedindex+2, digit.GetADCValue(1),meanPed[pedindex+2], dZEM2CorrLG);
168          ////printf("\t pedindex %d dZEMCorrHG = %1.0f dZEMCorrLG = %1.0f\n", pedindex, dZEM2CorrHG, dZEM2CorrLG);
169        }
170     }
171     else if(det == 4){  // *** ZN2
172        pedindex = quad+24;
173        tZN2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
174        if(tZN2CorrHG[quad]<0.) tZN2CorrHG[quad] = 0.;
175        tZN2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
176        if(tZN2CorrLG[quad]<0.) tZN2CorrLG[quad] = 0.;
177        //printf("\t pedindex %d tZN2CorrHG[%d] = %1.0f tZN2CorrLG[%d] = %1.0f\n", 
178        //       pedindex, quad, tZN2CorrHG[quad], quad, tZN2CorrLG[quad]);
179     }
180     else if(det == 5){  // *** ZP2 
181        pedindex = quad+34;
182        tZP2CorrHG[quad] = (Float_t) (digit.GetADCValue(0)-meanPed[pedindex]);
183        if(tZP2CorrHG[quad]<0.) tZP2CorrHG[quad] = 0.;
184        tZP2CorrLG[quad] = (Float_t) (digit.GetADCValue(1)-meanPed[pedindex+5]);
185        if(tZP2CorrLG[quad]<0.) tZP2CorrLG[quad] = 0.;
186        //printf("\t pedindex %d tZP2CorrHG[%d] = %1.0f tZP2CorrLG[%d] = %1.0f\n", 
187        //       pedindex, quad, tZP2CorrHG[quad], quad, tZP2CorrLG[quad]);
188     }
189   }
190
191   // reconstruct the event
192     ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG, 
193         tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG, 
194         tZP2CorrLG, dZEM1CorrHG, dZEM2CorrHG);
195
196 }
197
198 //_____________________________________________________________________________
199 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
200 {
201   // *** ZDC raw data reconstruction
202   // Works on the current event
203   
204   // Retrieving calibration data  
205   Float_t meanPed[47];
206   for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
207
208   rawReader->Reset();
209
210   // loop over raw data rawDatas
211   Float_t tZN1CorrHG[]={0.,0.,0.,0.,0.}, tZP1CorrHG[]={0.,0.,0.,0.,0.};
212   Float_t dZEM1CorrHG=0., dZEM2CorrHG=0.;
213   Float_t tZN2CorrHG[]={0.,0.,0.,0.,0.}, tZP2CorrHG[]={0.,0.,0.,0.,0.};
214   Float_t tZN1CorrLG[]={0.,0.,0.,0.,0.}, tZP1CorrLG[]={0.,0.,0.,0.,0.};
215   Float_t dZEM1CorrLG=0., dZEM2CorrLG=0.; 
216   Float_t tZN2CorrLG[]={0.,0.,0.,0.,0.}, tZP2CorrLG[]={0.,0.,0.,0.,0.};
217   //
218   AliZDCRawStream rawData(rawReader);
219   while (rawData.Next()) {
220     if(rawData.IsADCDataWord()){
221       Int_t det = rawData.GetSector(0);
222       Int_t quad = rawData.GetSector(1);
223       Int_t gain = rawData.GetADCGain();
224       Int_t pedindex;
225       //
226       if(det == 1){    
227         pedindex = quad;
228         if(gain == 0) tZN1CorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
229         else tZN1CorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]); 
230       }
231       else if(det == 2){ 
232         pedindex = quad+10;
233         if(gain == 0) tZP1CorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
234         else tZP1CorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]); 
235       }
236       else if(det == 3){ 
237         if(quad==1){     
238           pedindex = quad+20;
239           if(gain == 0) dZEM1CorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
240           else dZEM1CorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]); 
241         }
242         else if(quad==2){ 
243           pedindex = rawData.GetSector(1)+21;
244           if(gain == 0) dZEM2CorrHG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
245           else dZEM2CorrLG += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]); 
246         }
247       }
248       else if(det == 4){       
249         pedindex = rawData.GetSector(1)+24;
250         if(gain == 0) tZN2CorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
251         else tZN2CorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+2]); 
252       }
253       else if(det == 5){
254         pedindex = rawData.GetSector(1)+34;
255         if(gain == 0) tZP2CorrHG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]); 
256         else tZP2CorrLG[quad]  += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+5]); 
257       }
258     }
259   }
260     
261   // reconstruct the event
262     ReconstructEvent(clustersTree, tZN1CorrHG, tZP1CorrHG, tZN2CorrHG, 
263         tZP2CorrHG, tZN1CorrLG, tZP1CorrLG, tZN2CorrLG, 
264         tZP2CorrLG, dZEM1CorrHG, dZEM2CorrHG);
265
266 }
267
268 //_____________________________________________________________________________
269 void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree, 
270                 Float_t* ZN1ADCCorrHG, Float_t* ZP1ADCCorrHG, 
271                 Float_t* ZN2ADCCorrHG, Float_t* ZP2ADCCorrHG, 
272                 Float_t* ZN1ADCCorrLG, Float_t* ZP1ADCCorrLG, 
273                 Float_t* ZN2ADCCorrLG, Float_t* ZP2ADCCorrLG, 
274                 Float_t corrADCZEM1HG, Float_t corrADCZEM2HG) const
275 {
276   // ***** Reconstruct one event
277   
278   // *** RECONSTRUCTION FROM SIMULATED DATA
279   // It passes trhough the no. of phe which is known from simulations
280   //  ---      ADCchannel -> photoelectrons
281   // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
282   // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
283   //Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
284   //zn1phe  = ZN1Corr/convFactor;
285   //zp1phe  = ZP1Corr/convFactor;
286   //zemphe = ZEMCorr/convFactor;
287   //zn2phe  = ZN2Corr/convFactor;
288   //zp2phe  = ZP2Corr/convFactor;
289   ////if AliDebug(1,Form("\n    znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
290   //
291   ////  ---      Energy calibration
292   //// Conversion factors for hadronic ZDCs goes from phe yield to TRUE 
293   //// incident energy (conversion from GeV to TeV is included); while for EM 
294   //// calos conversion is from light yield to detected energy calculated by
295   //// GEANT NB -> ZN and ZP conversion factors are constant since incident
296   //// spectators have all the same energy, ZEM energy is obtained through a
297   //// fit over the whole range of incident particle energies 
298   //// (obtained with full HIJING simulations) 
299   //Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
300   //Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
301   //zn1energy  = zn1phe/zn1phexTeV;
302   //zp1energy  = zp1phe/zp1phexTeV;
303   //zdc1energy = zn1energy+zp1energy;
304   //zn2energy  = zn2phe/zn2phexTeV;
305   //zp2energy  = zp2phe/zp2phexTeV;
306   //zdc2energy = zn2energy+zp2energy;
307   //zemenergy = -4.81+0.3238*zemphe;
308   //if(zemenergy<0) zemenergy=0;
309   ////  if AliDebug(1,Form("    znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
310   ////                     "\n          zemenergy = %f TeV\n", znenergy, zpenergy, 
311   ////                     zdcenergy, zemenergy);
312   ////  if(zdcenergy==0)
313   ////    if AliDebug(1,Form("\n\n      ###     ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
314   ////                       " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy); 
315   
316   //
317   // *** RECONSTRUCTION FROM "REAL" DATA
318   //
319   // Retrieving calibration data
320   Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
321   for(Int_t ji=0; ji<5; ji++){
322      equalCoeffZN1[ji] = fCalibData->GetZN1EqualCoeff(ji);
323      equalCoeffZP1[ji] = fCalibData->GetZP1EqualCoeff(ji); 
324      equalCoeffZN2[ji] = fCalibData->GetZN2EqualCoeff(ji); 
325      equalCoeffZP2[ji] = fCalibData->GetZP2EqualCoeff(ji); 
326   }
327   //
328   Float_t calibEne[4];
329   for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fCalibData->GetEnCalib(ij);
330   //
331   Float_t endPointZEM = fCalibData->GetZEMEndValue();
332   Float_t cutFractionZEM = fCalibData->GetZEMCutFraction();
333   Float_t dZEMSup = fCalibData->GetDZEMSup();
334   Float_t dZEMInf = fCalibData->GetDZEMInf();
335   //
336   Float_t cutValueZEM = endPointZEM*cutFractionZEM;
337   Float_t supValueZEM = cutValueZEM+(endPointZEM*dZEMSup);
338   Float_t infValueZEM = cutValueZEM-(endPointZEM*dZEMInf);
339   //
340   Float_t maxValEZN1 = fCalibData->GetEZN1MaxValue();
341   Float_t maxValEZP1 = fCalibData->GetEZP1MaxValue();
342   Float_t maxValEZDC1 = fCalibData->GetEZDC1MaxValue();
343   Float_t maxValEZN2 = fCalibData->GetEZN2MaxValue();
344   Float_t maxValEZP2 = fCalibData->GetEZP2MaxValue();
345   Float_t maxValEZDC2 = fCalibData->GetEZDC2MaxValue();
346   //
347   //printf("\n\t AliZDCReconstructor -> ZEMEndPoint %1.0f, ZEMCutValue %1.0f,"
348   //   " ZEMSupValue %1.0f, ZEMInfValue %1.0f\n",endPointZEM,cutValueZEM,supValueZEM,infValueZEM);
349   
350   // Equalization of detector responses
351   Float_t equalTowZN1HG[5], equalTowZN2HG[5], equalTowZP1HG[5], equalTowZP2HG[5];
352   Float_t equalTowZN1LG[5], equalTowZN2LG[5], equalTowZP1LG[5], equalTowZP2LG[5];
353   for(Int_t gi=0; gi<5; gi++){
354      equalTowZN1HG[gi] = ZN1ADCCorrHG[gi]*equalCoeffZN1[gi];
355      equalTowZP1HG[gi] = ZP1ADCCorrHG[gi]*equalCoeffZP1[gi];
356      equalTowZN2HG[gi] = ZN2ADCCorrHG[gi]*equalCoeffZN2[gi];
357      equalTowZP2HG[gi] = ZP2ADCCorrHG[gi]*equalCoeffZP2[gi];
358      //
359      equalTowZN1LG[gi] = ZN1ADCCorrLG[gi]*equalCoeffZN1[gi];
360      equalTowZP1LG[gi] = ZP1ADCCorrLG[gi]*equalCoeffZP1[gi];
361      equalTowZN2LG[gi] = ZN2ADCCorrLG[gi]*equalCoeffZN2[gi];
362      equalTowZP2LG[gi] = ZP2ADCCorrLG[gi]*equalCoeffZP2[gi];
363   }
364   
365   // Energy calibration of detector responses
366   Float_t calibTowZN1HG[5], calibTowZN2HG[5], calibTowZP1HG[5], calibTowZP2HG[5];
367   Float_t calibSumZN1HG=0., calibSumZN2HG=0., calibSumZP1HG=0., calibSumZP2HG=0.;
368   Float_t calibTowZN1LG[5], calibTowZN2LG[5], calibTowZP1LG[5], calibTowZP2LG[5];
369   Float_t calibSumZN1LG=0., calibSumZN2LG=0., calibSumZ12LG=0., calibSumZP2LG=0.;
370   for(Int_t gi=0; gi<5; gi++){
371      calibTowZN1HG[gi] = equalTowZN1HG[gi]*calibEne[0];
372      calibTowZP1HG[gi] = equalTowZP1HG[gi]*calibEne[1];
373      calibTowZN2HG[gi] = equalTowZN2HG[gi]*calibEne[2];
374      calibTowZP2HG[gi] = equalTowZP2HG[gi]*calibEne[3];
375      calibSumZN1HG += calibTowZN1HG[gi];
376      calibSumZP1HG += calibTowZP1HG[gi];
377      calibSumZN2HG += calibTowZN2HG[gi];
378      calibSumZP2HG += calibTowZP2HG[gi];
379      //
380      calibTowZN1LG[gi] = equalTowZN1LG[gi]*calibEne[0];
381      calibTowZP1LG[gi] = equalTowZP1LG[gi]*calibEne[1];
382      calibTowZN2LG[gi] = equalTowZN2LG[gi]*calibEne[2];
383      calibTowZP2LG[gi] = equalTowZP2LG[gi]*calibEne[3];
384      calibSumZN1LG += calibTowZN1LG[gi];
385      calibSumZ12LG += calibTowZP1LG[gi];
386      calibSumZN2LG += calibTowZN2LG[gi];
387      calibSumZP2LG += calibTowZP2LG[gi];
388   }
389   
390   //  ---      Number of detected spectator nucleons
391   //  *** N.B. -> It works only in Pb-Pb
392   Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
393   nDetSpecNLeft = (Int_t) (calibSumZN1HG/2.760);
394   nDetSpecPLeft = (Int_t) (calibSumZP1HG/2.760);
395   nDetSpecNRight = (Int_t) (calibSumZN2HG/2.760);
396   nDetSpecPRight = (Int_t) (calibSumZP2HG/2.760);
397   /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
398     " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft, 
399     nDetSpecNRight, nDetSpecPRight);*/
400
401   //  ---      Number of generated spectator nucleons (from HIJING parameterization)
402   Int_t nGenSpecNLeft=0, nGenSpecPLeft=0, nGenSpecLeft=0;
403   Int_t nGenSpecNRight=0, nGenSpecPRight=0, nGenSpecRight=0;
404   Double_t impPar=0.;
405   //
406   // *** RECONSTRUCTION FROM SIMULATED DATA
407   // Cut value for Ezem (GeV)
408   // ### Results from production  -> 0<b<18 fm (Apr 2002)
409   /*Float_t eZEMCut = 420.;
410   Float_t deltaEZEMSup = 690.; 
411   Float_t deltaEZEMInf = 270.; 
412   if(zemenergy > (eZEMCut+deltaEZEMSup)){
413     nGenSpecNLeft  = (Int_t) (fZNCen->Eval(ZN1CalibSum));
414     nGenSpecPLeft  = (Int_t) (fZPCen->Eval(ZP1CalibSum));
415     nGenSpecLeft   = (Int_t) (fZDCCen->Eval(ZN1CalibSum+ZP1CalibSum));
416     nGenSpecNRight = (Int_t) (fZNCen->Eval(ZN2CalibSum));
417     nGenSpecPRight = (Int_t) (fZNCen->Eval(ZP2CalibSum));
418     nGenSpecRight  = (Int_t) (fZNCen->Eval(ZN2CalibSum+ZP2CalibSum));
419     impPar  = fbCen->Eval(ZN1CalibSum+ZP1CalibSum);
420   }
421   else if(zemenergy < (eZEMCut-deltaEZEMInf)){
422     nGenSpecNLeft = (Int_t) (fZNPer->Eval(ZN1CalibSum)); 
423     nGenSpecPLeft = (Int_t) (fZPPer->Eval(ZP1CalibSum));
424     nGenSpecLeft  = (Int_t) (fZDCPer->Eval(ZN1CalibSum+ZP1CalibSum));
425     impPar   = fbPer->Eval(ZN1CalibSum+ZP1CalibSum);
426   }
427   else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
428     nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
429     nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
430     nGenSpecLeft  = (Int_t)(fZEMsp->Eval(zemenergy));
431     impPar   =  fZEMb->Eval(zemenergy);
432   }
433   // ### Results from production  -> 0<b<18 fm (Apr 2002)
434   if(ZN1CalibSum>162.)  nGenSpecNLeft = (Int_t) (fZEMn->Eval(zemenergy));
435   if(ZP1CalibSum>59.75)  nGenSpecPLeft = (Int_t) (fZEMp->Eval(zemenergy));
436   if(ZN1CalibSum+ZP1CalibSum>221.5) nGenSpecLeft  = (Int_t)(fZEMsp->Eval(zemenergy));
437   if(ZN1CalibSum+ZP1CalibSum>220.)  impPar    =  fZEMb->Eval(zemenergy);
438   */
439   //
440   //
441   // *** RECONSTRUCTION FROM REAL DATA
442   //
443   Float_t corrADCZEMHG = corrADCZEM1HG + corrADCZEM2HG;
444   //
445   if(corrADCZEMHG > supValueZEM){
446     nGenSpecNLeft  = (Int_t) (fZNCen->Eval(calibSumZN1HG));
447     nGenSpecPLeft  = (Int_t) (fZPCen->Eval(calibSumZP1HG));
448     nGenSpecLeft   = (Int_t) (fZDCCen->Eval(calibSumZN1HG+calibSumZP1HG));
449     nGenSpecNRight = (Int_t) (fZNCen->Eval(calibSumZN2HG));
450     nGenSpecPRight = (Int_t) (fZNCen->Eval(calibSumZP2HG));
451     nGenSpecRight  = (Int_t) (fZNCen->Eval(calibSumZN2HG+calibSumZP2HG));
452     impPar  = fbCen->Eval(calibSumZN1HG+calibSumZP1HG);
453   }
454   else if(corrADCZEMHG < infValueZEM){
455     nGenSpecNLeft = (Int_t) (fZNPer->Eval(calibSumZN1HG)); 
456     nGenSpecPLeft = (Int_t) (fZPPer->Eval(calibSumZP1HG));
457     nGenSpecLeft  = (Int_t) (fZDCPer->Eval(calibSumZN1HG+calibSumZP1HG));
458     impPar   = fbPer->Eval(calibSumZN1HG+calibSumZP1HG);
459   }
460   else if(corrADCZEMHG >= infValueZEM && corrADCZEMHG <= supValueZEM){
461     nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
462     nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
463     nGenSpecLeft  = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
464     impPar   =  fZEMb->Eval(corrADCZEMHG);
465   }
466   // 
467   if(calibSumZN1HG/maxValEZN1>1.)  nGenSpecNLeft = (Int_t) (fZEMn->Eval(corrADCZEMHG));
468   if(calibSumZP1HG/maxValEZP1>1.)  nGenSpecPLeft = (Int_t) (fZEMp->Eval(corrADCZEMHG));
469   if((calibSumZN1HG+calibSumZP1HG/maxValEZDC1)>1.){
470      nGenSpecLeft = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
471      impPar = fZEMb->Eval(corrADCZEMHG);
472   }
473   if(calibSumZN2HG/maxValEZN2>1.)  nGenSpecNRight = (Int_t) (fZEMn->Eval(corrADCZEMHG));
474   if(calibSumZP2HG/maxValEZP2>1.)  nGenSpecPRight = (Int_t) (fZEMp->Eval(corrADCZEMHG));
475   if((calibSumZN2HG+calibSumZP2HG/maxValEZDC2)>1.) nGenSpecRight = (Int_t)(fZEMsp->Eval(corrADCZEMHG));
476   //
477   if(nGenSpecNLeft>125)    nGenSpecNLeft=125;
478   else if(nGenSpecNLeft<0) nGenSpecNLeft=0;
479   if(nGenSpecPLeft>82)     nGenSpecPLeft=82;
480   else if(nGenSpecPLeft<0) nGenSpecPLeft=0;
481   if(nGenSpecLeft>207)     nGenSpecLeft=207;
482   else if(nGenSpecLeft<0)  nGenSpecLeft=0;
483   
484   //  ---      Number of generated participants (from HIJING parameterization)
485   Int_t nPart, nPartTotLeft, nPartTotRight;
486   nPart = 207-nGenSpecNLeft-nGenSpecPLeft;
487   nPartTotLeft = 207-nGenSpecLeft;
488   nPartTotRight = 207-nGenSpecRight;
489   if(nPart<0) nPart=0;
490   if(nPartTotLeft<0) nPartTotLeft=0;
491   if(nPartTotRight<0) nPartTotRight=0;
492   //
493   // *** DEBUG ***
494   printf("\n\t AliZDCReconstructor -> calibSumZN1HG %1.0f, calibSumZP1HG %1.0f,"
495       "  calibSumZN2HG %1.0f, calibSumZP2HG %1.0f, corrADCZEMHG %1.0f\n", 
496       calibSumZN1HG,calibSumZP1HG,calibSumZN2HG,calibSumZP2HG,corrADCZEMHG);
497   printf("\t AliZDCReconstructor -> nGenSpecNLeft %d, nGenSpecPLeft %d, nGenSpecLeft %d\n"
498       "\t\t nGenSpecNRight %d, nGenSpecPRight %d, nGenSpecRight %d\n", 
499       nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, 
500       nGenSpecNRight, nGenSpecPRight, nGenSpecRight);
501   printf("\t AliZDCReconstructor ->  NpartL %d,  NpartR %d,  b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
502
503   // create the output tree
504   AliZDCReco reco(calibSumZN1HG, calibSumZP1HG, calibSumZN2HG, calibSumZP2HG, 
505                   calibTowZN1LG, calibTowZN2LG, calibTowZP1LG, calibTowZP2LG, 
506                   corrADCZEM1HG, corrADCZEM2HG,
507                   nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight, 
508                   nGenSpecNLeft, nGenSpecPLeft, nGenSpecLeft, nGenSpecNRight, 
509                   nGenSpecPRight, nGenSpecRight,
510                   nPartTotLeft, nPartTotRight, impPar);
511                   
512   AliZDCReco* preco = &reco;
513   const Int_t kBufferSize = 4000;
514   clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
515
516   // write the output tree
517   clustersTree->Fill();
518 }
519
520 //_____________________________________________________________________________
521 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
522 {
523   // fill energies and number of participants to the ESD
524
525   AliZDCReco reco;
526   AliZDCReco* preco = &reco;
527   clustersTree->SetBranchAddress("ZDC", &preco);
528
529   clustersTree->GetEntry(0);
530   /*Double_t tZN1Ene[4], tZN2Ene[4];
531   for(Int_t i=0; i<4; i++){
532      tZN1Ene[i] = reco.GetZN1EnTow(i);
533      tZN2Ene[i] = reco.GetZN2EnTow(i);
534   }
535   esd->SetZDC(tZN1Ene, tZN2Ene, reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEMsignal(),
536               reco.GetZN2Energy(), reco.GetZP2Energy(), 
537               reco.GetNPartLeft());
538   */
539   esd->SetZDC(reco.GetZN1Energy(), reco.GetZP1Energy(), reco.GetZEMsignal(), 
540               reco.GetZN2Energy(), reco.GetZP2Energy(), 
541               reco.GetNPartLeft());
542   
543   /*Double_t tZN1Ene[4], tZN2Ene[4];
544   for(Int_t i=0; i<4; i++){
545      tZN1Ene[i] = reco.GetZN1EnTow(i);
546      tZN2Ene[i] = reco.GetZN2EnTow(i);
547   }*/
548   //esd->SetZN1TowerEnergy(tZN1Ene);
549   //esd->SetZN2TowerEnergy(tZN2Ene);
550   
551 }
552
553 //_____________________________________________________________________________
554 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri) 
555 {
556   // Setting the storage
557
558   Bool_t deleteManager = kFALSE;
559   
560   AliCDBManager *manager = AliCDBManager::Instance();
561   AliCDBStorage *defstorage = manager->GetDefaultStorage();
562   
563   if(!defstorage || !(defstorage->Contains("ZDC"))){ 
564      AliWarning("No default storage set or default storage doesn't contain ZDC!");
565      manager->SetDefaultStorage(uri);
566      deleteManager = kTRUE;
567   }
568  
569   AliCDBStorage *storage = manager->GetDefaultStorage();
570
571   if(deleteManager){
572     AliCDBManager::Instance()->UnsetDefaultStorage();
573     defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
574   }
575
576   return storage; 
577 }
578
579 //_____________________________________________________________________________
580 AliZDCCalibData* AliZDCReconstructor::GetCalibData() const
581 {
582
583   // Getting calibration object for ZDC set
584
585   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Data");
586   if(!entry) AliFatal("No calibration data loaded!");  
587
588   AliZDCCalibData *calibdata = dynamic_cast<AliZDCCalibData*>  (entry->GetObject());
589   if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
590
591   return calibdata;
592 }