]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ZDC/AliZDCReconstructor.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
RecoParam object provided to comply with standars reconstruction
[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// ************** Class for ZDC reconstruction ************** //
21// Author: Chiara.Oppedisano@to.infn.it //
22// //
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) //
26// //
27///////////////////////////////////////////////////////////////////////////////
28
29
30#include <TH2F.h>
31#include <TH1D.h>
32#include <TAxis.h>
33#include <TMap.h>
34
35#include "AliRawReader.h"
36#include "AliGRPObject.h"
37#include "AliESDEvent.h"
38#include "AliESDZDC.h"
39#include "AliZDCDigit.h"
40#include "AliZDCRawStream.h"
41#include "AliZDCReco.h"
42#include "AliZDCReconstructor.h"
43#include "AliZDCPedestals.h"
44#include "AliZDCEnCalib.h"
45#include "AliZDCTowerCalib.h"
46#include "AliZDCRecoParam.h"
47#include "AliZDCRecoParampp.h"
48#include "AliZDCRecoParamPbPb.h"
49
50
51ClassImp(AliZDCReconstructor)
52AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0; //reconstruction parameters
53
54//_____________________________________________________________________________
55AliZDCReconstructor:: AliZDCReconstructor() :
56 fPedData(GetPedestalData()),
57 fEnCalibData(GetEnergyCalibData()),
58 fTowCalibData(GetTowerCalibData()),
59 fRecoMode(0),
60 fBeamEnergy(0.),
61 fNRun(0),
62 fIsCalibrationMB(kFALSE),
63 fPedSubMode(0),
64 fSignalThreshold(7)
65{
66 // **** Default constructor
67
68}
69
70
71//_____________________________________________________________________________
72AliZDCReconstructor::~AliZDCReconstructor()
73{
74// destructor
75 if(fRecoParam) delete fRecoParam;
76 if(fPedData) delete fPedData;
77 if(fEnCalibData) delete fEnCalibData;
78 if(fTowCalibData) delete fTowCalibData;
79}
80
81//____________________________________________________________________________
82void AliZDCReconstructor::Init()
83{
84 // Setting reconstruction mode
85 // Getting beam type and beam energy from GRP calibration object
86
87 if(fRecoMode==0 && fBeamEnergy==0.){
88 // Initialization of the GRP entry
89 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
90 AliGRPObject* grpData = 0x0;
91 if(entry){
92 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
93 if(m){
94 //m->Print();
95 grpData = new AliGRPObject();
96 grpData->ReadValuesFromMap(m);
97 }
98 else{
99 grpData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
100 }
101 entry->SetOwner(0);
102 AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
103 }
104 if(!grpData) AliError("No GRP entry found in OCDB!");
105
106 TString runType = grpData->GetRunType();
107 if(runType==AliGRPObject::GetInvalidString()){
108 AliWarning("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN");
109 runType = "UNKNOWN";
110 }
111 if((runType.CompareTo("CALIBRATION_MB")) == 0){
112 fIsCalibrationMB = kTRUE;
113 }
114
115 TString beamType = grpData->GetBeamType();
116 if(beamType==AliGRPObject::GetInvalidString()){
117 AliWarning("GRP/GRP/Data entry: missing value for the beam energy !");
118 AliError("\t ZDC does not reconstruct event 4 UNKNOWN beam type\n");
119 return;
120 }
121 if((beamType.CompareTo("P-P")) == 0){
122 fRecoMode=1;
123 fRecoParam = (AliZDCRecoParampp*) GetppRecoParamFromOCDB();
124 AliInfo(" Getting AliZDCRecoParampp object from OCDB \n");
125 }
126 else if((beamType.CompareTo("A-A")) == 0){
127 fRecoMode=2;
128 if(fIsCalibrationMB == kFALSE){
129 fRecoParam = (AliZDCRecoParamPbPb*) GetPbPbRecoParamFromOCDB();
130 AliInfo(" Getting AliZDCRecoParamPbPb object from OCDB\n");
131 }
132 else{
133 fRecoParam = new AliZDCRecoParamPbPb();
134 //
135 TH2F* hZDCvsZEM = new TH2F("hZDCvsZEM","hZDCvsZEM",100,0.,10.,100,0.,1000.);
136 hZDCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCvsZEM->SetYTitle("E_{ZDC} (TeV)");
137 fRecoParam->SetZDCvsZEM(hZDCvsZEM);
138 //
139 TH2F* hZDCCvsZEM = new TH2F("hZDCCvsZEM","hZDCCvsZEM",100,0.,10.,100,0.,500.);
140 hZDCCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCCvsZEM->SetYTitle("E_{ZDCC} (TeV)");
141 fRecoParam->SetZDCCvsZEM(hZDCCvsZEM);
142 //
143 TH2F* hZDCAvsZEM = new TH2F("hZDCAvsZEM","hZDCAvsZEM",100,0.,10.,100,0.,500.);
144 hZDCAvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCAvsZEM->SetYTitle("E_{ZDCA} (TeV)");
145 fRecoParam->SetZDCAvsZEM(hZDCAvsZEM);
146 //
147 AliInfo("\n ***** CALIBRATION_MB data -> building AliZDCRecoParamPbPb object *****");
148 }
149 }
150
151 fBeamEnergy = grpData->GetBeamEnergy();
152 if(fBeamEnergy==AliGRPObject::GetInvalidFloat()){
153 AliWarning("GRP/GRP/Data entry: missing value for the beam energy ! Using 0.");
154 fBeamEnergy = 0.;
155 }
156
157 if(fIsCalibrationMB==kFALSE)
158 AliInfo(Form("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f GeV *****\n",beamType.Data(), fBeamEnergy));
159 }
160 else{
161 AliError(" ATTENTION!!!!!! No beam type nor beam energy has been set!!!!!!\n");
162 }
163
164}
165
166//_____________________________________________________________________________
167void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
168{
169 // *** Local ZDC reconstruction for digits
170 // Works on the current event
171
172 // Retrieving calibration data
173 // Parameters for mean value pedestal subtraction
174 int const kNch = 24;
175 Float_t meanPed[2*kNch];
176 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
177 // Parameters pedestal subtraction through correlation with out-of-time signals
178 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
179 for(Int_t jj=0; jj<2*kNch; jj++){
180 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
181 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
182 }
183
184 // get digits
185 AliZDCDigit digit;
186 AliZDCDigit* pdigit = &digit;
187 digitsTree->SetBranchAddress("ZDC", &pdigit);
188 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
189
190 // loop over digits
191 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
192 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
193 for(Int_t i=0; i<10; i++){
194 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
195 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
196 }
197
198 Int_t digNentries = digitsTree->GetEntries();
199 Float_t ootDigi[kNch];
200 // -- Reading out-of-time signals (last kNch entries) for current event
201 if(fPedSubMode==1){
202 for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
203 ootDigi[iDigit] = digitsTree->GetEntry(iDigit);
204 }
205 }
206
207 for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
208 digitsTree->GetEntry(iDigit);
209 if (!pdigit) continue;
210 //
211 Int_t det = digit.GetSector(0);
212 Int_t quad = digit.GetSector(1);
213 Int_t pedindex = -1;
214 Float_t ped2SubHg=0., ped2SubLg=0.;
215 if(quad!=5){
216 if(det==1) pedindex = quad;
217 else if(det==2) pedindex = quad+5;
218 else if(det==3) pedindex = quad+9;
219 else if(det==4) pedindex = quad+12;
220 else if(det==5) pedindex = quad+17;
221 }
222 else pedindex = (det-1)/3+22;
223 //
224 if(fPedSubMode==0){
225 ped2SubHg = meanPed[pedindex];
226 ped2SubLg = meanPed[pedindex+kNch];
227 }
228 else if(fPedSubMode==1){
229 ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
230 ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
231 }
232
233 if(quad != 5){ // ZDC (not reference PTMs!)
234 if(det == 1){ // *** ZNC
235 tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
236 tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
237 if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
238 if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
239 // Ch. debug
240 //printf("\t pedindex %d tZN1Corr[%d] = %1.0f tZN1Corr[%d] = %1.0f",
241 // pedindex, quad, tZN1Corr[quad], quad+5, tZN1Corr[quad+5]);
242 }
243 else if(det == 2){ // *** ZP1
244 tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
245 tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
246 if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
247 if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
248 // Ch. debug
249 //printf("\t pedindex %d tZP1Corr[%d] = %1.0f tZP1Corr[%d] = %1.0f",
250 // pedindex, quad, tZP1Corr[quad], quad+5, tZP1Corr[quad+5]);
251 }
252 else if(det == 3){
253 if(quad == 1){ // *** ZEM1
254 dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
255 dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
256 if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
257 if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
258 // Ch. debug
259 //printf("\t pedindex %d tZEM1Corr[%d] = %1.0f tZEM1Corr[%d] = %1.0f",
260 // pedindex, quad, tZEM1Corr[quad], quad+1, tZEM1Corr[quad+1]);
261 }
262 else if(quad == 2){ // *** ZEM2
263 dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
264 dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
265 if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
266 if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
267 // Ch. debug
268 //printf("\t pedindex %d tZEM2Corr[%d] = %1.0f tZEM2Corr[%d] = %1.0f",
269 // pedindex, quad, tZEM2Corr[quad], quad+1, tZEM2Corr[quad+1]);
270 }
271 }
272 else if(det == 4){ // *** ZN2
273 tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
274 tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
275 if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
276 if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
277 // Ch. debug
278 //printf("\t pedindex %d tZN2Corr[%d] = %1.0f tZN2Corr[%d] = %1.0f\n",
279 // pedindex, quad, tZN2Corr[quad], quad+5, tZN2Corr[quad+5]);
280 }
281 else if(det == 5){ // *** ZP2
282 tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
283 tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
284 if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
285 if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
286 // Ch. debug
287 //printf("\t pedindex %d tZP2Corr[%d] = %1.0f tZP2Corr[%d] = %1.0f\n",
288 // pedindex, quad, tZP2Corr[quad], quad+5, tZP2Corr[quad+5]);
289 }
290 }
291 else{ // Reference PMs
292 if(det == 1){
293 sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
294 sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
295 // Ch. debug
296 if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
297 if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
298 }
299 else if(det == 4){
300 sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
301 sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
302 // Ch. debug
303 if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
304 if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
305 }
306 }
307
308 // Ch. debug
309 /*printf(" - AliZDCReconstructor -> digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
310 iDigit, det, quad, ped2SubHg, ped2SubLg);
311 printf(" HGChain -> RawDig %d DigCorr %1.2f\n", digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg);
312 printf(" LGChain -> RawDig %d DigCorr %1.2f\n", digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);
313 */
314 }//digits loop
315
316
317 // If CALIBRATION_MB run build the RecoParam object
318 if(fIsCalibrationMB){
319 Float_t ZDCC=0., ZDCA=0., ZEM=0;
320 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
321 for(Int_t jkl=0; jkl<5; jkl++){
322 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
323 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
324 }
325 // Using energies in TeV in fRecoParam object
326 BuildRecoParam(ZDCC/1000., ZDCA/1000., ZEM/1000.);
327 }
328 Bool_t recFlag1 = kFALSE, recFlag2 = kFALSE, recFlag3 = kFALSE;
329 // reconstruct the event
330 if(fRecoMode==1)
331 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
332 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, recFlag1, recFlag2, recFlag3);
333 else if(fRecoMode==2)
334 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
335 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, recFlag1, recFlag2, recFlag3);
336}
337
338//_____________________________________________________________________________
339void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
340{
341 // *** ZDC raw data reconstruction
342 // Works on the current event
343
344 // Retrieving calibration data
345 // Parameters for pedestal subtraction
346 int const kNch = 24;
347 Float_t meanPed[2*kNch];
348 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
349 // Parameters pedestal subtraction through correlation with out-of-time signals
350 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
351 for(Int_t jj=0; jj<2*kNch; jj++){
352 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
353 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
354 }
355
356 Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
357 Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
358 Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
359 Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
360 Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
361 Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
362 for(Int_t ich=0; ich<5; ich++){
363 adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
364 adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
365 adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
366 adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
367 if(ich<2){
368 adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
369 pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
370 }
371 }
372
373 // loop over raw data
374 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
375 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
376 for(Int_t i=0; i<10; i++){
377 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
378 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
379 }
380
381 //fNRun = (Int_t) rawReader->GetRunNumber();
382 Bool_t chOff=kFALSE, isUndflw=kFALSE, isOvflw=kFALSE;
383 //
384 rawReader->Reset();
385 AliZDCRawStream rawData(rawReader);
386 while(rawData.Next()){
387 if(rawData.IsCalibration() == kFALSE){ // Reading ADCs
388 //printf(" **** Reading ADC raw data **** \n");
389 //
390 if(rawData.GetNChannelsOn() < 48 ) chOff=kTRUE;
391 if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)) isUndflw=kTRUE;
392 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)) isOvflw=kTRUE;
393
394 if((rawData.IsADCDataWord()) && (isUndflw==kFALSE) && (isOvflw==kFALSE)){
395
396 Int_t adcMod = rawData.GetADCModule();
397 Int_t det = rawData.GetSector(0);
398 Int_t quad = rawData.GetSector(1);
399 Int_t gain = rawData.GetADCGain();
400 Int_t pedindex=0;
401 //
402 // Mean pedestal value subtraction -------------------------------------------------------
403 if(fPedSubMode == 0){
404 // Not interested in o.o.t. signals (ADC modules 2, 3)
405 if(adcMod == 2 || adcMod == 3) continue;
406 //
407 if(quad != 5){ // ZDCs (not reference PTMs)
408 if(det == 1){
409 pedindex = quad;
410 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
411 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
412 }
413 else if(det == 2){
414 pedindex = quad+5;
415 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
416 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
417 }
418 else if(det == 3){
419 pedindex = quad+9;
420 if(quad==1){
421 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
422 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
423 }
424 else if(quad==2){
425 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
426 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
427 }
428 }
429 else if(det == 4){
430 pedindex = quad+12;
431 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
432 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
433 }
434 else if(det == 5){
435 pedindex = quad+17;
436 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
437 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
438 }
439 }
440 else{ // reference PM
441 pedindex = (det-1)/3 + 22;
442 if(det == 1){
443 if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
444 else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
445 }
446 else if(det == 4){
447 if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
448 else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
449 }
450 }
451 // Ch. debug
452 /*printf(" -> AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f\n",
453 det,quad,gain, pedindex, meanPed[pedindex]);
454 printf(" -> AliZDCReconstructor: RawADC %1.0f ADCCorr %1.0f\n",
455 rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);*/
456
457 }// mean pedestal subtraction
458 // Pedestal subtraction from correlation ------------------------------------------------
459 else if(fPedSubMode == 1){
460 // In time signals
461 if(adcMod==0 || adcMod==1){
462 if(quad != 5){ // signals from ZDCs
463 if(det == 1){
464 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
465 else adcZN1lg[quad] = rawData.GetADCValue();
466 }
467 else if(det == 2){
468 if(gain==0) adcZP1[quad] = rawData.GetADCValue();
469 else adcZP1lg[quad] = rawData.GetADCValue();
470 }
471 else if(det == 3){
472 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
473 else adcZEMlg[quad-1] = rawData.GetADCValue();
474 }
475 else if(det == 4){
476 if(gain==0) adcZN2[quad] = rawData.GetADCValue();
477 else adcZN2lg[quad] = rawData.GetADCValue();
478 }
479 else if(det == 5){
480 if(gain==0) adcZP2[quad] = rawData.GetADCValue();
481 else adcZP2lg[quad] = rawData.GetADCValue();
482 }
483 }
484 else{ // signals from reference PM
485 if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
486 else pmReflg[quad-1] = rawData.GetADCValue();
487 }
488 }
489 // Out-of-time pedestals
490 else if(adcMod==2 || adcMod==3){
491 if(quad != 5){ // signals from ZDCs
492 if(det == 1){
493 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
494 else adcZN1ootlg[quad] = rawData.GetADCValue();
495 }
496 else if(det == 2){
497 if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
498 else adcZP1ootlg[quad] = rawData.GetADCValue();
499 }
500 else if(det == 3){
501 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
502 else adcZEMootlg[quad-1] = rawData.GetADCValue();
503 }
504 else if(det == 4){
505 if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
506 else adcZN2ootlg[quad] = rawData.GetADCValue();
507 }
508 else if(det == 5){
509 if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
510 else adcZP2ootlg[quad] = rawData.GetADCValue();
511 }
512 }
513 else{ // signals from reference PM
514 if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
515 else pmRefootlg[quad-1] = rawData.GetADCValue();
516 }
517 }
518 } // pedestal subtraction from correlation
519 // Ch. debug
520 //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
521 // det,quad,gain, pedindex, meanPed[pedindex]);
522 }//IsADCDataWord
523 }// Not raw data from calibration run!
524 else{
525 }
526 }//loop on raw data
527
528 if(fPedSubMode==1){
529 for(Int_t t=0; t<5; t++){
530 tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
531 tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
532 //
533 tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
534 tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
535 //
536 tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
537 tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
538 //
539 tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
540 tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
541 // 0---------0 Ch. debug 0---------0
542/* printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
543 printf("\tCorrCoeff0\tCorrCoeff1\n");
544 printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
545 printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
546 printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
547 printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
548 printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
549 printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
550 printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
551 printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
552
553 printf("ZN1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
554 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
555 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
556 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
557 //
558 printf("ZP1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
559 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
560 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
561 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
562 //
563 printf("ZN2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
564 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
565 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
566 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
567 //
568 printf("ZP2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
569 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
570 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
571 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
572*/
573 }
574 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
575 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
576 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
577 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
578 //
579 sPMRef1[0] = pmRef[0] - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
580 sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
581 sPMRef2[0] = pmRef[0] - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
582 sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
583 }
584
585 // If CALIBRATION_MB run build the RecoParam object
586 if(fIsCalibrationMB){
587 Float_t ZDCC=0., ZDCA=0., ZEM=0;
588 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
589 for(Int_t jkl=0; jkl<5; jkl++){
590 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
591 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
592 }
593 BuildRecoParam(ZDCC/100., ZDCA/100., ZEM/100.);
594 }
595 // reconstruct the event
596 else{
597 if(fRecoMode==1) // p-p data
598 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
599 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, chOff, isUndflw, isOvflw);
600 else if(fRecoMode==2) // Pb-Pb data
601 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
602 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, chOff, isUndflw, isOvflw);
603 }
604}
605
606//_____________________________________________________________________________
607void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* corrADCZN1,
608 Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
609 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
610 Bool_t channelsOff, Bool_t chUnderflow, Bool_t chOverflow) const
611{
612 // ****************** Reconstruct one event ******************
613
614 // ---- Setting reco flags
615 UInt_t recoFlag=0;
616 UInt_t rFlags[32];
617 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
618
619 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
620 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
621 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
622 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
623 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
624 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
625 //
626 if(channelsOff==kTRUE) rFlags[8] = 0x1;
627 if(chUnderflow == kTRUE) rFlags[9] = 0x1;
628 if(chOverflow==kTRUE) rFlags[10] = 0x1;
629 recoFlag = rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
630 rFlags[5] << 5 | rFlags[4] << 4 | rFlags[3] << 3 |
631 rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
632
633 // ****** Retrieving calibration data
634 // --- Equalization coefficients ---------------------------------------------
635 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
636 for(Int_t ji=0; ji<5; ji++){
637 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
638 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
639 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
640 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
641 }
642 // --- Energy calibration factors ------------------------------------
643 Float_t calibEne[4];
644 // **** Energy calibration coefficient set to 1
645 // **** (no trivial way to calibrate in p-p runs)
646 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
647
648 // ****** Equalization of detector responses
649 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
650 for(Int_t gi=0; gi<10; gi++){
651 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
652 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
653 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
654 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
655 }
656
657 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
658 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
659 for(Int_t gi=0; gi<5; gi++){
660 calibSumZN1[0] += equalTowZN1[gi];
661 calibSumZP1[0] += equalTowZP1[gi];
662 calibSumZN2[0] += equalTowZN2[gi];
663 calibSumZP2[0] += equalTowZP2[gi];
664 //
665 calibSumZN1[1] += equalTowZN1[gi+5];
666 calibSumZP1[1] += equalTowZP1[gi+5];
667 calibSumZN2[1] += equalTowZN2[gi+5];
668 calibSumZP2[1] += equalTowZP2[gi+5];
669 }
670 // High gain chain
671 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
672 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
673 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
674 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
675 // Low gain chain
676 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
677 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
678 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
679 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
680
681 // ****** Energy calibration of detector responses
682 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
683 for(Int_t gi=0; gi<5; gi++){
684 // High gain chain
685 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
686 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
687 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
688 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
689 // Low gain chain
690 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
691 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
692 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
693 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
694 }
695 //
696 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
697 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
698 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
699 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
700 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
701 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
702
703 // ****** No. of spectator and participants nucleons
704 // Variables calculated to comply with ESD structure
705 // *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
706 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
707 Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
708 Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
709 Double_t impPar=0., impPar1=0., impPar2=0.;
710
711 // create the output tree
712 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
713 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
714 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
715 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
716 nGenSpec, nGenSpecLeft, nGenSpecRight,
717 nPart, nPartTotLeft, nPartTotRight,
718 impPar, impPar1, impPar2,
719 recoFlag);
720
721 //AliZDCReco* preco = &reco;
722 const Int_t kBufferSize = 8000;
723 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
724
725 // write the output tree
726 clustersTree->Fill();
727 delete reco;
728}
729
730//_____________________________________________________________________________
731void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
732 Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
733 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
734 Bool_t channelsOff, Bool_t chUnderflow, Bool_t chOverflow) const
735{
736 // ****************** Reconstruct one event ******************
737
738 // ---- Setting reco flags
739 UInt_t recoFlag=0;
740 UInt_t rFlags[32];
741 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
742
743 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
744 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
745 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
746 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
747 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
748 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
749 //
750 if(channelsOff==kTRUE) rFlags[8] = 0x1;
751 if(chUnderflow == kTRUE) rFlags[9] = 0x1;
752 if(chOverflow==kTRUE) rFlags[10] = 0x1;
753 recoFlag = rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
754 rFlags[5] << 5 | rFlags[4] << 4 | rFlags[3] << 3 |
755 rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
756
757 // ****** Retrieving calibration data
758 // --- Equalization coefficients ---------------------------------------------
759 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
760 for(Int_t ji=0; ji<5; ji++){
761 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
762 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
763 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
764 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
765 }
766 // --- Energy calibration factors ------------------------------------
767 Float_t valFromOCDB[6], calibEne[6];
768 for(Int_t ij=0; ij<6; ij++){
769 valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
770 if(ij<4){
771 if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
772 else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
773 }
774 else calibEne[ij] = valFromOCDB[ij];
775 }
776
777 // ****** Equalization of detector responses
778 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
779 for(Int_t gi=0; gi<10; gi++){
780 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
781 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
782 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
783 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
784 }
785
786 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
787 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
788 for(Int_t gi=0; gi<5; gi++){
789 calibSumZN1[0] += equalTowZN1[gi];
790 calibSumZP1[0] += equalTowZP1[gi];
791 calibSumZN2[0] += equalTowZN2[gi];
792 calibSumZP2[0] += equalTowZP2[gi];
793 //
794 calibSumZN1[1] += equalTowZN1[gi+5];
795 calibSumZP1[1] += equalTowZP1[gi+5];
796 calibSumZN2[1] += equalTowZN2[gi+5];
797 calibSumZP2[1] += equalTowZP2[gi+5];
798 }
799 // High gain chain
800 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
801 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
802 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
803 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
804 // Low gain chain
805 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
806 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
807 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
808 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
809 //
810 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
811 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
812 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
813 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
814 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
815 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
816
817 // ****** Energy calibration of detector responses
818 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
819 for(Int_t gi=0; gi<5; gi++){
820 // High gain chain
821 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
822 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
823 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
824 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
825 // Low gain chain
826 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
827 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
828 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
829 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
830 }
831
832 // ****** Number of detected spectator nucleons
833 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
834 if(fBeamEnergy!=0){
835 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
836 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
837 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
838 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
839 }
840 else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
841 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
842 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
843 nDetSpecNRight, nDetSpecPRight);*/
844
845 Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
846 Int_t nPart=0, nPartA=0, nPartC=0;
847 Double_t b=0., bA=0., bC=0.;
848
849 if(fIsCalibrationMB == kFALSE){
850 // ****** Reconstruction parameters ------------------
851 // Ch. debug
852 //fRecoParam->Print("");
853 //
854 TH2F *hZDCvsZEM = fRecoParam->GethZDCvsZEM();
855 TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
856 TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
857 TH1D *hNpartDist = fRecoParam->GethNpartDist();
858 TH1D *hbDist = fRecoParam->GethbDist();
859 Float_t ClkCenter = fRecoParam->GetClkCenter();
860 //
861 Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
862 Double_t origin = xHighEdge*ClkCenter;
863 // Ch. debug
864 //printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
865 //
866 // ====> Summed ZDC info (sideA+side C)
867 TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
868 Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
869 Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
870 line->SetParameter(0, y/(x-origin));
871 line->SetParameter(1, -origin*y/(x-origin));
872 // Ch. debug
873 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
874 //printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
875 //
876 Double_t countPerc=0;
877 Double_t xBinCenter=0, yBinCenter=0;
878 for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
879 for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
880 xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
881 yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
882 //
883 if(line->GetParameter(0)>0){
884 if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
885 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
886 // Ch. debug
887 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
888 xBinCenter, yBinCenter, countPerc);*/
889 }
890 }
891 else{
892 if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
893 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
894 // Ch. debug
895 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
896 xBinCenter, yBinCenter, countPerc);*/
897 }
898 }
899 }
900 }
901 //
902 Double_t xSecPerc = 0.;
903 if(hZDCvsZEM->GetEntries()!=0){
904 xSecPerc = countPerc/hZDCvsZEM->GetEntries();
905 }
906 else{
907 AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
908 }
909 // Ch. debug
910 //printf(" xSecPerc %1.4f \n", xSecPerc);
911
912 // ====> side C
913 TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
914 Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
915 lineC->SetParameter(0, yC/(x-origin));
916 lineC->SetParameter(1, -origin*yC/(x-origin));
917 // Ch. debug
918 //printf(" ***************** Side C \n");
919 //printf(" E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
920 //
921 Double_t countPercC=0;
922 Double_t xBinCenterC=0, yBinCenterC=0;
923 for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
924 for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
925 xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
926 yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
927 if(lineC->GetParameter(0)>0){
928 if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
929 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
930 }
931 }
932 else{
933 if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
934 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
935 }
936 }
937 }
938 }
939 //
940 Double_t xSecPercC = 0.;
941 if(hZDCCvsZEM->GetEntries()!=0){
942 xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
943 }
944 else{
945 AliWarning(" Histogram hZDCCvsZEM from OCDB has no entries!!!");
946 }
947 // Ch. debug
948 //printf(" xSecPercC %1.4f \n", xSecPercC);
949
950 // ====> side A
951 TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
952 Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
953 lineA->SetParameter(0, yA/(x-origin));
954 lineA->SetParameter(1, -origin*yA/(x-origin));
955 //
956 // Ch. debug
957 //printf(" ***************** Side A \n");
958 //printf(" E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
959 //
960 Double_t countPercA=0;
961 Double_t xBinCenterA=0, yBinCenterA=0;
962 for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
963 for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
964 xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
965 yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
966 if(lineA->GetParameter(0)>0){
967 if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
968 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
969 }
970 }
971 else{
972 if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
973 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
974 }
975 }
976 }
977 }
978 //
979 Double_t xSecPercA = 0.;
980 if(hZDCAvsZEM->GetEntries()!=0){
981 xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
982 }
983 else{
984 AliWarning(" Histogram hZDCAvsZEM from OCDB has no entries!!!");
985 }
986 // Ch. debug
987 //printf(" xSecPercA %1.4f \n", xSecPercA);
988
989 // ****** Number of participants (from E_ZDC vs. E_ZEM correlation)
990 Int_t nPart=0, nPartC=0, nPartA=0;
991 Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
992 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
993 nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
994 if((1.-nPartFrac) < xSecPerc){
995 nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
996 // Ch. debug
997 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
998 //printf(" nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
999 break;
1000 }
1001 }
1002 if(nPart<0) nPart=0;
1003 //
1004 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1005 nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1006 if((1.-nPartFracC) < xSecPercC){
1007 nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1008 // Ch. debug
1009 //printf(" ***************** Side C \n");
1010 //printf(" nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1011 break;
1012 }
1013 }
1014 if(nPartC<0) nPartC=0;
1015 //
1016 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1017 nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1018 if((1.-nPartFracA) < xSecPercA){
1019 nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1020 // Ch. debug
1021 //printf(" ***************** Side A \n");
1022 //printf(" nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1023 break;
1024 }
1025 }
1026 if(nPartA<0) nPartA=0;
1027
1028 // ****** Impact parameter (from E_ZDC vs. E_ZEM correlation)
1029 Float_t b=0, bC=0, bA=0;
1030 Double_t bFrac=0., bFracC=0., bFracA=0.;
1031 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1032 bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1033 if((1.-bFrac) < xSecPerc){
1034 b = hbDist->GetBinLowEdge(ibbin);
1035 break;
1036 }
1037 }
1038 //
1039 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1040 bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1041 if((1.-bFracC) < xSecPercC){
1042 bC = hbDist->GetBinLowEdge(ibbin);
1043 break;
1044 }
1045 }
1046 //
1047 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1048 bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1049 if((1.-bFracA) < xSecPercA){
1050 bA = hbDist->GetBinLowEdge(ibbin);
1051 break;
1052 }
1053 }
1054
1055 // ****** Number of spectator nucleons
1056 Int_t nGenSpec=0, nGenSpecC=0, nGenSpecA=0;
1057 //
1058 nGenSpec = 416 - nPart;
1059 nGenSpecC = 416 - nPartC;
1060 nGenSpecA = 416 - nPartA;
1061 if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1062 if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1063 if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1064
1065 // Ch. debug
1066 /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1067 " calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
1068 calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1069 printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n",
1070 nGenSpecLeft, nGenSpecRight);
1071 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1072 */
1073
1074 delete lineC; delete lineA;
1075
1076 } // ONLY IF fIsCalibrationMB==kFALSE
1077
1078 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1079 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1080 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1081 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1082 nGenSpec, nGenSpecA, nGenSpecC,
1083 nPart, nPartA, nPartC, b, bA, bC,
1084 recoFlag);
1085
1086 const Int_t kBufferSize = 4000;
1087 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1088 // write the output tree
1089 clustersTree->Fill();
1090 delete reco;
1091}
1092
1093//_____________________________________________________________________________
1094void AliZDCReconstructor::BuildRecoParam(Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1095{
1096 // Calculate RecoParam object for Pb-Pb data
1097 (fRecoParam->GethZDCvsZEM())->Fill(ZDCC+ZDCA, ZEM);
1098 (fRecoParam->GethZDCCvsZEM())->Fill(ZDCC, ZEM);
1099 (fRecoParam->GethZDCAvsZEM())->Fill(ZDCA, ZEM);
1100
1101}
1102
1103//_____________________________________________________________________________
1104void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1105{
1106 // fill energies and number of participants to the ESD
1107
1108 if(fIsCalibrationMB==kTRUE) WritePbPbRecoParamInOCDB();
1109
1110 AliZDCReco reco;
1111 AliZDCReco* preco = &reco;
1112 clustersTree->SetBranchAddress("ZDC", &preco);
1113
1114 clustersTree->GetEntry(0);
1115 //
1116 AliESDZDC * esdzdc = esd->GetESDZDC();
1117 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1118 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1119 for(Int_t i=0; i<5; i++){
1120 tZN1Ene[i] = reco.GetZN1HREnTow(i);
1121 tZN2Ene[i] = reco.GetZN2HREnTow(i);
1122 tZP1Ene[i] = reco.GetZP1HREnTow(i);
1123 tZP2Ene[i] = reco.GetZP2HREnTow(i);
1124 //
1125 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1126 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1127 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1128 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1129 }
1130 //
1131 esdzdc->SetZN1TowerEnergy(tZN1Ene);
1132 esdzdc->SetZN2TowerEnergy(tZN2Ene);
1133 esdzdc->SetZP1TowerEnergy(tZP1Ene);
1134 esdzdc->SetZP2TowerEnergy(tZP2Ene);
1135 //
1136 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1137 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1138 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1139 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1140 //
1141 Int_t nPart = reco.GetNParticipants();
1142 Int_t nPartA = reco.GetNPartSideA();
1143 Int_t nPartC = reco.GetNPartSideC();
1144 Double_t b = reco.GetImpParameter();
1145 Double_t bA = reco.GetImpParSideA();
1146 Double_t bC = reco.GetImpParSideC();
1147 UInt_t recoFlag = reco.GetRecoFlag();
1148 //
1149 esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
1150 reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1151 nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1152
1153}
1154
1155//_____________________________________________________________________________
1156AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1157{
1158 // Setting the storage
1159
1160 Bool_t deleteManager = kFALSE;
1161
1162 AliCDBManager *manager = AliCDBManager::Instance();
1163 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1164
1165 if(!defstorage || !(defstorage->Contains("ZDC"))){
1166 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1167 manager->SetDefaultStorage(uri);
1168 deleteManager = kTRUE;
1169 }
1170
1171 AliCDBStorage *storage = manager->GetDefaultStorage();
1172
1173 if(deleteManager){
1174 AliCDBManager::Instance()->UnsetDefaultStorage();
1175 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1176 }
1177
1178 return storage;
1179}
1180
1181//_____________________________________________________________________________
1182AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1183{
1184
1185 // Getting pedestal calibration object for ZDC set
1186
1187 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1188 if(!entry) AliFatal("No calibration data loaded!");
1189
1190 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
1191 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1192
1193 return calibdata;
1194}
1195
1196//_____________________________________________________________________________
1197AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1198{
1199
1200 // Getting energy and equalization calibration object for ZDC set
1201
1202 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1203 if(!entry) AliFatal("No calibration data loaded!");
1204
1205 AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1206 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1207
1208 return calibdata;
1209}
1210
1211//_____________________________________________________________________________
1212AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1213{
1214
1215 // Getting energy and equalization calibration object for ZDC set
1216
1217 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1218 if(!entry) AliFatal("No calibration data loaded!");
1219
1220 AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1221 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1222
1223 return calibdata;
1224}
1225
1226//_____________________________________________________________________________
1227AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1228{
1229
1230 // Getting reconstruction parameters from OCDB
1231
1232 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1233 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1234
1235 AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1236 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1237
1238 return param;
1239
1240}
1241
1242//_____________________________________________________________________________
1243AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1244{
1245
1246 // Getting reconstruction parameters from OCDB
1247
1248 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1249 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1250
1251 AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1252 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1253
1254 return param;
1255
1256}
1257
1258//_____________________________________________________________________________
1259void AliZDCReconstructor::WritePbPbRecoParamInOCDB() const
1260{
1261
1262 // Writing Pb-Pb reconstruction parameters from OCDB
1263
1264 AliCDBManager *man = AliCDBManager::Instance();
1265 AliCDBMetaData *md= new AliCDBMetaData();
1266 md->SetResponsible("Chiara Oppedisano");
1267 md->SetComment("ZDC Pb-Pb reconstruction parameters");
1268 md->SetObjectClassName("AliZDCRecoParamPbPb");
1269 AliCDBId id("ZDC/Calib/RecoParamPbPb",fNRun,AliCDBRunRange::Infinity());
1270 man->Put(fRecoParam, id, md);
1271
1272}
1273