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 "AliESDEvent.h"
37 #include "AliESDZDC.h"
38 #include "AliZDCDigit.h"
39 #include "AliZDCRawStream.h"
40 #include "AliZDCReco.h"
41 #include "AliZDCReconstructor.h"
42 #include "AliZDCPedestals.h"
43 #include "AliZDCEnCalib.h"
44 #include "AliZDCTowerCalib.h"
45 #include "AliZDCMBCalib.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
54 AliZDCMBCalib *AliZDCReconstructor::fMBCalibData=0; //calibration parameters for A-A reconstruction
56 //_____________________________________________________________________________
57 AliZDCReconstructor:: AliZDCReconstructor() :
58 fPedData(GetPedestalData()),
59 fEnCalibData(GetEnergyCalibData()),
60 fTowCalibData(GetTowerCalibData()),
64 fIsCalibrationMB(kFALSE),
68 // **** 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;
80 if(fMBCalibData) delete fMBCalibData;
83 //____________________________________________________________________________
84 void AliZDCReconstructor::Init()
86 // Setting reconstruction mode
87 // Getting beam type and beam energy from GRP calibration object
89 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) &&
97 ((runType.CompareTo("PHYSICS"))==0 || (runType.CompareTo("CALIBRATION_BC"))==0)){
100 /*else if((beamType.CompareTo("UNKNOWN"))==0){
101 AliError("\t 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)){
109 else if((beamType.CompareTo("A-A")) == 0){
113 fBeamEnergy = GetRunInfo()->GetBeamEnergy();
114 if(fBeamEnergy==0.) AliWarning(" Beam energy value missing -> E_beam = 0");
116 if(fIsCalibrationMB==kFALSE)
117 printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f GeV *****\n\n",beamType.Data(), fBeamEnergy);
121 //_____________________________________________________________________________
122 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
124 // *** Local ZDC reconstruction for digits
125 // Works on the current event
127 // Retrieving calibration data
128 // Parameters for mean value pedestal subtraction
130 Float_t meanPed[2*kNch];
131 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
132 // Parameters pedestal subtraction through correlation with out-of-time signals
133 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
134 for(Int_t jj=0; jj<2*kNch; jj++){
135 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
136 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
141 AliZDCDigit* pdigit = &digit;
142 digitsTree->SetBranchAddress("ZDC", &pdigit);
143 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
146 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
147 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
148 for(Int_t i=0; i<10; i++){
149 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
150 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
153 Int_t digNentries = digitsTree->GetEntries();
154 Float_t ootDigi[kNch];
155 // -- Reading out-of-time signals (last kNch entries) for current event
157 for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
158 ootDigi[iDigit] = digitsTree->GetEntry(iDigit);
162 for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
163 digitsTree->GetEntry(iDigit);
164 if (!pdigit) continue;
166 Int_t det = digit.GetSector(0);
167 Int_t quad = digit.GetSector(1);
169 Float_t ped2SubHg=0., ped2SubLg=0.;
171 if(det==1) pedindex = quad;
172 else if(det==2) pedindex = quad+5;
173 else if(det==3) pedindex = quad+9;
174 else if(det==4) pedindex = quad+12;
175 else if(det==5) pedindex = quad+17;
177 else pedindex = (det-1)/3+22;
180 ped2SubHg = meanPed[pedindex];
181 ped2SubLg = meanPed[pedindex+kNch];
183 else if(fPedSubMode==1){
184 ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
185 ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
188 if(quad != 5){ // ZDC (not reference PTMs!)
189 if(det == 1){ // *** ZNC
190 tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
191 tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
192 if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
193 if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
195 else if(det == 2){ // *** ZP1
196 tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
197 tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
198 if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
199 if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
202 if(quad == 1){ // *** ZEM1
203 dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
204 dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
205 if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
206 if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
208 else if(quad == 2){ // *** ZEM2
209 dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
210 dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
211 if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
212 if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
215 else if(det == 4){ // *** ZN2
216 tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
217 tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
218 if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
219 if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
221 else if(det == 5){ // *** ZP2
222 tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
223 tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
224 if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
225 if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
228 else{ // Reference PMs
230 sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
231 sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
233 if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
234 if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
237 sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
238 sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
240 if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
241 if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
246 /*printf("AliZDCReconstructor: digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
247 iDigit, det, quad, ped2SubHg, ped2SubLg);
248 printf(" -> pedindex %d\n", pedindex);
249 printf(" HGChain -> RawDig %d DigCorr %1.2f",
250 digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg);
251 printf(" LGChain -> RawDig %d DigCorr %1.2f\n",
252 digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);*/
257 for(Int_t jj=0; jj<32; jj++) counts[jj]=0;
259 Int_t evQualityBlock[4] = {1,0,0,0};
260 Int_t triggerBlock[4] = {0,0,0,0};
261 Int_t chBlock[3] = {0,0,0};
264 // reconstruct the event
266 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
267 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
269 evQualityBlock, triggerBlock, chBlock, puBits);
270 else if(fRecoMode==2)
271 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
272 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
274 evQualityBlock, triggerBlock, chBlock, puBits);
277 //_____________________________________________________________________________
278 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
280 // *** ZDC raw data reconstruction
281 // Works on the current event
283 // Retrieving calibration data
284 // Parameters for pedestal subtraction
286 Float_t meanPed[2*kNch];
287 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
288 // Parameters pedestal subtraction through correlation with out-of-time signals
289 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
290 for(Int_t jj=0; jj<2*kNch; jj++){
291 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
292 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
295 Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
296 Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
297 Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
298 Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
299 Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
300 Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
301 for(Int_t ich=0; ich<5; ich++){
302 adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
303 adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
304 adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
305 adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
307 adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
308 pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
312 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
313 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
314 for(Int_t i=0; i<10; i++){
315 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
316 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
319 Bool_t isScalerOn=kFALSE;
321 UInt_t scalerData[32];
322 for(Int_t k=0; k<32; k++) scalerData[k]=0;
324 Int_t evQualityBlock[4] = {1,0,0,0};
325 Int_t triggerBlock[4] = {0,0,0,0};
326 Int_t chBlock[3] = {0,0,0};
329 //fNRun = (Int_t) rawReader->GetRunNumber();
330 Int_t kFirstADCGeo=0, kLastADCGeo=3, kScalerGeo=8, kPUGeo=29;
331 //Int_t kTrigScales=30, kTrigHistory=31;
333 // loop over raw data
335 AliZDCRawStream rawData(rawReader);
336 while(rawData.Next()){
338 // ***************************** Reading ADCs
339 if((rawData.GetADCModule()>=kFirstADCGeo) && (rawData.GetADCModule()<=kLastADCGeo)){
340 //printf(" **** Reading ADC raw data from module %d **** \n",rawData.GetADCModule());
342 if((rawData.IsADCDataWord()) && (rawData.GetNChannelsOn()<48)) chBlock[0] = kTRUE;
343 if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)) chBlock[1] = kTRUE;
344 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)) chBlock[2] = kTRUE;
345 if((rawData.IsADCDataWord()) && (rawData.IsADCEventGood() == kTRUE)) evQualityBlock[0] = kTRUE;
347 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow()==kFALSE)
348 && (rawData.IsOverflow()==kFALSE) && (rawData.IsADCEventGood()==kTRUE)){
350 Int_t adcMod = rawData.GetADCModule();
351 Int_t det = rawData.GetSector(0);
352 Int_t quad = rawData.GetSector(1);
353 Int_t gain = rawData.GetADCGain();
356 // Mean pedestal value subtraction -------------------------------------------------------
357 if(fPedSubMode == 0){
358 // Not interested in o.o.t. signals (ADC modules 2, 3)
359 if(adcMod == 2 || adcMod == 3) continue;
361 if(quad != 5){ // ZDCs (not reference PTMs)
364 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
365 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
369 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
370 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
375 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
376 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
379 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
380 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
385 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
386 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
390 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
391 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
394 else{ // reference PM
395 pedindex = (det-1)/3 + 22;
397 if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
398 else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
401 if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
402 else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
407 printf(" AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f",
408 det,quad,gain, pedindex, meanPed[pedindex]);
409 printf(" RawADC %d ADCCorr %1.0f\n",
410 rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);
412 }// mean pedestal subtraction
413 // Pedestal subtraction from correlation ------------------------------------------------
414 else if(fPedSubMode == 1){
416 if(adcMod==0 || adcMod==1){
417 if(quad != 5){ // signals from ZDCs
419 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
420 else adcZN1lg[quad] = rawData.GetADCValue();
423 if(gain==0) adcZP1[quad] = rawData.GetADCValue();
424 else adcZP1lg[quad] = rawData.GetADCValue();
427 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
428 else adcZEMlg[quad-1] = rawData.GetADCValue();
431 if(gain==0) adcZN2[quad] = rawData.GetADCValue();
432 else adcZN2lg[quad] = rawData.GetADCValue();
435 if(gain==0) adcZP2[quad] = rawData.GetADCValue();
436 else adcZP2lg[quad] = rawData.GetADCValue();
439 else{ // signals from reference PM
440 if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
441 else pmReflg[quad-1] = rawData.GetADCValue();
444 // Out-of-time pedestals
445 else if(adcMod==2 || adcMod==3){
446 if(quad != 5){ // signals from ZDCs
448 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
449 else adcZN1ootlg[quad] = rawData.GetADCValue();
452 if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
453 else adcZP1ootlg[quad] = rawData.GetADCValue();
456 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
457 else adcZEMootlg[quad-1] = rawData.GetADCValue();
460 if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
461 else adcZN2ootlg[quad] = rawData.GetADCValue();
464 if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
465 else adcZP2ootlg[quad] = rawData.GetADCValue();
468 else{ // signals from reference PM
469 if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
470 else pmRefootlg[quad-1] = rawData.GetADCValue();
473 } // pedestal subtraction from correlation
475 /*printf("\t AliZDCReconstructor: det %d quad %d res %d -> Ped[%d] = %1.0f\n",
476 det,quad,gain, pedindex, meanPed[pedindex]);*/
479 // ***************************** Reading Scaler
480 else if(rawData.GetADCModule()==kScalerGeo){
481 if(rawData.IsScalerWord()==kTRUE && rawData.IsScEventGood()==kTRUE){
483 scalerData[jsc] = rawData.GetTriggerCount();
485 //printf(" Reconstructed VME Scaler: %d %d ",jsc,scalerData[jsc]);
490 // ***************************** Reading PU
491 else if(rawData.GetADCModule()==kPUGeo){
492 puBits = rawData.GetDetectorPattern();
494 // ***************************** Reading trigger history
495 else if(rawData.IstriggerHistoryWord()==kTRUE){
496 triggerBlock[0] = rawData.IsCPTInputEMDTrigger();
497 triggerBlock[1] = rawData.IsCPTInputSemiCentralTrigger();
498 triggerBlock[2] = rawData.IsCPTInputCentralTrigger();
499 triggerBlock[3] = rawData.IsCPTInputMBTrigger();
505 for(Int_t t=0; t<5; t++){
506 tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
507 tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
509 tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
510 tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
512 tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
513 tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
515 tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
516 tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
517 // 0---------0 Ch. debug 0---------0
518 /* printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
519 printf("\tCorrCoeff0\tCorrCoeff1\n");
520 printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
521 printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
522 printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
523 printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
524 printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
525 printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
526 printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
527 printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
529 printf("ZN1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
530 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
531 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
532 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
534 printf("ZP1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
535 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
536 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
537 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
539 printf("ZN2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
540 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
541 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
542 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
544 printf("ZP2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
545 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
546 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
547 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
550 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
551 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
552 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
553 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
555 sPMRef1[0] = pmRef[0] - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
556 sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
557 sPMRef2[0] = pmRef[0] - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
558 sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
561 if(fRecoMode==1) // p-p data
562 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
563 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
564 isScalerOn, scalerData,
565 evQualityBlock, triggerBlock, chBlock, puBits);
566 else if(fRecoMode==2) // Pb-Pb data
567 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
568 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
569 isScalerOn, scalerData,
570 evQualityBlock, triggerBlock, chBlock, puBits);
573 //_____________________________________________________________________________
574 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* corrADCZN1,
575 Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
576 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
577 Bool_t isScalerOn, UInt_t* scaler,
578 Int_t* evQualityBlock, Int_t* triggerBlock, Int_t* chBlock, UInt_t puBits) const
580 // ****************** Reconstruct one event ******************
583 /*printf("\n*************************************************\n");
584 printf(" ReconstructEventpp -> values after pedestal subtraction:\n");
585 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
586 corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
587 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
588 corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
589 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
590 corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
591 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
592 corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
593 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
594 printf("*************************************************\n");*/
596 // ---------------------- Setting reco flags for ESD
598 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
600 if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
601 else rFlags[31] = 0x1;
603 if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
604 if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
605 if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
607 if(triggerBlock[0] == 1) rFlags[27] = 0x1;
608 if(triggerBlock[1] == 1) rFlags[26] = 0x1;
609 if(triggerBlock[2] == 1) rFlags[25] = 0x1;
610 if(triggerBlock[3] == 1) rFlags[24] = 0x1;
612 if(chBlock[0] == 1) rFlags[18] = 0x1;
613 if(chBlock[1] == 1) rFlags[17] = 0x1;
614 if(chBlock[2] == 1) rFlags[16] = 0x1;
617 rFlags[13] = puBits & 0x00000020;
618 rFlags[12] = puBits & 0x00000010;
619 rFlags[11] = puBits & 0x00000080;
620 rFlags[10] = puBits & 0x00000040;
621 rFlags[9] = puBits & 0x00000020;
622 rFlags[8] = puBits & 0x00000010;
624 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
625 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
626 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
627 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
628 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
629 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
631 UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
632 rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
633 0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
634 0x0 << 19 | rFlags[18] << 18 | rFlags[17] << 17 | rFlags[16] << 16 |
635 0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
636 rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
637 0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
638 rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
639 // --------------------------------------------------
641 // ****** Retrieving calibration data
642 // --- Equalization coefficients ---------------------------------------------
643 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
644 for(Int_t ji=0; ji<5; ji++){
645 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
646 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
647 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
648 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
650 // --- Energy calibration factors ------------------------------------
652 // **** Energy calibration coefficient set to 1
653 // **** (no trivial way to calibrate in p-p runs)
654 for(Int_t ij=0; ij<6; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
656 // ****** Equalization of detector responses
657 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
658 for(Int_t gi=0; gi<10; gi++){
660 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
661 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
662 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
663 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
666 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
667 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
668 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
669 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
673 /*printf("\n ------------- EQUALIZATION -------------\n");
674 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
675 equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
676 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
677 equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
678 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
679 equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
680 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
681 equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
682 printf(" ----------------------------------------\n");*/
684 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
685 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
686 for(Int_t gi=0; gi<5; gi++){
687 calibSumZN1[0] += equalTowZN1[gi];
688 calibSumZP1[0] += equalTowZP1[gi];
689 calibSumZN2[0] += equalTowZN2[gi];
690 calibSumZP2[0] += equalTowZP2[gi];
692 calibSumZN1[1] += equalTowZN1[gi+5];
693 calibSumZP1[1] += equalTowZP1[gi+5];
694 calibSumZN2[1] += equalTowZN2[gi+5];
695 calibSumZP2[1] += equalTowZP2[gi+5];
698 calibSumZN1[0] = calibSumZN1[0]*calibEne[0];
699 calibSumZP1[0] = calibSumZP1[0]*calibEne[1];
700 calibSumZN2[0] = calibSumZN2[0]*calibEne[2];
701 calibSumZP2[0] = calibSumZP2[0]*calibEne[3];
703 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
704 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
705 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
706 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
708 // ****** Energy calibration of detector responses
709 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
710 for(Int_t gi=0; gi<5; gi++){
712 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
713 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
714 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
715 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
717 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
718 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
719 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
720 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
723 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
724 calibZEM1[0] = corrADCZEM1[0]*calibEne[5];
725 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
726 calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
727 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
728 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
730 /*printf("\n ------------- CALIBRATION -------------\n");
731 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
732 calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
733 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
734 calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
735 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
736 calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
737 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
738 calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
739 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
740 printf(" ----------------------------------------\n");*/
742 // ****** No. of spectator and participants nucleons
743 // Variables calculated to comply with ESD structure
744 // *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
745 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
746 Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
747 Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
748 Double_t impPar=0., impPar1=0., impPar2=0.;
750 // create the output tree
751 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
752 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
753 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
754 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
755 nGenSpec, nGenSpecLeft, nGenSpecRight,
756 nPart, nPartTotLeft, nPartTotRight,
757 impPar, impPar1, impPar2,
758 recoFlag, isScalerOn, scaler);
760 const Int_t kBufferSize = 4000;
761 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
762 // write the output tree
763 clustersTree->Fill();
766 //_____________________________________________________________________________
767 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
768 Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
769 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
770 Bool_t isScalerOn, UInt_t* scaler,
771 Int_t* evQualityBlock, Int_t* triggerBlock, Int_t* chBlock, UInt_t puBits) const
773 // ****************** Reconstruct one event ******************
774 // ---------------------- Setting reco flags for ESD
776 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
778 if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
779 else rFlags[31] = 0x1;
781 if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
782 if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
783 if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
785 if(triggerBlock[0] == 1) rFlags[27] = 0x1;
786 if(triggerBlock[1] == 1) rFlags[26] = 0x1;
787 if(triggerBlock[2] == 1) rFlags[25] = 0x1;
788 if(triggerBlock[3] == 1) rFlags[24] = 0x1;
790 if(chBlock[0] == 1) rFlags[18] = 0x1;
791 if(chBlock[1] == 1) rFlags[17] = 0x1;
792 if(chBlock[2] == 1) rFlags[16] = 0x1;
794 rFlags[13] = puBits & 0x00000020;
795 rFlags[12] = puBits & 0x00000010;
796 rFlags[11] = puBits & 0x00000080;
797 rFlags[10] = puBits & 0x00000040;
798 rFlags[9] = puBits & 0x00000020;
799 rFlags[8] = puBits & 0x00000010;
801 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
802 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
803 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
804 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
805 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
806 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
808 UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
809 rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
810 0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
811 0x0 << 19 | rFlags[18] << 18 | rFlags[17] << 17 | rFlags[16] << 16 |
812 0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
813 rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
814 0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
815 rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
816 // --------------------------------------------------
819 // ****** Retrieving calibration data
820 // --- Equalization coefficients ---------------------------------------------
821 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
822 for(Int_t ji=0; ji<5; ji++){
823 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
824 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
825 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
826 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
828 // --- Energy calibration factors ------------------------------------
829 Float_t valFromOCDB[6], calibEne[6];
830 for(Int_t ij=0; ij<6; ij++){
831 valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
833 if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
834 else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
836 else calibEne[ij] = valFromOCDB[ij];
839 // ****** Equalization of detector responses
840 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
841 for(Int_t gi=0; gi<10; gi++){
843 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
844 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
845 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
846 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
849 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
850 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
851 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
852 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
856 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
857 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
858 for(Int_t gi=0; gi<5; gi++){
859 calibSumZN1[0] += equalTowZN1[gi];
860 calibSumZP1[0] += equalTowZP1[gi];
861 calibSumZN2[0] += equalTowZN2[gi];
862 calibSumZP2[0] += equalTowZP2[gi];
864 calibSumZN1[1] += equalTowZN1[gi+5];
865 calibSumZP1[1] += equalTowZP1[gi+5];
866 calibSumZN2[1] += equalTowZN2[gi+5];
867 calibSumZP2[1] += equalTowZP2[gi+5];
870 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]*8.;
871 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]*8.;
872 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]*8.;
873 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]*8.;
875 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
876 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
877 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
878 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
880 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
881 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]*8.;
882 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
883 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]*8.;
884 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
885 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
887 // ****** Energy calibration of detector responses
888 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
889 for(Int_t gi=0; gi<5; gi++){
891 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]*8.;
892 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]*8.;
893 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]*8.;
894 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]*8.;
896 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
897 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
898 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
899 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
902 // ****** Number of detected spectator nucleons
903 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
905 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
906 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
907 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
908 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
910 else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
911 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
912 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
913 nDetSpecNRight, nDetSpecPRight);*/
915 Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
916 Int_t nPart=0, nPartA=0, nPartC=0;
917 Double_t b=0., bA=0., bC=0.;
919 if(fIsCalibrationMB == kFALSE){
920 // ****** Reconstruction parameters ------------------
921 if (!fRecoParam) fRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
922 if(!fMBCalibData) fMBCalibData = const_cast<AliZDCMBCalib*>(GetMBCalibData());
924 TH2F *hZDCvsZEM = fMBCalibData->GethZDCvsZEM();
925 TH2F *hZDCCvsZEM = fMBCalibData->GethZDCCvsZEM();
926 TH2F *hZDCAvsZEM = fMBCalibData->GethZDCAvsZEM();
928 TH1D *hNpartDist = fRecoParam->GethNpartDist();
929 TH1D *hbDist = fRecoParam->GethbDist();
930 Float_t ClkCenter = fRecoParam->GetClkCenter();
932 Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
933 Double_t origin = xHighEdge*ClkCenter;
935 //printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
937 // ====> Summed ZDC info (sideA+side C)
938 TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
939 Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
940 Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
941 line->SetParameter(0, y/(x-origin));
942 line->SetParameter(1, -origin*y/(x-origin));
944 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
945 //printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
947 Double_t countPerc=0;
948 Double_t xBinCenter=0, yBinCenter=0;
949 for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
950 for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
951 xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
952 yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
954 if(line->GetParameter(0)>0){
955 if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
956 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
958 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
959 xBinCenter, yBinCenter, countPerc);*/
963 if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
964 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
966 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
967 xBinCenter, yBinCenter, countPerc);*/
973 Double_t xSecPerc = 0.;
974 if(hZDCvsZEM->GetEntries()!=0){
975 xSecPerc = countPerc/hZDCvsZEM->GetEntries();
978 AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
981 //printf(" xSecPerc %1.4f \n", xSecPerc);
984 TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
985 Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
986 lineC->SetParameter(0, yC/(x-origin));
987 lineC->SetParameter(1, -origin*yC/(x-origin));
989 //printf(" ***************** Side C \n");
990 //printf(" E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
992 Double_t countPercC=0;
993 Double_t xBinCenterC=0, yBinCenterC=0;
994 for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
995 for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
996 xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
997 yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
998 if(lineC->GetParameter(0)>0){
999 if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1000 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1004 if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1005 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1011 Double_t xSecPercC = 0.;
1012 if(hZDCCvsZEM->GetEntries()!=0){
1013 xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
1016 AliWarning(" Histogram hZDCCvsZEM from OCDB has no entries!!!");
1019 //printf(" xSecPercC %1.4f \n", xSecPercC);
1022 TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
1023 Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
1024 lineA->SetParameter(0, yA/(x-origin));
1025 lineA->SetParameter(1, -origin*yA/(x-origin));
1028 //printf(" ***************** Side A \n");
1029 //printf(" E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
1031 Double_t countPercA=0;
1032 Double_t xBinCenterA=0, yBinCenterA=0;
1033 for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
1034 for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
1035 xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
1036 yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
1037 if(lineA->GetParameter(0)>0){
1038 if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1039 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1043 if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1044 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1050 Double_t xSecPercA = 0.;
1051 if(hZDCAvsZEM->GetEntries()!=0){
1052 xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
1055 AliWarning(" Histogram hZDCAvsZEM from OCDB has no entries!!!");
1058 //printf(" xSecPercA %1.4f \n", xSecPercA);
1060 // ****** Number of participants (from E_ZDC vs. E_ZEM correlation)
1061 Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1062 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1063 nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1064 if((1.-nPartFrac) < xSecPerc){
1065 nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1067 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
1068 //printf(" nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1072 if(nPart<0) nPart=0;
1074 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1075 nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1076 if((1.-nPartFracC) < xSecPercC){
1077 nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1079 //printf(" ***************** Side C \n");
1080 //printf(" nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1084 if(nPartC<0) nPartC=0;
1086 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1087 nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1088 if((1.-nPartFracA) < xSecPercA){
1089 nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1091 //printf(" ***************** Side A \n");
1092 //printf(" nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1096 if(nPartA<0) nPartA=0;
1098 // ****** Impact parameter (from E_ZDC vs. E_ZEM correlation)
1099 Double_t bFrac=0., bFracC=0., bFracA=0.;
1100 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1101 bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1102 if((1.-bFrac) < xSecPerc){
1103 b = hbDist->GetBinLowEdge(ibbin);
1108 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1109 bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1110 if((1.-bFracC) < xSecPercC){
1111 bC = hbDist->GetBinLowEdge(ibbin);
1116 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1117 bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1118 if((1.-bFracA) < xSecPercA){
1119 bA = hbDist->GetBinLowEdge(ibbin);
1124 // ****** Number of spectator nucleons
1125 nGenSpec = 416 - nPart;
1126 nGenSpecC = 416 - nPartC;
1127 nGenSpecA = 416 - nPartA;
1128 if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1129 if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1130 if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1133 /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1134 " calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
1135 calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1136 printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n",
1137 nGenSpecLeft, nGenSpecRight);
1138 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1141 delete lineC; delete lineA;
1143 } // ONLY IF fIsCalibrationMB==kFALSE
1145 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1146 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1147 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1148 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1149 nGenSpec, nGenSpecA, nGenSpecC,
1150 nPart, nPartA, nPartC, b, bA, bC,
1151 recoFlag, isScalerOn, scaler);
1153 const Int_t kBufferSize = 4000;
1154 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1156 // write the output tree
1157 clustersTree->Fill();
1161 //_____________________________________________________________________________
1162 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1164 // fill energies and number of participants to the ESD
1167 AliZDCReco* preco = &reco;
1168 clustersTree->SetBranchAddress("ZDC", &preco);
1170 clustersTree->GetEntry(0);
1172 AliESDZDC * esdzdc = esd->GetESDZDC();
1173 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1174 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1175 for(Int_t i=0; i<5; i++){
1176 tZN1Ene[i] = reco.GetZN1HREnTow(i);
1177 tZN2Ene[i] = reco.GetZN2HREnTow(i);
1178 tZP1Ene[i] = reco.GetZP1HREnTow(i);
1179 tZP2Ene[i] = reco.GetZP2HREnTow(i);
1181 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1182 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1183 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1184 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1187 esdzdc->SetZN1TowerEnergy(tZN1Ene);
1188 esdzdc->SetZN2TowerEnergy(tZN2Ene);
1189 esdzdc->SetZP1TowerEnergy(tZP1Ene);
1190 esdzdc->SetZP2TowerEnergy(tZP2Ene);
1192 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1193 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1194 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1195 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1197 Int_t nPart = reco.GetNParticipants();
1198 Int_t nPartA = reco.GetNPartSideA();
1199 Int_t nPartC = reco.GetNPartSideC();
1200 Double_t b = reco.GetImpParameter();
1201 Double_t bA = reco.GetImpParSideA();
1202 Double_t bC = reco.GetImpParSideC();
1203 UInt_t recoFlag = reco.GetRecoFlag();
1205 esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
1206 reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1207 nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1209 // Writing ZDC scaler for cross section calculation
1210 // ONLY IF the scaler has been read during the event
1211 if(reco.IsScalerOn()==kTRUE){
1213 for(Int_t jk=0; jk<32; jk++) counts[jk] = reco.GetZDCScaler(jk);
1214 esd->SetZDCScaler(counts);
1218 //_____________________________________________________________________________
1219 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1221 // Setting the storage
1223 Bool_t deleteManager = kFALSE;
1225 AliCDBManager *manager = AliCDBManager::Instance();
1226 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1228 if(!defstorage || !(defstorage->Contains("ZDC"))){
1229 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1230 manager->SetDefaultStorage(uri);
1231 deleteManager = kTRUE;
1234 AliCDBStorage *storage = manager->GetDefaultStorage();
1237 AliCDBManager::Instance()->UnsetDefaultStorage();
1238 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1244 //_____________________________________________________________________________
1245 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1248 // Getting pedestal calibration object for ZDC set
1250 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1251 if(!entry) AliFatal("No calibration data loaded!");
1252 entry->SetOwner(kFALSE);
1254 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
1255 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1260 //_____________________________________________________________________________
1261 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1264 // Getting energy and equalization calibration object for ZDC set
1266 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1267 if(!entry) AliFatal("No calibration data loaded!");
1268 entry->SetOwner(kFALSE);
1270 AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1271 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1276 //_____________________________________________________________________________
1277 AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1280 // Getting energy and equalization calibration object for ZDC set
1282 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1283 if(!entry) AliFatal("No calibration data loaded!");
1284 entry->SetOwner(kFALSE);
1286 AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1287 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1292 //_____________________________________________________________________________
1293 AliZDCMBCalib* AliZDCReconstructor::GetMBCalibData() const
1296 // Getting energy and equalization calibration object for ZDC set
1298 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/MBCalib");
1299 if(!entry) AliFatal("No calibration data loaded!");
1300 entry->SetOwner(kFALSE);
1302 AliZDCMBCalib *calibdata = dynamic_cast<AliZDCMBCalib*> (entry->GetObject());
1303 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1308 //_____________________________________________________________________________
1309 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1312 // Getting reconstruction parameters from OCDB
1314 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1315 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1317 AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1318 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1324 //_____________________________________________________________________________
1325 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1328 // Getting reconstruction parameters from OCDB
1330 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1331 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1333 AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1334 if(!param) AliFatal("No RecoParam object in OCDB entry!");