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