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