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