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