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