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 "AliZDCRecoParam.h"
46 #include "AliZDCRecoParampp.h"
47 #include "AliZDCRecoParamPbPb.h"
48 #include "AliRunInfo.h"
51 ClassImp(AliZDCReconstructor)
52 AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0; //reconstruction parameters
54 //_____________________________________________________________________________
55 AliZDCReconstructor:: AliZDCReconstructor() :
56 fPedData(GetPedestalData()),
57 fEnCalibData(GetEnergyCalibData()),
58 fTowCalibData(GetTowerCalibData()),
62 fIsCalibrationMB(kFALSE),
66 // **** Default constructor
70 //_____________________________________________________________________________
71 AliZDCReconstructor::~AliZDCReconstructor()
74 // if(fRecoParam) delete fRecoParam;
75 if(fPedData) delete fPedData;
76 if(fEnCalibData) delete fEnCalibData;
77 if(fTowCalibData) delete fTowCalibData;
80 //____________________________________________________________________________
81 void AliZDCReconstructor::Init()
83 // Setting reconstruction mode
84 // Getting beam type and beam energy from GRP calibration object
86 TString runType = GetRunInfo()->GetRunType();
87 if((runType.CompareTo("CALIBRATION_MB")) == 0){
88 fIsCalibrationMB = kTRUE;
91 TString beamType = GetRunInfo()->GetBeamType();
92 // This is a temporary solution to allow reconstruction in tests without beam
93 if(((beamType.CompareTo("UNKNOWN"))==0) && ((runType.CompareTo("PHYSICS")) == 0)){
96 /*else if((beamType.CompareTo("UNKNOWN"))==0){
97 AliError("\t UNKNOWN beam type\n");
101 if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
102 ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)){
105 else if((beamType.CompareTo("A-A")) == 0){
109 fBeamEnergy = GetRunInfo()->GetBeamEnergy();
111 AliWarning(" Beam energy value missing -> E_beam = 0");
115 if(fIsCalibrationMB==kFALSE)
116 printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f GeV *****\n\n",beamType.Data(), fBeamEnergy);
120 //_____________________________________________________________________________
121 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
123 // *** Local ZDC reconstruction for digits
124 // Works on the current event
126 // Retrieving calibration data
127 // Parameters for mean value pedestal subtraction
129 Float_t meanPed[2*kNch];
130 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
131 // Parameters pedestal subtraction through correlation with out-of-time signals
132 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
133 for(Int_t jj=0; jj<2*kNch; jj++){
134 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
135 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
140 AliZDCDigit* pdigit = &digit;
141 digitsTree->SetBranchAddress("ZDC", &pdigit);
142 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
145 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
146 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
147 for(Int_t i=0; i<10; i++){
148 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
149 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
152 Int_t digNentries = digitsTree->GetEntries();
153 Float_t ootDigi[kNch];
154 // -- Reading out-of-time signals (last kNch entries) for current event
156 for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
157 ootDigi[iDigit] = digitsTree->GetEntry(iDigit);
161 for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
162 digitsTree->GetEntry(iDigit);
163 if (!pdigit) continue;
165 Int_t det = digit.GetSector(0);
166 Int_t quad = digit.GetSector(1);
168 Float_t ped2SubHg=0., ped2SubLg=0.;
170 if(det==1) pedindex = quad;
171 else if(det==2) pedindex = quad+5;
172 else if(det==3) pedindex = quad+9;
173 else if(det==4) pedindex = quad+12;
174 else if(det==5) pedindex = quad+17;
176 else pedindex = (det-1)/3+22;
179 ped2SubHg = meanPed[pedindex];
180 ped2SubLg = meanPed[pedindex+kNch];
182 else if(fPedSubMode==1){
183 ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
184 ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
187 if(quad != 5){ // ZDC (not reference PTMs!)
188 if(det == 1){ // *** ZNC
189 tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
190 tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
191 if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
192 if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
194 else if(det == 2){ // *** ZP1
195 tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
196 tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
197 if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
198 if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
201 if(quad == 1){ // *** ZEM1
202 dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
203 dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
204 if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
205 if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
207 else if(quad == 2){ // *** ZEM2
208 dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
209 dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
210 if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
211 if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
214 else if(det == 4){ // *** ZN2
215 tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
216 tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
217 if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
218 if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
220 else if(det == 5){ // *** ZP2
221 tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
222 tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
223 if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
224 if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
227 else{ // Reference PMs
229 sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
230 sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
232 if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
233 if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
236 sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
237 sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
239 if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
240 if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
245 /*printf("AliZDCReconstructor: digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
246 iDigit, det, quad, ped2SubHg, ped2SubLg);
247 printf(" -> pedindex %d\n", pedindex);
248 printf(" HGChain -> RawDig %d DigCorr %1.2f",
249 digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg);
250 printf(" LGChain -> RawDig %d DigCorr %1.2f\n",
251 digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);*/
256 for(Int_t jj=0; jj<32; jj++) counts[jj]=0;
258 Int_t evQualityBlock[4] = {1,0,0,0};
259 Int_t triggerBlock[4] = {0,0,0,0};
260 Int_t chBlock[3] = {0,0,0};
263 // reconstruct the event
265 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
266 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
268 evQualityBlock, triggerBlock, chBlock, puBits);
269 else if(fRecoMode==2)
270 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
271 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
273 evQualityBlock, triggerBlock, chBlock, puBits);
276 //_____________________________________________________________________________
277 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
279 // *** ZDC raw data reconstruction
280 // Works on the current event
282 // Retrieving calibration data
283 // Parameters for pedestal subtraction
285 Float_t meanPed[2*kNch];
286 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
287 // Parameters pedestal subtraction through correlation with out-of-time signals
288 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
289 for(Int_t jj=0; jj<2*kNch; jj++){
290 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
291 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
294 Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
295 Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
296 Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
297 Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
298 Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
299 Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
300 for(Int_t ich=0; ich<5; ich++){
301 adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
302 adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
303 adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
304 adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
306 adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
307 pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
311 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
312 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
313 for(Int_t i=0; i<10; i++){
314 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
315 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
318 Bool_t isScalerOn=kFALSE;
320 UInt_t scalerData[32];
321 for(Int_t k=0; k<32; k++) scalerData[k]=0;
323 Int_t evQualityBlock[4] = {1,0,0,0};
324 Int_t triggerBlock[4] = {0,0,0,0};
325 Int_t chBlock[3] = {0,0,0};
328 //fNRun = (Int_t) rawReader->GetRunNumber();
329 Int_t kFirstADCGeo=0, kLastADCGeo=3, kScalerGeo=8, kPUGeo=29;
330 //Int_t kTrigScales=30, kTrigHistory=31;
332 // loop over raw data
334 AliZDCRawStream rawData(rawReader);
335 while(rawData.Next()){
337 // ***************************** Reading ADCs
338 if((rawData.GetADCModule()>=kFirstADCGeo) && (rawData.GetADCModule()<=kLastADCGeo)){
339 //printf(" **** Reading ADC raw data from module %d **** \n",rawData.GetADCModule());
341 if((rawData.IsADCDataWord()) && (rawData.GetNChannelsOn()<48)) chBlock[0] = kTRUE;
342 if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)) chBlock[1] = kTRUE;
343 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)) chBlock[2] = kTRUE;
344 if((rawData.IsADCDataWord()) && (rawData.IsADCEventGood() == kTRUE)) evQualityBlock[0] = kTRUE;
346 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow()==kFALSE)
347 && (rawData.IsOverflow()==kFALSE) && (rawData.IsADCEventGood()==kTRUE)){
349 Int_t adcMod = rawData.GetADCModule();
350 Int_t det = rawData.GetSector(0);
351 Int_t quad = rawData.GetSector(1);
352 Int_t gain = rawData.GetADCGain();
355 // Mean pedestal value subtraction -------------------------------------------------------
356 if(fPedSubMode == 0){
357 // Not interested in o.o.t. signals (ADC modules 2, 3)
358 if(adcMod == 2 || adcMod == 3) continue;
360 if(quad != 5){ // ZDCs (not reference PTMs)
363 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
364 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
368 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
369 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
374 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
375 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
378 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
379 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
384 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
385 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
389 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
390 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
393 else{ // reference PM
394 pedindex = (det-1)/3 + 22;
396 if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
397 else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
400 if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
401 else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
406 printf(" AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f",
407 det,quad,gain, pedindex, meanPed[pedindex]);
408 printf(" RawADC %d ADCCorr %1.0f\n",
409 rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);
411 }// mean pedestal subtraction
412 // Pedestal subtraction from correlation ------------------------------------------------
413 else if(fPedSubMode == 1){
415 if(adcMod==0 || adcMod==1){
416 if(quad != 5){ // signals from ZDCs
418 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
419 else adcZN1lg[quad] = rawData.GetADCValue();
422 if(gain==0) adcZP1[quad] = rawData.GetADCValue();
423 else adcZP1lg[quad] = rawData.GetADCValue();
426 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
427 else adcZEMlg[quad-1] = rawData.GetADCValue();
430 if(gain==0) adcZN2[quad] = rawData.GetADCValue();
431 else adcZN2lg[quad] = rawData.GetADCValue();
434 if(gain==0) adcZP2[quad] = rawData.GetADCValue();
435 else adcZP2lg[quad] = rawData.GetADCValue();
438 else{ // signals from reference PM
439 if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
440 else pmReflg[quad-1] = rawData.GetADCValue();
443 // Out-of-time pedestals
444 else if(adcMod==2 || adcMod==3){
445 if(quad != 5){ // signals from ZDCs
447 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
448 else adcZN1ootlg[quad] = rawData.GetADCValue();
451 if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
452 else adcZP1ootlg[quad] = rawData.GetADCValue();
455 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
456 else adcZEMootlg[quad-1] = rawData.GetADCValue();
459 if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
460 else adcZN2ootlg[quad] = rawData.GetADCValue();
463 if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
464 else adcZP2ootlg[quad] = rawData.GetADCValue();
467 else{ // signals from reference PM
468 if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
469 else pmRefootlg[quad-1] = rawData.GetADCValue();
472 } // pedestal subtraction from correlation
474 /*printf("\t AliZDCReconstructor: det %d quad %d res %d -> Ped[%d] = %1.0f\n",
475 det,quad,gain, pedindex, meanPed[pedindex]);*/
478 // ***************************** Reading Scaler
479 else if(rawData.GetADCModule()==kScalerGeo){
480 if(rawData.IsScHeaderRead()==kTRUE && rawData.IsScEventGood()==kTRUE){
482 scalerData[jsc] = rawData.GetTriggerCount();
486 // ***************************** Reading PU
487 else if(rawData.GetADCModule()==kPUGeo){
488 puBits = rawData.GetDetectorPattern();
490 // ***************************** Reading trigger history
491 else if(rawData.IstriggerHistoryWord()==kTRUE){
492 triggerBlock[0] = rawData.IsCPTInputEMDTrigger();
493 triggerBlock[1] = rawData.IsCPTInputSemiCentralTrigger();
494 triggerBlock[2] = rawData.IsCPTInputCentralTrigger();
495 triggerBlock[3] = rawData.IsCPTInputMBTrigger();
501 for(Int_t t=0; t<5; t++){
502 tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
503 tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
505 tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
506 tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
508 tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
509 tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
511 tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
512 tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
513 // 0---------0 Ch. debug 0---------0
514 /* printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
515 printf("\tCorrCoeff0\tCorrCoeff1\n");
516 printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
517 printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
518 printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
519 printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
520 printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
521 printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
522 printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
523 printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
525 printf("ZN1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
526 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
527 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
528 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
530 printf("ZP1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
531 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
532 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
533 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
535 printf("ZN2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
536 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
537 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
538 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
540 printf("ZP2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
541 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
542 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
543 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
546 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
547 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
548 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
549 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
551 sPMRef1[0] = pmRef[0] - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
552 sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
553 sPMRef2[0] = pmRef[0] - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
554 sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
557 if(fRecoMode==1) // p-p data
558 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
559 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
560 isScalerOn, scalerData,
561 evQualityBlock, triggerBlock, chBlock, puBits);
562 else if(fRecoMode==2) // Pb-Pb data
563 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
564 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
565 isScalerOn, scalerData,
566 evQualityBlock, triggerBlock, chBlock, puBits);
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 isScalerOn, UInt_t* scaler,
574 Int_t* evQualityBlock, Int_t* triggerBlock, Int_t* chBlock, UInt_t puBits) const
576 // ****************** Reconstruct one event ******************
579 /*printf("\n*************************************************\n");
580 printf(" ReconstructEventpp -> values after pedestal subtraction:\n");
581 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
582 corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
583 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
584 corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
585 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
586 corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
587 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
588 corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
589 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
590 printf("*************************************************\n");*/
592 // ---------------------- Setting reco flags for ESD
594 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
596 if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
597 else rFlags[31] = 0x1;
599 if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
600 if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
601 if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
603 if(triggerBlock[0] == 1) rFlags[27] = 0x1;
604 if(triggerBlock[1] == 1) rFlags[26] = 0x1;
605 if(triggerBlock[2] == 1) rFlags[25] = 0x1;
606 if(triggerBlock[3] == 1) rFlags[24] = 0x1;
608 if(chBlock[0] == 1) rFlags[18] = 0x1;
609 if(chBlock[1] == 1) rFlags[17] = 0x1;
610 if(chBlock[2] == 1) rFlags[16] = 0x1;
613 rFlags[13] = puBits & 0x00000020;
614 rFlags[12] = puBits & 0x00000010;
615 rFlags[11] = puBits & 0x00000080;
616 rFlags[10] = puBits & 0x00000040;
617 rFlags[9] = puBits & 0x00000020;
618 rFlags[8] = puBits & 0x00000010;
620 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
621 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
622 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
623 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
624 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
625 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
627 UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
628 rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
629 0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
630 0x0 << 19 | rFlags[18] << 18 | rFlags[17] << 17 | rFlags[16] << 16 |
631 0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
632 rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
633 0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
634 rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
635 // --------------------------------------------------
637 // ****** Retrieving calibration data
638 // --- Equalization coefficients ---------------------------------------------
639 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
640 for(Int_t ji=0; ji<5; ji++){
641 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
642 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
643 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
644 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
646 // --- Energy calibration factors ------------------------------------
648 // **** Energy calibration coefficient set to 1
649 // **** (no trivial way to calibrate in p-p runs)
650 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
652 // ****** Equalization of detector responses
653 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
654 for(Int_t gi=0; gi<10; gi++){
656 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
657 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
658 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
659 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
662 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
663 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
664 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
665 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
669 /*printf("\n ------------- EQUALIZATION -------------\n");
670 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
671 equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
672 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
673 equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
674 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
675 equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
676 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
677 equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
678 printf(" ----------------------------------------\n");*/
680 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
681 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
682 for(Int_t gi=0; gi<5; gi++){
683 calibSumZN1[0] += equalTowZN1[gi];
684 calibSumZP1[0] += equalTowZP1[gi];
685 calibSumZN2[0] += equalTowZN2[gi];
686 calibSumZP2[0] += equalTowZP2[gi];
688 calibSumZN1[1] += equalTowZN1[gi+5];
689 calibSumZP1[1] += equalTowZP1[gi+5];
690 calibSumZN2[1] += equalTowZN2[gi+5];
691 calibSumZP2[1] += equalTowZP2[gi+5];
694 calibSumZN1[0] = calibSumZN1[0]*calibEne[0];
695 calibSumZP1[0] = calibSumZP1[0]*calibEne[1];
696 calibSumZN2[0] = calibSumZN2[0]*calibEne[2];
697 calibSumZP2[0] = calibSumZP2[0]*calibEne[3];
699 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
700 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
701 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
702 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
704 // ****** Energy calibration of detector responses
705 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
706 for(Int_t gi=0; gi<5; gi++){
708 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
709 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
710 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
711 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
713 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
714 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
715 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
716 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
719 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
720 calibZEM1[0] = corrADCZEM1[0]*calibEne[5];
721 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
722 calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
723 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
724 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
726 /*printf("\n ------------- CALIBRATION -------------\n");
727 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
728 calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
729 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
730 calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
731 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
732 calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
733 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
734 calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
735 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
736 printf(" ----------------------------------------\n");*/
738 // ****** No. of spectator and participants nucleons
739 // Variables calculated to comply with ESD structure
740 // *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
741 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
742 Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
743 Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
744 Double_t impPar=0., impPar1=0., impPar2=0.;
746 // create the output tree
747 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
748 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
749 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
750 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
751 nGenSpec, nGenSpecLeft, nGenSpecRight,
752 nPart, nPartTotLeft, nPartTotRight,
753 impPar, impPar1, impPar2,
754 recoFlag, isScalerOn, scaler);
756 const Int_t kBufferSize = 4000;
757 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
758 // write the output tree
759 clustersTree->Fill();
762 //_____________________________________________________________________________
763 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
764 Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
765 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
766 Bool_t isScalerOn, UInt_t* scaler,
767 Int_t* evQualityBlock, Int_t* triggerBlock, Int_t* chBlock, UInt_t puBits) const
769 // ****************** Reconstruct one event ******************
770 // ---------------------- Setting reco flags for ESD
772 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
774 if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
775 else rFlags[31] = 0x1;
777 if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
778 if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
779 if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
781 if(triggerBlock[0] == 1) rFlags[27] = 0x1;
782 if(triggerBlock[1] == 1) rFlags[26] = 0x1;
783 if(triggerBlock[2] == 1) rFlags[25] = 0x1;
784 if(triggerBlock[3] == 1) rFlags[24] = 0x1;
786 if(chBlock[0] == 1) rFlags[18] = 0x1;
787 if(chBlock[1] == 1) rFlags[17] = 0x1;
788 if(chBlock[2] == 1) rFlags[16] = 0x1;
790 rFlags[13] = puBits & 0x00000020;
791 rFlags[12] = puBits & 0x00000010;
792 rFlags[11] = puBits & 0x00000080;
793 rFlags[10] = puBits & 0x00000040;
794 rFlags[9] = puBits & 0x00000020;
795 rFlags[8] = puBits & 0x00000010;
797 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
798 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
799 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
800 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
801 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
802 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
804 UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
805 rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
806 0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
807 0x0 << 19 | rFlags[18] << 18 | rFlags[17] << 17 | rFlags[16] << 16 |
808 0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
809 rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
810 0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
811 rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
812 // --------------------------------------------------
815 // ****** Retrieving calibration data
816 // --- Equalization coefficients ---------------------------------------------
817 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
818 for(Int_t ji=0; ji<5; ji++){
819 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
820 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
821 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
822 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
824 // --- Energy calibration factors ------------------------------------
825 Float_t valFromOCDB[6], calibEne[6];
826 for(Int_t ij=0; ij<6; ij++){
827 valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
829 if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
830 else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
832 else calibEne[ij] = valFromOCDB[ij];
835 // ****** Equalization of detector responses
836 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
837 for(Int_t gi=0; gi<10; gi++){
839 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
840 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
841 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
842 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
845 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
846 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
847 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
848 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
852 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
853 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
854 for(Int_t gi=0; gi<5; gi++){
855 calibSumZN1[0] += equalTowZN1[gi];
856 calibSumZP1[0] += equalTowZP1[gi];
857 calibSumZN2[0] += equalTowZN2[gi];
858 calibSumZP2[0] += equalTowZP2[gi];
860 calibSumZN1[1] += equalTowZN1[gi+5];
861 calibSumZP1[1] += equalTowZP1[gi+5];
862 calibSumZN2[1] += equalTowZN2[gi+5];
863 calibSumZP2[1] += equalTowZP2[gi+5];
866 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]*8.;
867 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]*8.;
868 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]*8.;
869 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]*8.;
871 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
872 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
873 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
874 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
876 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
877 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]*8.;
878 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
879 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]*8.;
880 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
881 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
883 // ****** Energy calibration of detector responses
884 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
885 for(Int_t gi=0; gi<5; gi++){
887 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]*8.;
888 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]*8.;
889 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]*8.;
890 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]*8.;
892 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
893 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
894 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
895 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
898 // ****** Number of detected spectator nucleons
899 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
901 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
902 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
903 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
904 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
906 else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
907 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
908 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
909 nDetSpecNRight, nDetSpecPRight);*/
911 Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
912 Int_t nPart=0, nPartA=0, nPartC=0;
913 Double_t b=0., bA=0., bC=0.;
915 if(fIsCalibrationMB == kFALSE){
916 // ****** Reconstruction parameters ------------------
918 //fRecoParam->Print("");
921 if (!fRecoParam) fRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
923 TH2F *hZDCvsZEM = fRecoParam->GethZDCvsZEM();
924 TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
925 TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
926 TH1D *hNpartDist = fRecoParam->GethNpartDist();
927 TH1D *hbDist = fRecoParam->GethbDist();
928 Float_t ClkCenter = fRecoParam->GetClkCenter();
930 Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
931 Double_t origin = xHighEdge*ClkCenter;
933 //printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
935 // ====> Summed ZDC info (sideA+side C)
936 TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
937 Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
938 Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
939 line->SetParameter(0, y/(x-origin));
940 line->SetParameter(1, -origin*y/(x-origin));
942 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
943 //printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
945 Double_t countPerc=0;
946 Double_t xBinCenter=0, yBinCenter=0;
947 for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
948 for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
949 xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
950 yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
952 if(line->GetParameter(0)>0){
953 if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
954 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
956 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
957 xBinCenter, yBinCenter, countPerc);*/
961 if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
962 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
964 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
965 xBinCenter, yBinCenter, countPerc);*/
971 Double_t xSecPerc = 0.;
972 if(hZDCvsZEM->GetEntries()!=0){
973 xSecPerc = countPerc/hZDCvsZEM->GetEntries();
976 AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
979 //printf(" xSecPerc %1.4f \n", xSecPerc);
982 TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
983 Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
984 lineC->SetParameter(0, yC/(x-origin));
985 lineC->SetParameter(1, -origin*yC/(x-origin));
987 //printf(" ***************** Side C \n");
988 //printf(" E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
990 Double_t countPercC=0;
991 Double_t xBinCenterC=0, yBinCenterC=0;
992 for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
993 for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
994 xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
995 yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
996 if(lineC->GetParameter(0)>0){
997 if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
998 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1002 if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1003 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1009 Double_t xSecPercC = 0.;
1010 if(hZDCCvsZEM->GetEntries()!=0){
1011 xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
1014 AliWarning(" Histogram hZDCCvsZEM from OCDB has no entries!!!");
1017 //printf(" xSecPercC %1.4f \n", xSecPercC);
1020 TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
1021 Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
1022 lineA->SetParameter(0, yA/(x-origin));
1023 lineA->SetParameter(1, -origin*yA/(x-origin));
1026 //printf(" ***************** Side A \n");
1027 //printf(" E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
1029 Double_t countPercA=0;
1030 Double_t xBinCenterA=0, yBinCenterA=0;
1031 for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
1032 for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
1033 xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
1034 yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
1035 if(lineA->GetParameter(0)>0){
1036 if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1037 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1041 if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1042 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1048 Double_t xSecPercA = 0.;
1049 if(hZDCAvsZEM->GetEntries()!=0){
1050 xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
1053 AliWarning(" Histogram hZDCAvsZEM from OCDB has no entries!!!");
1056 //printf(" xSecPercA %1.4f \n", xSecPercA);
1058 // ****** Number of participants (from E_ZDC vs. E_ZEM correlation)
1059 Int_t nPart=0, nPartC=0, nPartA=0;
1060 Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1061 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1062 nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1063 if((1.-nPartFrac) < xSecPerc){
1064 nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1066 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
1067 //printf(" nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1071 if(nPart<0) nPart=0;
1073 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1074 nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1075 if((1.-nPartFracC) < xSecPercC){
1076 nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1078 //printf(" ***************** Side C \n");
1079 //printf(" nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1083 if(nPartC<0) nPartC=0;
1085 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1086 nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1087 if((1.-nPartFracA) < xSecPercA){
1088 nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1090 //printf(" ***************** Side A \n");
1091 //printf(" nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1095 if(nPartA<0) nPartA=0;
1097 // ****** Impact parameter (from E_ZDC vs. E_ZEM correlation)
1098 Float_t b=0, bC=0, bA=0;
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 Int_t nGenSpec=0, nGenSpecC=0, nGenSpecA=0;
1127 nGenSpec = 416 - nPart;
1128 nGenSpecC = 416 - nPartC;
1129 nGenSpecA = 416 - nPartA;
1130 if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1131 if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1132 if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1135 /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1136 " calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
1137 calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1138 printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n",
1139 nGenSpecLeft, nGenSpecRight);
1140 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1143 delete lineC; delete lineA;
1145 } // ONLY IF fIsCalibrationMB==kFALSE
1147 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1148 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1149 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1150 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1151 nGenSpec, nGenSpecA, nGenSpecC,
1152 nPart, nPartA, nPartC, b, bA, bC,
1153 recoFlag, isScalerOn, scaler);
1155 const Int_t kBufferSize = 4000;
1156 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1158 // write the output tree
1159 clustersTree->Fill();
1162 //_____________________________________________________________________________
1163 void AliZDCReconstructor::BuildRecoParam(Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1165 // Calculate RecoParam object for Pb-Pb data
1166 (fRecoParam->GethZDCvsZEM())->Fill(ZDCC+ZDCA, ZEM);
1167 (fRecoParam->GethZDCCvsZEM())->Fill(ZDCC, ZEM);
1168 (fRecoParam->GethZDCAvsZEM())->Fill(ZDCA, ZEM);
1172 //_____________________________________________________________________________
1173 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1175 // fill energies and number of participants to the ESD
1178 AliZDCReco* preco = &reco;
1179 clustersTree->SetBranchAddress("ZDC", &preco);
1181 clustersTree->GetEntry(0);
1183 AliESDZDC * esdzdc = esd->GetESDZDC();
1184 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1185 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1186 for(Int_t i=0; i<5; i++){
1187 tZN1Ene[i] = reco.GetZN1HREnTow(i);
1188 tZN2Ene[i] = reco.GetZN2HREnTow(i);
1189 tZP1Ene[i] = reco.GetZP1HREnTow(i);
1190 tZP2Ene[i] = reco.GetZP2HREnTow(i);
1192 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1193 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1194 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1195 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1198 esdzdc->SetZN1TowerEnergy(tZN1Ene);
1199 esdzdc->SetZN2TowerEnergy(tZN2Ene);
1200 esdzdc->SetZP1TowerEnergy(tZP1Ene);
1201 esdzdc->SetZP2TowerEnergy(tZP2Ene);
1203 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1204 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1205 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1206 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1208 Int_t nPart = reco.GetNParticipants();
1209 Int_t nPartA = reco.GetNPartSideA();
1210 Int_t nPartC = reco.GetNPartSideC();
1211 Double_t b = reco.GetImpParameter();
1212 Double_t bA = reco.GetImpParSideA();
1213 Double_t bC = reco.GetImpParSideC();
1214 UInt_t recoFlag = reco.GetRecoFlag();
1216 esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
1217 reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1218 nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1220 // Writing ZDC scaler for cross section calculation
1221 // ONLY IF the scaler has been read during the event
1222 if(reco.IsScalerOn()==kTRUE){
1224 for(Int_t jk=0; jk<32; jk++) counts[jk] = reco.GetZDCScaler(jk);
1225 esd->SetZDCScaler(counts);
1229 //_____________________________________________________________________________
1230 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1232 // Setting the storage
1234 Bool_t deleteManager = kFALSE;
1236 AliCDBManager *manager = AliCDBManager::Instance();
1237 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1239 if(!defstorage || !(defstorage->Contains("ZDC"))){
1240 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1241 manager->SetDefaultStorage(uri);
1242 deleteManager = kTRUE;
1245 AliCDBStorage *storage = manager->GetDefaultStorage();
1248 AliCDBManager::Instance()->UnsetDefaultStorage();
1249 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1255 //_____________________________________________________________________________
1256 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1259 // Getting pedestal calibration object for ZDC set
1261 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1262 if(!entry) AliFatal("No calibration data loaded!");
1264 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
1265 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1270 //_____________________________________________________________________________
1271 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1274 // Getting energy and equalization calibration object for ZDC set
1276 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1277 if(!entry) AliFatal("No calibration data loaded!");
1279 AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1280 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1285 //_____________________________________________________________________________
1286 AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1289 // Getting energy and equalization calibration object for ZDC set
1291 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1292 if(!entry) AliFatal("No calibration data loaded!");
1294 AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1295 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1300 //_____________________________________________________________________________
1301 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1304 // Getting reconstruction parameters from OCDB
1306 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1307 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1309 AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1310 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1316 //_____________________________________________________________________________
1317 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1320 // Getting reconstruction parameters from OCDB
1322 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1323 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1325 AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1326 if(!param) AliFatal("No RecoParam object in OCDB entry!");