]> 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
Adding a reminder for coders
[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 "AliESDEvent.h"
37#include "AliESDZDC.h"
38#include "AliZDCDigit.h"
39#include "AliZDCRawStream.h"
40#include "AliZDCReco.h"
41#include "AliZDCReconstructor.h"
42#include "AliZDCPedestals.h"
43#include "AliZDCEnCalib.h"
44#include "AliZDCTowerCalib.h"
45#include "AliZDCMBCalib.h"
46#include "AliZDCRecoParam.h"
47#include "AliZDCRecoParampp.h"
48#include "AliZDCRecoParamPbPb.h"
49#include "AliRunInfo.h"
50#include "AliLHCClockPhase.h"
51
52
53ClassImp(AliZDCReconstructor)
54AliZDCRecoParam *AliZDCReconstructor::fgRecoParam=0; //reconstruction parameters
55AliZDCMBCalib *AliZDCReconstructor::fgMBCalibData=0; //calibration parameters for A-A reconstruction
56
57//_____________________________________________________________________________
58AliZDCReconstructor:: AliZDCReconstructor() :
59 fPedData(GetPedestalData()),
60 fEnCalibData(GetEnergyCalibData()),
61 fTowCalibData(GetTowerCalibData()),
62 fRecoMode(0),
63 fBeamEnergy(0.),
64 fNRun(0),
65 fIsCalibrationMB(kFALSE),
66 fPedSubMode(0),
67 fSignalThreshold(7),
68 fMeanPhase(0),
69 fESDZDC(NULL)
70{
71 // **** Default constructor
72}
73
74
75//_____________________________________________________________________________
76AliZDCReconstructor::~AliZDCReconstructor()
77{
78// destructor
79// if(fgRecoParam) delete fgRecoParam;
80 if(fPedData) delete fPedData;
81 if(fEnCalibData) delete fEnCalibData;
82 if(fTowCalibData) delete fTowCalibData;
83 if(fgMBCalibData) delete fgMBCalibData;
84 if(fESDZDC) delete fESDZDC;
85}
86
87//____________________________________________________________________________
88void AliZDCReconstructor::Init()
89{
90 // Setting reconstruction parameters
91
92 TString runType = GetRunInfo()->GetRunType();
93 if((runType.CompareTo("CALIBRATION_MB")) == 0){
94 fIsCalibrationMB = kTRUE;
95 }
96
97 TString beamType = GetRunInfo()->GetBeamType();
98 // This is a temporary solution to allow reconstruction in tests without beam
99 if(((beamType.CompareTo("UNKNOWN"))==0) &&
100 ((runType.CompareTo("PHYSICS"))==0 || (runType.CompareTo("CALIBRATION_BC"))==0)){
101 fRecoMode=1;
102 }
103 /*else if((beamType.CompareTo("UNKNOWN"))==0){
104 AliError("\t UNKNOWN beam type\n");
105 return;
106 }*/
107
108 fBeamEnergy = GetRunInfo()->GetBeamEnergy();
109 if(fBeamEnergy<0.01){
110 AliWarning(" Beam energy value missing -> setting it to 1380 GeV ");
111 fBeamEnergy = 1380.;
112 }
113
114 if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
115 ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)){
116 fRecoMode=1;
117 }
118 else if((beamType.CompareTo("A-A")) == 0 || (beamType.CompareTo("AA")) == 0){
119 fRecoMode=2;
120 if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
121 if(fgRecoParam){
122 fgRecoParam->SetGlauberMCDist(fBeamEnergy);
123 }
124 }
125
126 AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase");
127 if (!entry) AliFatal("LHC clock-phase shift is not found in OCDB !");
128 AliLHCClockPhase *phaseLHC = (AliLHCClockPhase*)entry->GetObject();
129 // 4/2/2011 According to A. Di Mauro BEAM1 measurement is more reliable
130 // than BEAM2 and therefore also than the average of the 2
131 fMeanPhase = phaseLHC->GetMeanPhaseB1();
132
133 if(fIsCalibrationMB==kFALSE)
134 AliInfo(Form("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f + %1.0f GeV *****\n\n",
135 beamType.Data(), fBeamEnergy, fBeamEnergy));
136
137 // if EMD calibration run NO ENERGY CALIBRATION should be performed
138 // pp-like reconstruction must be performed (E cailb. coeff. = 1)
139 if((runType.CompareTo("CALIBRATION_EMD")) == 0){
140 fRecoMode=1;
141 fBeamEnergy = 1380.;
142 }
143
144 fESDZDC = new AliESDZDC();
145
146}
147
148
149//____________________________________________________________________________
150void AliZDCReconstructor::Init(TString beamType, Float_t beamEnergy)
151{
152 // Setting reconstruction mode
153 // Needed to work in the HLT framework
154
155 fIsCalibrationMB = kFALSE;
156
157 fBeamEnergy = beamEnergy;
158
159 if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
160 ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)){
161 fRecoMode=1;
162 }
163 else if((beamType.CompareTo("A-A")) == 0 || (beamType.CompareTo("AA")) == 0){
164 fRecoMode=2;
165 if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
166 if( fgRecoParam ) fgRecoParam->SetGlauberMCDist(fBeamEnergy);
167 }
168
169 AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase");
170 if (!entry) AliFatal("LHC clock-phase shift is not found in OCDB !");
171 AliLHCClockPhase *phaseLHC = (AliLHCClockPhase*)entry->GetObject();
172 fMeanPhase = phaseLHC->GetMeanPhase();
173
174 fESDZDC = new AliESDZDC();
175
176 AliInfo(Form("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f + %1.0f GeV *****\n\n",
177 beamType.Data(), fBeamEnergy, fBeamEnergy));
178
179}
180
181//_____________________________________________________________________________
182void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
183{
184 // *** Local ZDC reconstruction for digits
185 // Works on the current event
186
187 // Retrieving calibration data
188 // Parameters for mean value pedestal subtraction
189 int const kNch = 24;
190 Float_t meanPed[2*kNch];
191 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
192 // Parameters pedestal subtraction through correlation with out-of-time signals
193 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
194 for(Int_t jj=0; jj<2*kNch; jj++){
195 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
196 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
197 }
198
199 // get digits
200 AliZDCDigit digit;
201 AliZDCDigit* pdigit = &digit;
202 digitsTree->SetBranchAddress("ZDC", &pdigit);
203 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
204
205 // loop over digits
206 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
207 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
208 for(Int_t i=0; i<10; i++){
209 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
210 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
211 }
212
213 Int_t digNentries = digitsTree->GetEntries();
214 Float_t ootDigi[kNch]; Int_t i=0;
215 // -- Reading out-of-time signals (last kNch entries) for current event
216 if(fPedSubMode==1){
217 for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
218 if(i<=kNch) ootDigi[i] = digitsTree->GetEntry(iDigit);
219 else AliWarning(" Can't read more out of time values: index>kNch !!!\n");
220 i++;
221 }
222 }
223
224 for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
225 digitsTree->GetEntry(iDigit);
226 if (!pdigit) continue;
227 //
228 Int_t det = digit.GetSector(0);
229 Int_t quad = digit.GetSector(1);
230 Int_t pedindex = -1;
231 Float_t ped2SubHg=0., ped2SubLg=0.;
232 if(quad!=5){
233 if(det==1) pedindex = quad;
234 else if(det==2) pedindex = quad+5;
235 else if(det==3) pedindex = quad+9;
236 else if(det==4) pedindex = quad+12;
237 else if(det==5) pedindex = quad+17;
238 }
239 else pedindex = (det-1)/3+22;
240 //
241 if(fPedSubMode==0){
242 ped2SubHg = meanPed[pedindex];
243 ped2SubLg = meanPed[pedindex+kNch];
244 }
245 else if(fPedSubMode==1){
246 ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
247 ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
248 }
249
250 if(quad != 5){ // ZDC (not reference PTMs!)
251 if(det == 1){ // *** ZNC
252 tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
253 tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
254 }
255 else if(det == 2){ // *** ZP1
256 tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
257 tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
258 }
259 else if(det == 3){
260 if(quad == 1){ // *** ZEM1
261 dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
262 dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
263 }
264 else if(quad == 2){ // *** ZEM2
265 dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
266 dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
267 }
268 }
269 else if(det == 4){ // *** ZN2
270 tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
271 tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
272 }
273 else if(det == 5){ // *** ZP2
274 tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
275 tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
276 }
277 }
278 else{ // Reference PMs
279 if(det == 1){
280 sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
281 sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
282 }
283 else if(det == 4){
284 sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
285 sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
286 }
287 }
288
289 // Ch. debug
290 /*printf("AliZDCReconstructor: digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
291 iDigit, det, quad, ped2SubHg, ped2SubLg);
292 printf(" -> pedindex %d\n", pedindex);
293 printf(" HGChain -> RawDig %d DigCorr %1.2f",
294 digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg);
295 printf(" LGChain -> RawDig %d DigCorr %1.2f\n",
296 digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);*/
297
298 }//digits loop
299
300 UInt_t counts[32];
301 Int_t tdc[32][4];
302 for(Int_t jj=0; jj<32; jj++){
303 counts[jj]=0;
304 for(Int_t ii=0; ii<4; ii++) tdc[jj][ii]=0;
305 }
306
307 Int_t evQualityBlock[4] = {1,0,0,0};
308 Int_t triggerBlock[4] = {0,0,0,0};
309 Int_t chBlock[3] = {0,0,0};
310 UInt_t puBits=0;
311
312 // reconstruct the event
313 if(fRecoMode==1)
314 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
315 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
316 kFALSE, counts, tdc,
317 evQualityBlock, triggerBlock, chBlock, puBits);
318 else if(fRecoMode==2)
319 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
320 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
321 kFALSE, counts, tdc,
322 evQualityBlock, triggerBlock, chBlock, puBits);
323}
324
325//_____________________________________________________________________________
326void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
327{
328 // *** ZDC raw data reconstruction
329 // Works on the current event
330
331 // Retrieving calibration data
332 // Parameters for pedestal subtraction
333 int const kNch = 24;
334 Float_t meanPed[2*kNch];
335 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
336 // Parameters pedestal subtraction through correlation with out-of-time signals
337 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
338 for(Int_t jj=0; jj<2*kNch; jj++){
339 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
340 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
341 //printf(" %d %1.4f %1.4f\n", jj,corrCoeff0[jj],corrCoeff1[jj]);
342 }
343
344 Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
345 Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
346 Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
347 Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
348 Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
349 Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
350 for(Int_t ich=0; ich<5; ich++){
351 adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
352 adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
353 adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
354 adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
355 if(ich<2){
356 adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
357 pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
358 }
359 }
360
361 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
362 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
363 for(Int_t i=0; i<10; i++){
364 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
365 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
366 }
367
368 Bool_t isScalerOn=kFALSE;
369 Int_t jsc=0, itdc=0, iprevtdc=-1, ihittdc=0;
370 UInt_t scalerData[32];
371 Int_t tdcData[32][4];
372 for(Int_t k=0; k<32; k++){
373 scalerData[k]=0;
374 for(Int_t i=0; i<4; i++) tdcData[k][i]=0;
375 }
376
377
378 Int_t evQualityBlock[4] = {1,0,0,0};
379 Int_t triggerBlock[4] = {0,0,0,0};
380 Int_t chBlock[3] = {0,0,0};
381 UInt_t puBits=0;
382
383 Int_t kFirstADCGeo=0, kLastADCGeo=3, kScalerGeo=8, kZDCTDCGeo=4, kPUGeo=29;
384 //Int_t kTrigScales=30, kTrigHistory=31;
385
386 // loop over raw data
387 //rawReader->Reset();
388 AliZDCRawStream rawData(rawReader);
389 while(rawData.Next()){
390
391 // ***************************** Reading ADCs
392 if((rawData.GetADCModule()>=kFirstADCGeo) && (rawData.GetADCModule()<=kLastADCGeo)){
393 //printf(" **** Reading ADC raw data from module %d **** \n",rawData.GetADCModule());
394 //
395 if((rawData.IsADCDataWord()) && (rawData.GetNChannelsOn()<48)) chBlock[0] = kTRUE;
396 if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)) chBlock[1] = kTRUE;
397 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)) chBlock[2] = kTRUE;
398 if((rawData.IsADCDataWord()) && (rawData.IsADCEventGood() == kTRUE)) evQualityBlock[0] = kTRUE;
399
400 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow()==kFALSE)
401 && (rawData.IsOverflow()==kFALSE) && (rawData.IsADCEventGood()==kTRUE)){
402
403 Int_t adcMod = rawData.GetADCModule();
404 Int_t det = rawData.GetSector(0);
405 Int_t quad = rawData.GetSector(1);
406 Int_t gain = rawData.GetADCGain();
407 Int_t pedindex=0;
408 //
409 // Mean pedestal value subtraction -------------------------------------------------------
410 if(fPedSubMode == 0){
411 // **** Pb-Pb data taking 2010 -> subtracting some ch. from correlation ****
412 // Not interested in o.o.t. signals (ADC modules 2, 3)
413 //if(adcMod == 2 || adcMod == 3) continue;
414 if(((det==1 && quad==0) || (det==3))){
415 if(det == 1){
416 if(adcMod==0 || adcMod==1){
417 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
418 else adcZN1lg[quad] = rawData.GetADCValue();
419 }
420 else if(adcMod==2 || adcMod==3){
421 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
422 else adcZN1ootlg[quad] = rawData.GetADCValue();
423 }
424 }
425 else if(det == 3){
426 if(adcMod==0 || adcMod==1){
427 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
428 else adcZEMlg[quad-1] = rawData.GetADCValue();
429 }
430 else if(adcMod==2 || adcMod==3){
431 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
432 else adcZEMootlg[quad-1] = rawData.GetADCValue();
433 }
434 }
435 }
436 // When oot values are read the ADC modules 2, 3 can be skipped!!!
437 if(adcMod == 2 || adcMod == 3) continue;
438
439 // *************************************************************************
440 if(quad != 5){ // ZDCs (not reference PTMs)
441 if(det==1 && quad!=0){
442 pedindex = quad;
443 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
444 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
445 }
446 else if(det==2){
447 pedindex = quad+5;
448 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
449 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
450 }
451 /*else if(det == 3){
452 pedindex = quad+9;
453 if(quad==1){
454 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
455 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
456 }
457 else if(quad==2){
458 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
459 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
460 }
461 }*/
462 else if(det == 4){
463 pedindex = quad+12;
464 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
465 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
466 }
467 else if(det == 5){
468 pedindex = quad+17;
469 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
470 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
471 }
472 }
473 else{ // reference PM
474 pedindex = (det-1)/3 + 22;
475 if(det == 1){
476 if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
477 else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
478 }
479 else if(det == 4){
480 if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
481 else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
482 }
483 }
484 // Ch. debug
485 /*if(gain==0){
486 printf(" AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f",
487 det,quad,gain, pedindex, meanPed[pedindex]);
488 printf(" RawADC %d ADCCorr %1.0f\n",
489 rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);
490 }*/
491 }// mean pedestal subtraction
492 // Pedestal subtraction from correlation ------------------------------------------------
493 else if(fPedSubMode == 1){
494 // In time signals
495 if(adcMod==0 || adcMod==1){
496 if(quad != 5){ // signals from ZDCs
497 if(det == 1){
498 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
499 else adcZN1lg[quad] = rawData.GetADCValue();
500 }
501 else if(det == 2){
502 if(gain==0) adcZP1[quad] = rawData.GetADCValue();
503 else adcZP1lg[quad] = rawData.GetADCValue();
504 }
505 else if(det == 3){
506 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
507 else adcZEMlg[quad-1] = rawData.GetADCValue();
508 }
509 else if(det == 4){
510 if(gain==0) adcZN2[quad] = rawData.GetADCValue();
511 else adcZN2lg[quad] = rawData.GetADCValue();
512 }
513 else if(det == 5){
514 if(gain==0) adcZP2[quad] = rawData.GetADCValue();
515 else adcZP2lg[quad] = rawData.GetADCValue();
516 }
517 }
518 else{ // signals from reference PM
519 if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
520 else pmReflg[quad-1] = rawData.GetADCValue();
521 }
522 }
523 // Out-of-time pedestals
524 else if(adcMod==2 || adcMod==3){
525 if(quad != 5){ // signals from ZDCs
526 if(det == 1){
527 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
528 else adcZN1ootlg[quad] = rawData.GetADCValue();
529 }
530 else if(det == 2){
531 if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
532 else adcZP1ootlg[quad] = rawData.GetADCValue();
533 }
534 else if(det == 3){
535 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
536 else adcZEMootlg[quad-1] = rawData.GetADCValue();
537 }
538 else if(det == 4){
539 if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
540 else adcZN2ootlg[quad] = rawData.GetADCValue();
541 }
542 else if(det == 5){
543 if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
544 else adcZP2ootlg[quad] = rawData.GetADCValue();
545 }
546 }
547 else{ // signals from reference PM
548 if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
549 else pmRefootlg[quad-1] = rawData.GetADCValue();
550 }
551 }
552 } // pedestal subtraction from correlation
553 // Ch. debug
554 /*printf("\t AliZDCReconstructor: det %d quad %d res %d -> Ped[%d] = %1.0f\n",
555 det,quad,gain, pedindex, meanPed[pedindex]);*/
556 }//IsADCDataWord
557 }// ADC DATA
558 // ***************************** Reading Scaler
559 else if(rawData.GetADCModule()==kScalerGeo){
560 if(rawData.IsScalerWord()==kTRUE && rawData.IsScEventGood()==kTRUE){
561 isScalerOn = kTRUE;
562 scalerData[jsc] = rawData.GetTriggerCount();
563 // Ch. debug
564 //printf(" Reconstructed VME Scaler: %d %d ",jsc,scalerData[jsc]);
565 //
566 jsc++;
567 }
568 }// VME SCALER DATA
569 // ***************************** Reading ZDC TDC
570 else if(rawData.GetADCModule()==kZDCTDCGeo && rawData.IsZDCTDCDatum()==kTRUE){
571 itdc = rawData.GetChannel();
572 if(itdc==iprevtdc) ihittdc++;
573 else ihittdc=0;
574 iprevtdc=itdc;
575 tdcData[itdc][ihittdc] = rawData.GetZDCTDCDatum();
576 // Ch. debug
577 //printf(" Reconstructed TDC[%d, %d] %d ",itdc, ihittdc, tdcData[itdc][ihittdc]);
578 }// ZDC TDC DATA
579 // ***************************** Reading PU
580 else if(rawData.GetADCModule()==kPUGeo){
581 puBits = rawData.GetDetectorPattern();
582 }
583 // ***************************** Reading trigger history
584 else if(rawData.IstriggerHistoryWord()==kTRUE){
585 triggerBlock[0] = rawData.IsCPTInputEMDTrigger();
586 triggerBlock[1] = rawData.IsCPTInputSemiCentralTrigger();
587 triggerBlock[2] = rawData.IsCPTInputCentralTrigger();
588 triggerBlock[3] = rawData.IsCPTInputMBTrigger();
589 }
590
591 }//loop on raw data
592
593 if(fPedSubMode==1){
594 for(Int_t t=0; t<5; t++){
595 tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
596 tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
597 //
598 tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
599 tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
600 //
601 tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
602 tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
603 //
604 tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
605 tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
606 }
607 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[10]*adcZEMoot[0]+corrCoeff0[10]);
608 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[10+kNch]*adcZEMootlg[0]+corrCoeff0[10+kNch]);
609 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[11]*adcZEMoot[1]+corrCoeff0[11]);
610 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[11+kNch]*adcZEMootlg[1]+corrCoeff0[11+kNch]);
611 //
612 sPMRef1[0] = pmRef[0] - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
613 sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
614 sPMRef2[0] = pmRef[0] - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
615 sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
616 }
617 if(fPedSubMode==0 && fRecoMode==2){
618 // **** Pb-Pb data taking 2010 -> subtracting some ch. from correlation ****
619 tZN1Corr[0] = adcZN1[0] - (corrCoeff1[0]*adcZN1oot[0]+corrCoeff0[0]);
620 tZN1Corr[5] = adcZN1lg[0] - (corrCoeff1[kNch]*adcZN1ootlg[0]+corrCoeff0[kNch]);
621 // Ch. debug
622 //printf(" adcZN1 %d adcZN1oot %d tZN1Corr %1.2f \n", adcZN1[0],adcZN1oot[0],tZN1Corr[0]);
623 //printf(" adcZN1lg %d adcZN1ootlg %d tZN1Corrlg %1.2f \n", adcZN1lg[0],adcZN1ootlg[0],tZN1Corr[5]);
624 //
625 //tZP1Corr[2] = adcZP1[2] - (corrCoeff1[2+5]*adcZP1oot[2]+corrCoeff0[2+5]);
626 //tZP1Corr[2+5] = adcZP1lg[2] - (corrCoeff1[2+5+kNch]*adcZP1ootlg[2]+corrCoeff0[2+5+kNch]);
627 //
628 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[10]*adcZEMoot[0]+corrCoeff0[10]);
629 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[10+kNch]*adcZEMootlg[0]+corrCoeff0[10+kNch]);
630 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[11]*adcZEMoot[1]+corrCoeff0[11]);
631 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[11+kNch]*adcZEMootlg[1]+corrCoeff0[11+kNch]);
632 // *************************************************************************
633 }
634 else if(fPedSubMode==0 && fRecoMode==1){
635 // **** p-p data taking 2011 -> temporary patch to overcome DA problem ****
636 tZN1Corr[0] = adcZN1[0] - meanPed[0];
637 tZN1Corr[5] = adcZN1lg[0] - meanPed[kNch];
638 //
639 dZEM1Corr[0] = adcZEM[0] - meanPed[10];
640 dZEM1Corr[1] = adcZEMlg[0] - meanPed[10+kNch];
641 dZEM2Corr[0] = adcZEM[1] - meanPed[11];
642 dZEM2Corr[1] = adcZEMlg[1] - meanPed[11+kNch];
643 // *************************************************************************
644 }
645
646 if(fRecoMode==1) // p-p data
647 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
648 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
649 isScalerOn, scalerData, tdcData,
650 evQualityBlock, triggerBlock, chBlock, puBits);
651 else if(fRecoMode==2) // Pb-Pb data
652 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
653 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2,
654 isScalerOn, scalerData, tdcData,
655 evQualityBlock, triggerBlock, chBlock, puBits);
656}
657
658//_____________________________________________________________________________
659void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree,
660 const Float_t* const corrADCZN1, const Float_t* const corrADCZP1,
661 const Float_t* const corrADCZN2, const Float_t* const corrADCZP2,
662 const Float_t* const corrADCZEM1, const Float_t* const corrADCZEM2,
663 Float_t* sPMRef1, Float_t* sPMRef2, Bool_t isScalerOn, UInt_t* scaler,
664 Int_t tdcData[32][4], const Int_t* const evQualityBlock,
665 const Int_t* const triggerBlock, const Int_t* const chBlock, UInt_t puBits) const
666{
667 // ****************** Reconstruct one event ******************
668
669 // CH. debug
670 /*printf("\n*************************************************\n");
671 printf(" ReconstructEventpp -> values after pedestal subtraction:\n");
672 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
673 corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
674 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
675 corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
676 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
677 corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
678 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
679 corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
680 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
681 printf("*************************************************\n");*/
682
683 // ---------------------- Setting reco flags for ESD
684 UInt_t rFlags[32];
685 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
686
687 if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
688 else rFlags[31] = 0x1;
689 //
690 if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
691 if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
692 if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
693
694 if(triggerBlock[0] == 1) rFlags[27] = 0x1;
695 if(triggerBlock[1] == 1) rFlags[26] = 0x1;
696 if(triggerBlock[2] == 1) rFlags[25] = 0x1;
697 if(triggerBlock[3] == 1) rFlags[24] = 0x1;
698
699 if(chBlock[0] == 1) rFlags[18] = 0x1;
700 if(chBlock[1] == 1) rFlags[17] = 0x1;
701 if(chBlock[2] == 1) rFlags[16] = 0x1;
702
703
704 rFlags[13] = puBits & 0x00000020;
705 rFlags[12] = puBits & 0x00000010;
706 rFlags[11] = puBits & 0x00000080;
707 rFlags[10] = puBits & 0x00000040;
708 rFlags[9] = puBits & 0x00000020;
709 rFlags[8] = puBits & 0x00000010;
710
711 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
712 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
713 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
714 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
715 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
716 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
717
718 UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
719 rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
720 0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
721 0x0 << 19 | rFlags[18] << 18 | rFlags[17] << 17 | rFlags[16] << 16 |
722 0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
723 rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
724 0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
725 rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
726 // --------------------------------------------------
727
728 // ****** Retrieving calibration data
729 // --- Equalization coefficients ---------------------------------------------
730 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
731 for(Int_t ji=0; ji<5; ji++){
732 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
733 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
734 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
735 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
736 }
737 // --- Energy calibration factors ------------------------------------
738 Float_t calibEne[6];
739 // **** Energy calibration coefficient set to 1
740 // **** (no trivial way to calibrate in p-p runs)
741 for(Int_t ij=0; ij<6; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
742
743 // ****** Equalization of detector responses
744 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
745 for(Int_t gi=0; gi<10; gi++){
746 if(gi<5){
747 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
748 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
749 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
750 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
751 }
752 else{
753 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
754 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
755 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
756 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
757 }
758 }
759 // Ch. debug
760 /*printf("\n ------------- EQUALIZATION -------------\n");
761 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
762 equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
763 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
764 equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
765 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
766 equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
767 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
768 equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
769 printf(" ----------------------------------------\n");*/
770
771 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
772 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
773 for(Int_t gi=0; gi<5; gi++){
774 calibSumZN1[0] += equalTowZN1[gi];
775 calibSumZP1[0] += equalTowZP1[gi];
776 calibSumZN2[0] += equalTowZN2[gi];
777 calibSumZP2[0] += equalTowZP2[gi];
778 //
779 calibSumZN1[1] += equalTowZN1[gi+5];
780 calibSumZP1[1] += equalTowZP1[gi+5];
781 calibSumZN2[1] += equalTowZN2[gi+5];
782 calibSumZP2[1] += equalTowZP2[gi+5];
783 }
784 // High gain chain
785 calibSumZN1[0] = calibSumZN1[0]*calibEne[0];
786 calibSumZP1[0] = calibSumZP1[0]*calibEne[1];
787 calibSumZN2[0] = calibSumZN2[0]*calibEne[2];
788 calibSumZP2[0] = calibSumZP2[0]*calibEne[3];
789 // Low gain chain
790 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
791 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
792 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
793 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
794
795 // ****** Energy calibration of detector responses
796 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
797 for(Int_t gi=0; gi<5; gi++){
798 // High gain chain
799 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0];
800 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1];
801 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2];
802 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3];
803 // Low gain chain
804 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
805 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
806 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
807 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
808 }
809 //
810 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
811 calibZEM1[0] = corrADCZEM1[0]*calibEne[4];
812 calibZEM1[1] = corrADCZEM1[1]*calibEne[4];
813 calibZEM2[0] = corrADCZEM2[0]*calibEne[5];
814 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
815 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
816 // Ch. debug
817 /*printf("\n ------------- CALIBRATION -------------\n");
818 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
819 calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
820 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
821 calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
822 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
823 calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
824 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
825 calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
826 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
827 printf(" ----------------------------------------\n");*/
828
829 // ****** No. of spectator and participants nucleons
830 // Variables calculated to comply with ESD structure
831 // *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
832 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
833 Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
834 Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
835 Double_t impPar=0., impPar1=0., impPar2=0.;
836
837 Bool_t energyFlag = kFALSE;
838 // create the output tree
839 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
840 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
841 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
842 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
843 nGenSpec, nGenSpecLeft, nGenSpecRight,
844 nPart, nPartTotLeft, nPartTotRight,
845 impPar, impPar1, impPar2,
846 recoFlag, energyFlag, isScalerOn, scaler, tdcData);
847
848 const Int_t kBufferSize = 4000;
849 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
850 // write the output tree
851 clustersTree->Fill();
852 delete reco;
853}
854
855//_____________________________________________________________________________
856void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
857 const Float_t* const corrADCZN1, const Float_t* const corrADCZP1,
858 const Float_t* const corrADCZN2, const Float_t* const corrADCZP2,
859 const Float_t* const corrADCZEM1, const Float_t* const corrADCZEM2,
860 Float_t* sPMRef1, Float_t* sPMRef2, Bool_t isScalerOn, UInt_t* scaler,
861 Int_t tdcData[32][4], const Int_t* const evQualityBlock,
862 const Int_t* const triggerBlock, const Int_t* const chBlock, UInt_t puBits) const
863{
864 // ****************** Reconstruct one event ******************
865 // ---------------------- Setting reco flags for ESD
866 UInt_t rFlags[32];
867 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
868
869 if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
870 else rFlags[31] = 0x1;
871 //
872 if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
873 if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
874 if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
875
876 if(triggerBlock[0] == 1) rFlags[27] = 0x1;
877 if(triggerBlock[1] == 1) rFlags[26] = 0x1;
878 if(triggerBlock[2] == 1) rFlags[25] = 0x1;
879 if(triggerBlock[3] == 1) rFlags[24] = 0x1;
880
881 if(chBlock[0] == 1) rFlags[18] = 0x1;
882 if(chBlock[1] == 1) rFlags[17] = 0x1;
883 if(chBlock[2] == 1) rFlags[16] = 0x1;
884
885 rFlags[13] = puBits & 0x00000020;
886 rFlags[12] = puBits & 0x00000010;
887 rFlags[11] = puBits & 0x00000080;
888 rFlags[10] = puBits & 0x00000040;
889 rFlags[9] = puBits & 0x00000020;
890 rFlags[8] = puBits & 0x00000010;
891
892 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
893 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
894 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
895 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
896 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
897 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
898
899 UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
900 rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
901 0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
902 0x0 << 19 | rFlags[18] << 18 | rFlags[17] << 17 | rFlags[16] << 16 |
903 0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
904 rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
905 0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
906 rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
907 // --------------------------------------------------
908
909
910 // CH. debug
911/* printf("\n*************************************************\n");
912 printf(" ReconstructEventPbPb -> values after pedestal subtraction:\n");
913 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
914 corrADCZN1[0],corrADCZN1[1],corrADCZN1[2],corrADCZN1[3],corrADCZN1[4]);
915 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
916 corrADCZP1[0],corrADCZP1[1],corrADCZP1[2],corrADCZP1[3],corrADCZP1[4]);
917 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
918 corrADCZN2[0],corrADCZN2[1],corrADCZN2[2],corrADCZN2[3],corrADCZN2[4]);
919 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
920 corrADCZP2[0],corrADCZP2[1],corrADCZP2[2],corrADCZP2[3],corrADCZP2[4]);
921 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",corrADCZEM1[0],corrADCZEM2[0]);
922 printf("*************************************************\n");
923*/
924 // ****** Retrieving calibration data
925 // --- Equalization coefficients ---------------------------------------------
926 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
927 for(Int_t ji=0; ji<5; ji++){
928 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
929 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
930 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
931 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
932 }
933 // --- Energy calibration factors ------------------------------------
934 Float_t calibEne[6];
935 // The energy calibration object already takes into account of E_beam
936 // -> the value from the OCDB can be directly used (Jul 2010)
937 for(Int_t ij=0; ij<6; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
938
939 // ****** Equalization of detector responses
940 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
941 for(Int_t gi=0; gi<10; gi++){
942 if(gi<5){
943 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
944 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
945 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
946 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
947 }
948 else{
949 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
950 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
951 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
952 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
953 }
954 }
955
956 // Ch. debug
957/* printf("\n ------------- EQUALIZATION -------------\n");
958 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
959 equalTowZN1[0],equalTowZN1[1],equalTowZN1[2],equalTowZN1[3],equalTowZN1[4]);
960 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
961 equalTowZP1[0],equalTowZP1[1],equalTowZP1[2],equalTowZP1[3],equalTowZP1[4]);
962 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
963 equalTowZN2[0],equalTowZN2[1],equalTowZN2[2],equalTowZN2[3],equalTowZN2[4]);
964 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
965 equalTowZP2[0],equalTowZP2[1],equalTowZP2[2],equalTowZP2[3],equalTowZP2[4]);
966 printf(" ----------------------------------------\n");
967*/
968
969 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
970 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
971 for(Int_t gi=0; gi<5; gi++){
972 calibSumZN1[0] += equalTowZN1[gi];
973 calibSumZP1[0] += equalTowZP1[gi];
974 calibSumZN2[0] += equalTowZN2[gi];
975 calibSumZP2[0] += equalTowZP2[gi];
976 //
977 calibSumZN1[1] += equalTowZN1[gi+5];
978 calibSumZP1[1] += equalTowZP1[gi+5];
979 calibSumZN2[1] += equalTowZN2[gi+5];
980 calibSumZP2[1] += equalTowZP2[gi+5];
981 }
982 //
983 //fEnCalibData->Print("");
984
985 // High gain chain
986 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]*8.;
987 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]*8.;
988 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]*8.;
989 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]*8.;
990 // Low gain chain
991 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
992 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
993 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
994 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
995 //
996 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
997 calibZEM1[0] = corrADCZEM1[0]*calibEne[4]*8.;
998 calibZEM1[1] = corrADCZEM1[1]*calibEne[4];
999 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]*8.;
1000 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
1001 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
1002
1003 // ****** Energy calibration of detector responses
1004 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
1005 for(Int_t gi=0; gi<5; gi++){
1006 // High gain chain
1007 calibTowZN1[gi] = equalTowZN1[gi]*2*calibEne[0]*8.;
1008 calibTowZP1[gi] = equalTowZP1[gi]*2*calibEne[1]*8.;
1009 calibTowZN2[gi] = equalTowZN2[gi]*2*calibEne[2]*8.;
1010 calibTowZP2[gi] = equalTowZP2[gi]*2*calibEne[3]*8.;
1011 // Low gain chain
1012 calibTowZN1[gi+5] = equalTowZN1[gi+5]*2*calibEne[0];
1013 calibTowZP1[gi+5] = equalTowZP1[gi+5]*2*calibEne[1];
1014 calibTowZN2[gi+5] = equalTowZN2[gi+5]*2*calibEne[2];
1015 calibTowZP2[gi+5] = equalTowZP2[gi+5]*2*calibEne[3];
1016 }
1017
1018 // Ch. debug
1019/* printf("\n ------------- CALIBRATION -------------\n");
1020 printf(" ADCZN1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1021 calibTowZN1[0],calibTowZN1[1],calibTowZN1[2],calibTowZN1[3],calibTowZN1[4]);
1022 printf(" ADCZP1 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1023 calibTowZP1[0],calibTowZP1[1],calibTowZP1[2],calibTowZP1[3],calibTowZP1[4]);
1024 printf(" ADCZN2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1025 calibTowZN2[0],calibTowZN2[1],calibTowZN2[2],calibTowZN2[3],calibTowZN2[4]);
1026 printf(" ADCZP2 [%1.2f %1.2f %1.2f %1.2f %1.2f]\n",
1027 calibTowZP2[0],calibTowZP2[1],calibTowZP2[2],calibTowZP2[3],calibTowZP2[4]);
1028 printf(" ADCZEM1 [%1.2f] ADCZEM2 [%1.2f] \n",calibZEM1[0],calibZEM2[0]);
1029 printf(" ----------------------------------------\n");
1030*/
1031 // ****** Number of detected spectator nucleons
1032 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
1033 if(fBeamEnergy>0.01){
1034 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
1035 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
1036 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
1037 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
1038 }
1039 else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
1040 /*printf("\n\t AliZDCReconstructor -> fBeamEnergy %1.0f: nDetSpecNsideA %d, nDetSpecPsideA %d,"
1041 " nDetSpecNsideC %d, nDetSpecPsideC %d\n",fBeamEnergy,nDetSpecNLeft, nDetSpecPLeft,
1042 nDetSpecNRight, nDetSpecPRight);*/
1043
1044 Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
1045 Int_t nPart=0, nPartA=0, nPartC=0;
1046 Double_t b=0., bA=0., bC=0.;
1047
1048 if(fIsCalibrationMB == kFALSE){
1049 // ****** Reconstruction parameters ------------------
1050 if(!fgRecoParam) fgRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
1051 if(!fgRecoParam){
1052 AliError(" RecoParam object not retrieved correctly: not reconstructing event!!!");
1053 return;
1054 }
1055 TH1D* hNpartDist = fgRecoParam->GethNpartDist();
1056 TH1D* hbDist = fgRecoParam->GethbDist();
1057 Float_t fClkCenter = fgRecoParam->GetClkCenter();
1058 if(!hNpartDist || !hbDist){
1059 AliError("Something wrong in Glauber MC histos got from AliZDCREcoParamPbPb: NO EVENT RECO FOR ZDC DATA!!!\n\n");
1060 return;
1061 }
1062
1063 if(!fgMBCalibData) fgMBCalibData = const_cast<AliZDCMBCalib*>(GetMBCalibData());
1064 TH2F *hZDCvsZEM = fgMBCalibData->GethZDCvsZEM();
1065 TH2F *hZDCCvsZEM = fgMBCalibData->GethZDCCvsZEM();
1066 TH2F *hZDCAvsZEM = fgMBCalibData->GethZDCAvsZEM();
1067 //
1068 Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
1069 Double_t origin = xHighEdge*fClkCenter;
1070 // Ch. debug
1071 //printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
1072 //
1073 // ====> Summed ZDC info (sideA+side C)
1074 TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
1075 Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
1076 Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
1077 line->SetParameter(0, y/(x-origin));
1078 line->SetParameter(1, -origin*y/(x-origin));
1079 // Ch. debug
1080 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
1081 //printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
1082 //
1083 Double_t countPerc=0;
1084 Double_t xBinCenter=0, yBinCenter=0;
1085 for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
1086 for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
1087 xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
1088 yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
1089 //
1090 if(line->GetParameter(0)>0){
1091 if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
1092 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
1093 // Ch. debug
1094 //printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
1095 //xBinCenter, yBinCenter, countPerc);
1096 }
1097 }
1098 else{
1099 if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
1100 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
1101 // Ch. debug
1102 //printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
1103 //xBinCenter, yBinCenter, countPerc);
1104 }
1105 }
1106 }
1107 }
1108 //
1109 Double_t xSecPerc = 0.;
1110 if(hZDCvsZEM->GetEntries()!=0){
1111 xSecPerc = countPerc/hZDCvsZEM->GetEntries();
1112 }
1113 else{
1114 AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
1115 }
1116 // Ch. debug
1117 //printf(" xSecPerc %1.4f \n", xSecPerc);
1118
1119 // ====> side C
1120 TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
1121 Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
1122 lineC->SetParameter(0, yC/(x-origin));
1123 lineC->SetParameter(1, -origin*yC/(x-origin));
1124 // Ch. debug
1125 //printf(" ***************** Side C \n");
1126 //printf(" E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
1127 //
1128 Double_t countPercC=0;
1129 Double_t xBinCenterC=0, yBinCenterC=0;
1130 for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
1131 for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
1132 xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
1133 yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
1134 if(lineC->GetParameter(0)>0){
1135 if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1136 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1137 }
1138 }
1139 else{
1140 if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
1141 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
1142 }
1143 }
1144 }
1145 }
1146 //
1147 Double_t xSecPercC = 0.;
1148 if(hZDCCvsZEM->GetEntries()!=0){
1149 xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
1150 }
1151 else{
1152 AliWarning(" Histogram hZDCCvsZEM from OCDB has no entries!!!");
1153 }
1154 // Ch. debug
1155 //printf(" xSecPercC %1.4f \n", xSecPercC);
1156
1157 // ====> side A
1158 TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
1159 Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
1160 lineA->SetParameter(0, yA/(x-origin));
1161 lineA->SetParameter(1, -origin*yA/(x-origin));
1162 //
1163 // Ch. debug
1164 //printf(" ***************** Side A \n");
1165 //printf(" E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
1166 //
1167 Double_t countPercA=0;
1168 Double_t xBinCenterA=0, yBinCenterA=0;
1169 for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
1170 for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
1171 xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
1172 yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
1173 if(lineA->GetParameter(0)>0){
1174 if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1175 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1176 }
1177 }
1178 else{
1179 if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
1180 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
1181 }
1182 }
1183 }
1184 }
1185 //
1186 Double_t xSecPercA = 0.;
1187 if(hZDCAvsZEM->GetEntries()!=0){
1188 xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
1189 }
1190 else{
1191 AliWarning(" Histogram hZDCAvsZEM from OCDB has no entries!!!");
1192 }
1193 // Ch. debug
1194 //printf(" xSecPercA %1.4f \n", xSecPercA);
1195
1196 // ****** Number of participants (from E_ZDC vs. E_ZEM correlation)
1197 Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1198 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1199 nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1200 if((1.-nPartFrac) < xSecPerc){
1201 nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1202 // Ch. debug
1203 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
1204 //printf(" nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1205 break;
1206 }
1207 }
1208 if(nPart<0) nPart=0;
1209 //
1210 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1211 nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1212 if((1.-nPartFracC) < xSecPercC){
1213 nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1214 // Ch. debug
1215 //printf(" ***************** Side C \n");
1216 //printf(" nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1217 break;
1218 }
1219 }
1220 if(nPartC<0) nPartC=0;
1221 //
1222 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1223 nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1224 if((1.-nPartFracA) < xSecPercA){
1225 nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1226 // Ch. debug
1227 //printf(" ***************** Side A \n");
1228 //printf(" nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1229 break;
1230 }
1231 }
1232 if(nPartA<0) nPartA=0;
1233
1234 // ****** Impact parameter (from E_ZDC vs. E_ZEM correlation)
1235 Double_t bFrac=0., bFracC=0., bFracA=0.;
1236 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1237 bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1238 if(bFrac > xSecPerc){
1239 b = hbDist->GetBinLowEdge(ibbin);
1240 break;
1241 }
1242 }
1243 //
1244 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1245 bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1246 if(bFracC > xSecPercC){
1247 bC = hbDist->GetBinLowEdge(ibbin);
1248 break;
1249 }
1250 }
1251 //
1252 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1253 bFracA += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1254 if(bFracA > xSecPercA){
1255 bA = hbDist->GetBinLowEdge(ibbin);
1256 break;
1257 }
1258 }
1259
1260 // ****** Number of spectator nucleons
1261 nGenSpec = 416 - nPart;
1262 nGenSpecC = 416 - nPartC;
1263 nGenSpecA = 416 - nPartA;
1264 if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1265 if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1266 if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1267
1268 delete line;
1269 delete lineC; delete lineA;
1270
1271 } // ONLY IF fIsCalibrationMB==kFALSE
1272
1273 Bool_t energyFlag = kTRUE;
1274 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1275 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1276 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1277 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1278 nGenSpec, nGenSpecA, nGenSpecC,
1279 nPart, nPartA, nPartC, b, bA, bC,
1280 recoFlag, energyFlag, isScalerOn, scaler, tdcData);
1281
1282 const Int_t kBufferSize = 4000;
1283 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1284 //reco->Print("");
1285 // write the output tree
1286 clustersTree->Fill();
1287 delete reco;
1288}
1289
1290
1291//_____________________________________________________________________________
1292void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1293{
1294 // fill energies and number of participants to the ESD
1295
1296 AliZDCReco reco;
1297 AliZDCReco* preco = &reco;
1298 clustersTree->SetBranchAddress("ZDC", &preco);
1299 clustersTree->GetEntry(0);
1300 //
1301 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1302 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1303 for(Int_t i=0; i<5; i++){
1304 tZN1Ene[i] = reco.GetZN1HREnTow(i);
1305 tZN2Ene[i] = reco.GetZN2HREnTow(i);
1306 tZP1Ene[i] = reco.GetZP1HREnTow(i);
1307 tZP2Ene[i] = reco.GetZP2HREnTow(i);
1308 //
1309 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1310 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1311 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1312 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1313 }
1314 //
1315 fESDZDC->SetZN1TowerEnergy(tZN1Ene);
1316 fESDZDC->SetZN2TowerEnergy(tZN2Ene);
1317 fESDZDC->SetZP1TowerEnergy(tZP1Ene);
1318 fESDZDC->SetZP2TowerEnergy(tZP2Ene);
1319 //
1320 fESDZDC->SetZN1TowerEnergyLR(tZN1EneLR);
1321 fESDZDC->SetZN2TowerEnergyLR(tZN2EneLR);
1322 fESDZDC->SetZP1TowerEnergyLR(tZP1EneLR);
1323 fESDZDC->SetZP2TowerEnergyLR(tZP2EneLR);
1324 //
1325 Int_t nPart = reco.GetNParticipants();
1326 Int_t nPartA = reco.GetNPartSideA();
1327 Int_t nPartC = reco.GetNPartSideC();
1328 Double_t b = reco.GetImpParameter();
1329 Double_t bA = reco.GetImpParSideA();
1330 Double_t bC = reco.GetImpParSideC();
1331 UInt_t recoFlag = reco.GetRecoFlag();
1332
1333 fESDZDC->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(),
1334 reco.GetZEM1HRsignal(), reco.GetZEM2HRsignal(),
1335 reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1336 nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1337
1338 // Writing ZDC scaler for cross section calculation
1339 // ONLY IF the scaler has been read during the event
1340 if(reco.IsScalerOn()==kTRUE){
1341 UInt_t counts[32];
1342 for(Int_t jk=0; jk<32; jk++) counts[jk] = reco.GetZDCScaler(jk);
1343 fESDZDC->SetZDCScaler(counts);
1344 }
1345
1346 // Writing TDC data into ZDC ESDs
1347 Int_t tdcValues[32][4];
1348 Float_t tdcCorrected[32][4];
1349 for(Int_t jk=0; jk<32; jk++){
1350 for(Int_t lk=0; lk<4; lk++){
1351 tdcValues[jk][lk] = reco.GetZDCTDCData(jk, lk);
1352 }
1353 }
1354 // 4/2/2011 -> Subtracting L0 (tdcValues[15]) instead of ADC gate
1355 // we try to keep the TDC oscillations as low as possible!
1356 for(Int_t jk=0; jk<32; jk++){
1357 for(Int_t lk=0; lk<4; lk++){
1358 if(tdcValues[jk][lk]!=0.) tdcCorrected[jk][lk] = 0.025*(tdcValues[jk][lk]-tdcValues[15][0])+fMeanPhase;
1359 }
1360 }
1361 fESDZDC->SetZDCTDCData(tdcValues);
1362 fESDZDC->SetZDCTDCCorrected(tdcCorrected);
1363 fESDZDC->AliESDZDC::SetBit(AliESDZDC::kCorrectedTDCFilled, reco.GetEnergyFlag());
1364 fESDZDC->AliESDZDC::SetBit(AliESDZDC::kEnergyCalibratedSignal, kTRUE);
1365
1366 if(esd) esd->SetZDCData(fESDZDC);
1367}
1368
1369//_____________________________________________________________________________
1370AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1371{
1372 // Setting the storage
1373
1374 Bool_t deleteManager = kFALSE;
1375
1376 AliCDBManager *manager = AliCDBManager::Instance();
1377 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1378
1379 if(!defstorage || !(defstorage->Contains("ZDC"))){
1380 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1381 manager->SetDefaultStorage(uri);
1382 deleteManager = kTRUE;
1383 }
1384
1385 AliCDBStorage *storage = manager->GetDefaultStorage();
1386
1387 if(deleteManager){
1388 AliCDBManager::Instance()->UnsetDefaultStorage();
1389 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1390 }
1391
1392 return storage;
1393}
1394
1395//_____________________________________________________________________________
1396AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1397{
1398
1399 // Getting pedestal calibration object for ZDC set
1400
1401 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1402 if(!entry) AliFatal("No calibration data loaded!");
1403 entry->SetOwner(kFALSE);
1404
1405 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
1406 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1407
1408 return calibdata;
1409}
1410
1411//_____________________________________________________________________________
1412AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1413{
1414
1415 // Getting energy and equalization calibration object for ZDC set
1416
1417 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1418 if(!entry) AliFatal("No calibration data loaded!");
1419 entry->SetOwner(kFALSE);
1420
1421 AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1422 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1423
1424 return calibdata;
1425}
1426
1427//_____________________________________________________________________________
1428AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1429{
1430
1431 // Getting energy and equalization calibration object for ZDC set
1432
1433 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1434 if(!entry) AliFatal("No calibration data loaded!");
1435 entry->SetOwner(kFALSE);
1436
1437 AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1438 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1439
1440 return calibdata;
1441}
1442
1443//_____________________________________________________________________________
1444AliZDCMBCalib* AliZDCReconstructor::GetMBCalibData() const
1445{
1446
1447 // Getting energy and equalization calibration object for ZDC set
1448
1449 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/MBCalib");
1450 if(!entry) AliFatal("No calibration data loaded!");
1451 entry->SetOwner(kFALSE);
1452
1453 AliZDCMBCalib *calibdata = dynamic_cast<AliZDCMBCalib*> (entry->GetObject());
1454 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1455
1456 return calibdata;
1457}