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