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