1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // ************** Class for ZDC reconstruction ************** //
21 // Author: Chiara.Oppedisano@to.infn.it //
23 // NOTATIONS ADOPTED TO IDENTIFY DETECTORS (used in different ages!): //
24 // (ZN1,ZP1) or (ZNC, ZPC) or RIGHT refers to side C (RB26) //
25 // (ZN2,ZP2) or (ZNA, ZPA) or LEFT refers to side A (RB24) //
27 ///////////////////////////////////////////////////////////////////////////////
35 #include "AliRawReader.h"
36 #include "AliGRPObject.h"
37 #include "AliESDEvent.h"
38 #include "AliESDZDC.h"
39 #include "AliZDCDigit.h"
40 #include "AliZDCRawStream.h"
41 #include "AliZDCReco.h"
42 #include "AliZDCReconstructor.h"
43 #include "AliZDCPedestals.h"
44 #include "AliZDCEnCalib.h"
45 #include "AliZDCTowerCalib.h"
46 #include "AliZDCRecoParam.h"
47 #include "AliZDCRecoParampp.h"
48 #include "AliZDCRecoParamPbPb.h"
49 #include "AliRunInfo.h"
52 ClassImp(AliZDCReconstructor)
53 AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0; //reconstruction parameters
55 //_____________________________________________________________________________
56 AliZDCReconstructor:: AliZDCReconstructor() :
57 fPedData(GetPedestalData()),
58 fEnCalibData(GetEnergyCalibData()),
59 fTowCalibData(GetTowerCalibData()),
63 fIsCalibrationMB(kFALSE),
67 // **** Default constructor
72 //_____________________________________________________________________________
73 AliZDCReconstructor::~AliZDCReconstructor()
76 // if(fRecoParam) delete fRecoParam;
77 if(fPedData) delete fPedData;
78 if(fEnCalibData) delete fEnCalibData;
79 if(fTowCalibData) delete fTowCalibData;
82 //____________________________________________________________________________
83 void AliZDCReconstructor::Init()
85 // Setting reconstruction mode
86 // Getting beam type and beam energy from GRP calibration object
88 TString runType = GetRunInfo()->GetRunType();
90 if((runType.CompareTo("CALIBRATION_MB")) == 0){
91 fIsCalibrationMB = kTRUE;
94 TString beamType = GetRunInfo()->GetBeamType();
95 // This is a temporary solution to allow reconstruction in tests without beam
96 if(((beamType.CompareTo("UNKNOWN"))==0) && ((runType.CompareTo("PHYSICS")) == 0)){
99 else if(beamType==AliGRPObject::GetInvalidString()){
100 AliWarning("GRP/GRP/Data entry: missing value for the beam type !");
101 AliError("\t ZDC does not reconstruct event 4 UNKNOWN beam type\n");
105 if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
106 ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)) fRecoMode=1;
107 else if((beamType.CompareTo("A-A")) == 0){
109 if(fIsCalibrationMB == kTRUE){
110 fRecoParam = new AliZDCRecoParamPbPb();
112 TH2F* hZDCvsZEM = new TH2F("hZDCvsZEM","hZDCvsZEM",100,0.,10.,100,0.,1000.);
113 hZDCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCvsZEM->SetYTitle("E_{ZDC} (TeV)");
114 fRecoParam->SetZDCvsZEM(hZDCvsZEM);
116 TH2F* hZDCCvsZEM = new TH2F("hZDCCvsZEM","hZDCCvsZEM",100,0.,10.,100,0.,500.);
117 hZDCCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCCvsZEM->SetYTitle("E_{ZDCC} (TeV)");
118 fRecoParam->SetZDCCvsZEM(hZDCCvsZEM);
120 TH2F* hZDCAvsZEM = new TH2F("hZDCAvsZEM","hZDCAvsZEM",100,0.,10.,100,0.,500.);
121 hZDCAvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCAvsZEM->SetYTitle("E_{ZDCA} (TeV)");
122 fRecoParam->SetZDCAvsZEM(hZDCAvsZEM);
124 AliInfo("\n ***** CALIBRATION_MB data -> building AliZDCRecoParamPbPb object *****");
128 fBeamEnergy = GetRunInfo()->GetBeamEnergy();
129 if(fBeamEnergy==AliGRPObject::GetInvalidFloat()){
130 AliWarning("GRP/GRP/Data entry: missing value for the beam energy ! Using 0.");
134 if(fIsCalibrationMB==kFALSE)
135 printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f GeV *****\n\n",beamType.Data(), fBeamEnergy);
139 //_____________________________________________________________________________
140 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
142 // *** Local ZDC reconstruction for digits
143 // Works on the current event
145 // Retrieving calibration data
146 // Parameters for mean value pedestal subtraction
148 Float_t meanPed[2*kNch];
149 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
150 // Parameters pedestal subtraction through correlation with out-of-time signals
151 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
152 for(Int_t jj=0; jj<2*kNch; jj++){
153 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
154 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
159 AliZDCDigit* pdigit = &digit;
160 digitsTree->SetBranchAddress("ZDC", &pdigit);
161 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
164 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
165 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
166 for(Int_t i=0; i<10; i++){
167 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
168 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
171 Int_t digNentries = digitsTree->GetEntries();
172 Float_t ootDigi[kNch];
173 // -- Reading out-of-time signals (last kNch entries) for current event
175 for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
176 ootDigi[iDigit] = digitsTree->GetEntry(iDigit);
180 for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
181 digitsTree->GetEntry(iDigit);
182 if (!pdigit) continue;
184 Int_t det = digit.GetSector(0);
185 Int_t quad = digit.GetSector(1);
187 Float_t ped2SubHg=0., ped2SubLg=0.;
189 if(det==1) pedindex = quad;
190 else if(det==2) pedindex = quad+5;
191 else if(det==3) pedindex = quad+9;
192 else if(det==4) pedindex = quad+12;
193 else if(det==5) pedindex = quad+17;
195 else pedindex = (det-1)/3+22;
198 ped2SubHg = meanPed[pedindex];
199 ped2SubLg = meanPed[pedindex+kNch];
201 else if(fPedSubMode==1){
202 ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
203 ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
206 if(quad != 5){ // ZDC (not reference PTMs!)
207 if(det == 1){ // *** ZNC
208 tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
209 tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
210 if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
211 if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
213 else if(det == 2){ // *** ZP1
214 tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
215 tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
216 if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
217 if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
220 if(quad == 1){ // *** ZEM1
221 dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
222 dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
223 if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
224 if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
226 else if(quad == 2){ // *** ZEM2
227 dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
228 dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
229 if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
230 if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
233 else if(det == 4){ // *** ZN2
234 tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
235 tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
236 if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
237 if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
239 else if(det == 5){ // *** ZP2
240 tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
241 tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
242 if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
243 if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
246 else{ // Reference PMs
248 sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
249 sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
251 if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
252 if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
255 sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
256 sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
258 if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
259 if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
264 /*printf("AliZDCReconstructor: digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
265 iDigit, det, quad, ped2SubHg, ped2SubLg);
266 printf(" -> pedindex %d\n", pedindex);
267 printf(" HGChain -> RawDig %d DigCorr %1.2f",
268 digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg);
269 printf(" LGChain -> RawDig %d DigCorr %1.2f\n",
270 digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);*/
275 // If CALIBRATION_MB run build the RecoParam object
276 if(fIsCalibrationMB){
277 Float_t ZDCC=0., ZDCA=0., ZEM=0;
278 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
279 for(Int_t jkl=0; jkl<5; jkl++){
280 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
281 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
283 // Using energies in TeV in fRecoParam object
284 BuildRecoParam(ZDCC/1000., ZDCA/1000., ZEM/1000.);
286 Bool_t recFlag1 = kFALSE, recFlag2 = kFALSE, recFlag3 = kFALSE;
287 // reconstruct the event
289 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
290 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, recFlag1, recFlag2, recFlag3);
291 else if(fRecoMode==2)
292 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
293 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, recFlag1, recFlag2, recFlag3);
296 //_____________________________________________________________________________
297 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
299 // *** ZDC raw data reconstruction
300 // Works on the current event
302 // Retrieving calibration data
303 // Parameters for pedestal subtraction
305 Float_t meanPed[2*kNch];
306 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
307 // Parameters pedestal subtraction through correlation with out-of-time signals
308 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
309 for(Int_t jj=0; jj<2*kNch; jj++){
310 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
311 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
314 Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
315 Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
316 Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
317 Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
318 Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
319 Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
320 for(Int_t ich=0; ich<5; ich++){
321 adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
322 adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
323 adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
324 adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
326 adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
327 pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
331 // loop over raw data
332 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
333 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
334 for(Int_t i=0; i<10; i++){
335 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
336 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
339 //fNRun = (Int_t) rawReader->GetRunNumber();
340 Bool_t chOff=kFALSE, isUndflw=kFALSE, isOvflw=kFALSE, isADCEvGood=kFALSE;
341 Int_t kFirstADCGeo=0, kLastADCGeo=3;
342 //Int_t kScalerGeo=8, kPUGeo=29, kTrigScales=30, kTrigHistory=31;
345 AliZDCRawStream rawData(rawReader);
346 while(rawData.Next()){
347 //if(rawData.IsCalibration() == kFALSE){ // Reading ADCs
348 if((rawData.GetADCModule()>=kFirstADCGeo) && (rawData.GetADCModule()<=kLastADCGeo)){
349 //printf(" **** Reading ADC raw data from module %d **** \n",rawData.GetADCModule());
351 if((rawData.IsADCDataWord()) && (rawData.GetNChannelsOn()<48)) chOff = kTRUE;
352 if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)) isUndflw = kTRUE;
353 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)) isOvflw = kTRUE;
354 if((rawData.IsADCDataWord()) && (rawData.IsADCEventGood() == kTRUE)) isADCEvGood = kTRUE;
356 if((rawData.IsADCDataWord()) && (isUndflw==kFALSE)
357 && (isOvflw==kFALSE) && (isADCEvGood==kTRUE)){
359 Int_t adcMod = rawData.GetADCModule();
360 Int_t det = rawData.GetSector(0);
361 Int_t quad = rawData.GetSector(1);
362 Int_t gain = rawData.GetADCGain();
365 // Mean pedestal value subtraction -------------------------------------------------------
366 if(fPedSubMode == 0){
367 // Not interested in o.o.t. signals (ADC modules 2, 3)
368 if(adcMod == 2 || adcMod == 3) continue;
370 if(quad != 5){ // ZDCs (not reference PTMs)
373 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
374 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
378 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
379 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
384 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
385 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
388 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
389 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
394 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
395 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
399 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
400 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
403 else{ // reference PM
404 pedindex = (det-1)/3 + 22;
406 if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
407 else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
410 if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
411 else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
415 /*printf(" AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f\n",
416 det,quad,gain, pedindex, meanPed[pedindex]);
417 printf(" RawADC %d ADCCorr %1.0f\n",
418 rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);*/
420 }// mean pedestal subtraction
421 // Pedestal subtraction from correlation ------------------------------------------------
422 else if(fPedSubMode == 1){
424 if(adcMod==0 || adcMod==1){
425 if(quad != 5){ // signals from ZDCs
427 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
428 else adcZN1lg[quad] = rawData.GetADCValue();
431 if(gain==0) adcZP1[quad] = rawData.GetADCValue();
432 else adcZP1lg[quad] = rawData.GetADCValue();
435 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
436 else adcZEMlg[quad-1] = rawData.GetADCValue();
439 if(gain==0) adcZN2[quad] = rawData.GetADCValue();
440 else adcZN2lg[quad] = rawData.GetADCValue();
443 if(gain==0) adcZP2[quad] = rawData.GetADCValue();
444 else adcZP2lg[quad] = rawData.GetADCValue();
447 else{ // signals from reference PM
448 if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
449 else pmReflg[quad-1] = rawData.GetADCValue();
452 // Out-of-time pedestals
453 else if(adcMod==2 || adcMod==3){
454 if(quad != 5){ // signals from ZDCs
456 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
457 else adcZN1ootlg[quad] = rawData.GetADCValue();
460 if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
461 else adcZP1ootlg[quad] = rawData.GetADCValue();
464 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
465 else adcZEMootlg[quad-1] = rawData.GetADCValue();
468 if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
469 else adcZN2ootlg[quad] = rawData.GetADCValue();
472 if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
473 else adcZP2ootlg[quad] = rawData.GetADCValue();
476 else{ // signals from reference PM
477 if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
478 else pmRefootlg[quad-1] = rawData.GetADCValue();
481 } // pedestal subtraction from correlation
483 /*printf("\t AliZDCReconstructor: det %d quad %d res %d -> Ped[%d] = %1.0f\n",
484 det,quad,gain, pedindex, meanPed[pedindex]);*/
486 }// Not raw data from calibration run!
492 for(Int_t t=0; t<5; t++){
493 tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
494 tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
496 tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
497 tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
499 tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
500 tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
502 tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
503 tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
504 // 0---------0 Ch. debug 0---------0
505 /* printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
506 printf("\tCorrCoeff0\tCorrCoeff1\n");
507 printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
508 printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
509 printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
510 printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
511 printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
512 printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
513 printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
514 printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
516 printf("ZN1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
517 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
518 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
519 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
521 printf("ZP1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
522 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
523 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
524 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
526 printf("ZN2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
527 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
528 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
529 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
531 printf("ZP2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
532 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
533 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
534 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
537 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
538 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
539 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
540 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
542 sPMRef1[0] = pmRef[0] - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
543 sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
544 sPMRef2[0] = pmRef[0] - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
545 sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
548 // If CALIBRATION_MB run build the RecoParam object
549 if(fIsCalibrationMB){
550 Float_t ZDCC=0., ZDCA=0., ZEM=0;
551 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
552 for(Int_t jkl=0; jkl<5; jkl++){
553 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
554 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
556 BuildRecoParam(ZDCC/100., ZDCA/100., ZEM/100.);
558 // reconstruct the event
560 if(fRecoMode==1) // p-p data
561 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
562 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, chOff, isUndflw, isOvflw);
563 else if(fRecoMode==2) // Pb-Pb data
564 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
565 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, chOff, isUndflw, isOvflw);
569 //_____________________________________________________________________________
570 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* corrADCZN1,
571 Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
572 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
573 Bool_t channelsOff, Bool_t chUnderflow, Bool_t chOverflow) const
575 // ****************** Reconstruct one event ******************
578 /*printf("\n*************************************************\n");
579 printf(" ReconstructEventpp -> values after pedestal subtraction:\n");
580 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
581 corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
582 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
583 corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
584 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
585 corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
586 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
587 corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
588 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
589 printf("*************************************************\n");*/
591 // ---- Setting reco flags
594 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
596 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
597 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
598 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
599 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
600 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
601 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
603 if(channelsOff==kTRUE) rFlags[8] = 0x1;
604 if(chUnderflow == kTRUE) rFlags[9] = 0x1;
605 if(chOverflow==kTRUE) rFlags[10] = 0x1;
606 recoFlag = rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
607 rFlags[5] << 5 | rFlags[4] << 4 | rFlags[3] << 3 |
608 rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
610 // ****** Retrieving calibration data
611 // --- Equalization coefficients ---------------------------------------------
612 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
613 for(Int_t ji=0; ji<5; ji++){
614 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
615 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
616 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
617 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
619 // --- Energy calibration factors ------------------------------------
621 // **** Energy calibration coefficient set to 1
622 // **** (no trivial way to calibrate in p-p runs)
623 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
625 // ****** Equalization of detector responses
626 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
627 for(Int_t gi=0; gi<10; gi++){
629 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
630 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
631 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
632 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
635 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
636 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
637 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
638 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
642 /*printf("\n ------------- EQUALIZATION -------------\n");
643 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
644 equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
645 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
646 equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
647 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
648 equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
649 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
650 equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
651 printf(" ----------------------------------------\n");*/
653 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
654 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
655 for(Int_t gi=0; gi<5; gi++){
656 calibSumZN1[0] += equalTowZN1[gi];
657 calibSumZP1[0] += equalTowZP1[gi];
658 calibSumZN2[0] += equalTowZN2[gi];
659 calibSumZP2[0] += equalTowZP2[gi];
661 calibSumZN1[1] += equalTowZN1[gi+5];
662 calibSumZP1[1] += equalTowZP1[gi+5];
663 calibSumZN2[1] += equalTowZN2[gi+5];
664 calibSumZP2[1] += equalTowZP2[gi+5];
667 calibSumZN1[0] = calibSumZN1[0]*calibEne[0];
668 calibSumZP1[0] = calibSumZP1[0]*calibEne[1];
669 calibSumZN2[0] = calibSumZN2[0]*calibEne[2];
670 calibSumZP2[0] = calibSumZP2[0]*calibEne[3];
672 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
673 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
674 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
675 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
677 // ****** Energy calibration of detector responses
678 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
679 for(Int_t gi=0; gi<5; gi++){
681 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
682 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
683 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
684 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
686 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
687 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
688 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
689 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
692 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
693 calibZEM1[0] = corrADCZEM1[0]*calibEne[5];
694 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
695 calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
696 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
697 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
699 /*printf("\n ------------- CALIBRATION -------------\n");
700 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
701 calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
702 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
703 calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
704 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
705 calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
706 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
707 calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
708 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
709 printf(" ----------------------------------------\n");*/
711 // ****** No. of spectator and participants nucleons
712 // Variables calculated to comply with ESD structure
713 // *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
714 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
715 Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
716 Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
717 Double_t impPar=0., impPar1=0., impPar2=0.;
719 // create the output tree
720 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
721 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
722 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
723 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
724 nGenSpec, nGenSpecLeft, nGenSpecRight,
725 nPart, nPartTotLeft, nPartTotRight,
726 impPar, impPar1, impPar2,
729 const Int_t kBufferSize = 4000;
730 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
731 // write the output tree
732 clustersTree->Fill();
735 //_____________________________________________________________________________
736 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
737 Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
738 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
739 Bool_t channelsOff, Bool_t chUnderflow, Bool_t chOverflow) const
741 // ****************** Reconstruct one event ******************
743 // ---- Setting reco flags
746 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
748 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
749 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
750 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
751 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
752 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
753 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
755 if(channelsOff==kTRUE) rFlags[8] = 0x1;
756 if(chUnderflow == kTRUE) rFlags[9] = 0x1;
757 if(chOverflow==kTRUE) rFlags[10] = 0x1;
758 recoFlag = rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
759 rFlags[5] << 5 | rFlags[4] << 4 | rFlags[3] << 3 |
760 rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
762 // ****** Retrieving calibration data
763 // --- Equalization coefficients ---------------------------------------------
764 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
765 for(Int_t ji=0; ji<5; ji++){
766 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
767 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
768 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
769 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
771 // --- Energy calibration factors ------------------------------------
772 Float_t valFromOCDB[6], calibEne[6];
773 for(Int_t ij=0; ij<6; ij++){
774 valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
776 if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
777 else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
779 else calibEne[ij] = valFromOCDB[ij];
782 // ****** Equalization of detector responses
783 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
784 for(Int_t gi=0; gi<10; gi++){
786 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
787 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
788 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
789 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
792 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
793 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
794 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
795 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
799 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
800 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
801 for(Int_t gi=0; gi<5; gi++){
802 calibSumZN1[0] += equalTowZN1[gi];
803 calibSumZP1[0] += equalTowZP1[gi];
804 calibSumZN2[0] += equalTowZN2[gi];
805 calibSumZP2[0] += equalTowZP2[gi];
807 calibSumZN1[1] += equalTowZN1[gi+5];
808 calibSumZP1[1] += equalTowZP1[gi+5];
809 calibSumZN2[1] += equalTowZN2[gi+5];
810 calibSumZP2[1] += equalTowZP2[gi+5];
813 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]*8.;
814 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]*8.;
815 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]*8.;
816 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]*8.;
818 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
819 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
820 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
821 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
823 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
824 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]*8.;
825 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
826 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]*8.;
827 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
828 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
830 // ****** Energy calibration of detector responses
831 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
832 for(Int_t gi=0; gi<5; gi++){
834 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]*8.;
835 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]*8.;
836 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]*8.;
837 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]*8.;
839 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
840 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
841 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
842 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
845 // ****** Number of detected spectator nucleons
846 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
848 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
849 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
850 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
851 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
853 else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
854 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
855 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
856 nDetSpecNRight, nDetSpecPRight);*/
858 Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
859 Int_t nPart=0, nPartA=0, nPartC=0;
860 Double_t b=0., bA=0., bC=0.;
862 if(fIsCalibrationMB == kFALSE){
863 // ****** Reconstruction parameters ------------------
865 //fRecoParam->Print("");
868 if (!fRecoParam) fRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
870 TH2F *hZDCvsZEM = fRecoParam->GethZDCvsZEM();
871 TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
872 TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
873 TH1D *hNpartDist = fRecoParam->GethNpartDist();
874 TH1D *hbDist = fRecoParam->GethbDist();
875 Float_t ClkCenter = fRecoParam->GetClkCenter();
877 Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
878 Double_t origin = xHighEdge*ClkCenter;
880 //printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
882 // ====> Summed ZDC info (sideA+side C)
883 TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
884 Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
885 Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
886 line->SetParameter(0, y/(x-origin));
887 line->SetParameter(1, -origin*y/(x-origin));
889 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
890 //printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
892 Double_t countPerc=0;
893 Double_t xBinCenter=0, yBinCenter=0;
894 for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
895 for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
896 xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
897 yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
899 if(line->GetParameter(0)>0){
900 if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
901 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
903 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
904 xBinCenter, yBinCenter, countPerc);*/
908 if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
909 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
911 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
912 xBinCenter, yBinCenter, countPerc);*/
918 Double_t xSecPerc = 0.;
919 if(hZDCvsZEM->GetEntries()!=0){
920 xSecPerc = countPerc/hZDCvsZEM->GetEntries();
923 AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
926 //printf(" xSecPerc %1.4f \n", xSecPerc);
929 TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
930 Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
931 lineC->SetParameter(0, yC/(x-origin));
932 lineC->SetParameter(1, -origin*yC/(x-origin));
934 //printf(" ***************** Side C \n");
935 //printf(" E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
937 Double_t countPercC=0;
938 Double_t xBinCenterC=0, yBinCenterC=0;
939 for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
940 for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
941 xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
942 yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
943 if(lineC->GetParameter(0)>0){
944 if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
945 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
949 if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
950 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
956 Double_t xSecPercC = 0.;
957 if(hZDCCvsZEM->GetEntries()!=0){
958 xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
961 AliWarning(" Histogram hZDCCvsZEM from OCDB has no entries!!!");
964 //printf(" xSecPercC %1.4f \n", xSecPercC);
967 TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
968 Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
969 lineA->SetParameter(0, yA/(x-origin));
970 lineA->SetParameter(1, -origin*yA/(x-origin));
973 //printf(" ***************** Side A \n");
974 //printf(" E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
976 Double_t countPercA=0;
977 Double_t xBinCenterA=0, yBinCenterA=0;
978 for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
979 for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
980 xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
981 yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
982 if(lineA->GetParameter(0)>0){
983 if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
984 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
988 if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
989 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
995 Double_t xSecPercA = 0.;
996 if(hZDCAvsZEM->GetEntries()!=0){
997 xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
1000 AliWarning(" Histogram hZDCAvsZEM from OCDB has no entries!!!");
1003 //printf(" xSecPercA %1.4f \n", xSecPercA);
1005 // ****** Number of participants (from E_ZDC vs. E_ZEM correlation)
1006 Int_t nPart=0, nPartC=0, nPartA=0;
1007 Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1008 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1009 nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1010 if((1.-nPartFrac) < xSecPerc){
1011 nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1013 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
1014 //printf(" nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1018 if(nPart<0) nPart=0;
1020 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1021 nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1022 if((1.-nPartFracC) < xSecPercC){
1023 nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1025 //printf(" ***************** Side C \n");
1026 //printf(" nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1030 if(nPartC<0) nPartC=0;
1032 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1033 nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1034 if((1.-nPartFracA) < xSecPercA){
1035 nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1037 //printf(" ***************** Side A \n");
1038 //printf(" nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1042 if(nPartA<0) nPartA=0;
1044 // ****** Impact parameter (from E_ZDC vs. E_ZEM correlation)
1045 Float_t b=0, bC=0, bA=0;
1046 Double_t bFrac=0., bFracC=0., bFracA=0.;
1047 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1048 bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1049 if((1.-bFrac) < xSecPerc){
1050 b = hbDist->GetBinLowEdge(ibbin);
1055 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1056 bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1057 if((1.-bFracC) < xSecPercC){
1058 bC = hbDist->GetBinLowEdge(ibbin);
1063 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1064 bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1065 if((1.-bFracA) < xSecPercA){
1066 bA = hbDist->GetBinLowEdge(ibbin);
1071 // ****** Number of spectator nucleons
1072 Int_t nGenSpec=0, nGenSpecC=0, nGenSpecA=0;
1074 nGenSpec = 416 - nPart;
1075 nGenSpecC = 416 - nPartC;
1076 nGenSpecA = 416 - nPartA;
1077 if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1078 if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1079 if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1082 /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1083 " calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
1084 calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1085 printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n",
1086 nGenSpecLeft, nGenSpecRight);
1087 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1090 delete lineC; delete lineA;
1092 } // ONLY IF fIsCalibrationMB==kFALSE
1094 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1095 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1096 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1097 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1098 nGenSpec, nGenSpecA, nGenSpecC,
1099 nPart, nPartA, nPartC, b, bA, bC,
1102 const Int_t kBufferSize = 4000;
1103 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1105 // write the output tree
1106 clustersTree->Fill();
1109 //_____________________________________________________________________________
1110 void AliZDCReconstructor::BuildRecoParam(Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1112 // Calculate RecoParam object for Pb-Pb data
1113 (fRecoParam->GethZDCvsZEM())->Fill(ZDCC+ZDCA, ZEM);
1114 (fRecoParam->GethZDCCvsZEM())->Fill(ZDCC, ZEM);
1115 (fRecoParam->GethZDCAvsZEM())->Fill(ZDCA, ZEM);
1119 //_____________________________________________________________________________
1120 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1122 // fill energies and number of participants to the ESD
1124 if(fIsCalibrationMB==kTRUE) WritePbPbRecoParamInOCDB();
1127 AliZDCReco* preco = &reco;
1128 clustersTree->SetBranchAddress("ZDC", &preco);
1130 clustersTree->GetEntry(0);
1132 AliESDZDC * esdzdc = esd->GetESDZDC();
1133 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1134 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1135 for(Int_t i=0; i<5; i++){
1136 tZN1Ene[i] = reco.GetZN1HREnTow(i);
1137 tZN2Ene[i] = reco.GetZN2HREnTow(i);
1138 tZP1Ene[i] = reco.GetZP1HREnTow(i);
1139 tZP2Ene[i] = reco.GetZP2HREnTow(i);
1141 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1142 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1143 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1144 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1147 esdzdc->SetZN1TowerEnergy(tZN1Ene);
1148 esdzdc->SetZN2TowerEnergy(tZN2Ene);
1149 esdzdc->SetZP1TowerEnergy(tZP1Ene);
1150 esdzdc->SetZP2TowerEnergy(tZP2Ene);
1152 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1153 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1154 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1155 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1157 Int_t nPart = reco.GetNParticipants();
1158 Int_t nPartA = reco.GetNPartSideA();
1159 Int_t nPartC = reco.GetNPartSideC();
1160 Double_t b = reco.GetImpParameter();
1161 Double_t bA = reco.GetImpParSideA();
1162 Double_t bC = reco.GetImpParSideC();
1163 UInt_t recoFlag = reco.GetRecoFlag();
1165 esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
1166 reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1167 nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1171 //_____________________________________________________________________________
1172 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1174 // Setting the storage
1176 Bool_t deleteManager = kFALSE;
1178 AliCDBManager *manager = AliCDBManager::Instance();
1179 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1181 if(!defstorage || !(defstorage->Contains("ZDC"))){
1182 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1183 manager->SetDefaultStorage(uri);
1184 deleteManager = kTRUE;
1187 AliCDBStorage *storage = manager->GetDefaultStorage();
1190 AliCDBManager::Instance()->UnsetDefaultStorage();
1191 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1197 //_____________________________________________________________________________
1198 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1201 // Getting pedestal calibration object for ZDC set
1203 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1204 if(!entry) AliFatal("No calibration data loaded!");
1206 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
1207 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1212 //_____________________________________________________________________________
1213 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1216 // Getting energy and equalization calibration object for ZDC set
1218 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1219 if(!entry) AliFatal("No calibration data loaded!");
1221 AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1222 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1227 //_____________________________________________________________________________
1228 AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1231 // Getting energy and equalization calibration object for ZDC set
1233 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1234 if(!entry) AliFatal("No calibration data loaded!");
1236 AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1237 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1242 //_____________________________________________________________________________
1243 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1246 // Getting reconstruction parameters from OCDB
1248 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1249 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1251 AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1252 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1258 //_____________________________________________________________________________
1259 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1262 // Getting reconstruction parameters from OCDB
1264 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1265 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1267 AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1268 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1274 //_____________________________________________________________________________
1275 void AliZDCReconstructor::WritePbPbRecoParamInOCDB() const
1278 // Writing Pb-Pb reconstruction parameters from OCDB
1280 AliCDBManager *man = AliCDBManager::Instance();
1281 AliCDBMetaData *md= new AliCDBMetaData();
1282 md->SetResponsible("Chiara Oppedisano");
1283 md->SetComment("ZDC Pb-Pb reconstruction parameters");
1284 md->SetObjectClassName("AliZDCRecoParamPbPb");
1285 AliCDBId id("ZDC/Calib/RecoParamPbPb",fNRun,AliCDBRunRange::Infinity());
1286 man->Put(fRecoParam, id, md);