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