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