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