1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // ************** Class for ZDC reconstruction ************** //
21 // Author: Chiara.Oppedisano@to.infn.it //
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) //
27 ///////////////////////////////////////////////////////////////////////////////
35 #include "AliRawReader.h"
36 #include "AliGRPObject.h"
37 #include "AliESDEvent.h"
38 #include "AliESDZDC.h"
39 #include "AliZDCDigit.h"
40 #include "AliZDCRawStream.h"
41 #include "AliZDCReco.h"
42 #include "AliZDCReconstructor.h"
43 #include "AliZDCPedestals.h"
44 #include "AliZDCEnCalib.h"
45 #include "AliZDCTowerCalib.h"
46 #include "AliZDCRecoParam.h"
47 #include "AliZDCRecoParampp.h"
48 #include "AliZDCRecoParamPbPb.h"
51 ClassImp(AliZDCReconstructor)
52 AliZDCRecoParam *AliZDCReconstructor::fRecoParam=0; //reconstruction parameters
54 //_____________________________________________________________________________
55 AliZDCReconstructor:: AliZDCReconstructor() :
56 fPedData(GetPedestalData()),
57 fEnCalibData(GetEnergyCalibData()),
58 fTowCalibData(GetTowerCalibData()),
62 fIsCalibrationMB(kFALSE),
66 // **** Default constructor
71 //_____________________________________________________________________________
72 AliZDCReconstructor::~AliZDCReconstructor()
75 if(fRecoParam) delete fRecoParam;
76 if(fPedData) delete fPedData;
77 if(fEnCalibData) delete fEnCalibData;
78 if(fTowCalibData) delete fTowCalibData;
81 //____________________________________________________________________________
82 void AliZDCReconstructor::Init()
84 // Setting reconstruction mode
85 // Getting beam type and beam energy from GRP calibration object
87 if(fRecoMode==0 && fBeamEnergy==0.){
88 // Initialization of the GRP entry
89 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
90 AliGRPObject* grpData = 0x0;
92 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
95 grpData = new AliGRPObject();
96 grpData->ReadValuesFromMap(m);
99 grpData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
102 AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
104 if(!grpData) AliError("No GRP entry found in OCDB!");
106 TString runType = grpData->GetRunType();
107 if(runType==AliGRPObject::GetInvalidString()){
108 AliWarning("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN");
111 if((runType.CompareTo("CALIBRATION_MB")) == 0){
112 fIsCalibrationMB = kTRUE;
115 TString beamType = grpData->GetBeamType();
116 // This is a temporary solution to allow reconstruction in tests without beam
117 if(((beamType.CompareTo("UNKNOWN"))==0) && ((runType.CompareTo("PHYSICS")) == 0)){
120 else if(beamType==AliGRPObject::GetInvalidString()){
121 AliWarning("GRP/GRP/Data entry: missing value for the beam type !");
122 AliError("\t ZDC does not reconstruct event 4 UNKNOWN beam type\n");
126 if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
127 ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0)) fRecoMode=1;
128 else if((beamType.CompareTo("A-A")) == 0){
130 if(fIsCalibrationMB == kTRUE){
131 fRecoParam = new AliZDCRecoParamPbPb();
133 TH2F* hZDCvsZEM = new TH2F("hZDCvsZEM","hZDCvsZEM",100,0.,10.,100,0.,1000.);
134 hZDCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCvsZEM->SetYTitle("E_{ZDC} (TeV)");
135 fRecoParam->SetZDCvsZEM(hZDCvsZEM);
137 TH2F* hZDCCvsZEM = new TH2F("hZDCCvsZEM","hZDCCvsZEM",100,0.,10.,100,0.,500.);
138 hZDCCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCCvsZEM->SetYTitle("E_{ZDCC} (TeV)");
139 fRecoParam->SetZDCCvsZEM(hZDCCvsZEM);
141 TH2F* hZDCAvsZEM = new TH2F("hZDCAvsZEM","hZDCAvsZEM",100,0.,10.,100,0.,500.);
142 hZDCAvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCAvsZEM->SetYTitle("E_{ZDCA} (TeV)");
143 fRecoParam->SetZDCAvsZEM(hZDCAvsZEM);
145 AliInfo("\n ***** CALIBRATION_MB data -> building AliZDCRecoParamPbPb object *****");
149 fBeamEnergy = grpData->GetBeamEnergy();
150 if(fBeamEnergy==AliGRPObject::GetInvalidFloat()){
151 AliWarning("GRP/GRP/Data entry: missing value for the beam energy ! Using 0.");
155 if(fIsCalibrationMB==kFALSE)
156 printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f GeV *****\n\n",beamType.Data(), fBeamEnergy);
159 AliError(" ATTENTION!!!!!! No beam type nor beam energy has been set!!!!!!\n");
164 //_____________________________________________________________________________
165 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
167 // *** Local ZDC reconstruction for digits
168 // Works on the current event
170 // Retrieving calibration data
171 // Parameters for mean value pedestal subtraction
173 Float_t meanPed[2*kNch];
174 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
175 // Parameters pedestal subtraction through correlation with out-of-time signals
176 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
177 for(Int_t jj=0; jj<2*kNch; jj++){
178 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
179 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
184 AliZDCDigit* pdigit = &digit;
185 digitsTree->SetBranchAddress("ZDC", &pdigit);
186 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
189 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
190 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
191 for(Int_t i=0; i<10; i++){
192 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
193 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
196 Int_t digNentries = digitsTree->GetEntries();
197 Float_t ootDigi[kNch];
198 // -- Reading out-of-time signals (last kNch entries) for current event
200 for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
201 ootDigi[iDigit] = digitsTree->GetEntry(iDigit);
205 for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
206 digitsTree->GetEntry(iDigit);
207 if (!pdigit) continue;
209 Int_t det = digit.GetSector(0);
210 Int_t quad = digit.GetSector(1);
212 Float_t ped2SubHg=0., ped2SubLg=0.;
214 if(det==1) pedindex = quad;
215 else if(det==2) pedindex = quad+5;
216 else if(det==3) pedindex = quad+9;
217 else if(det==4) pedindex = quad+12;
218 else if(det==5) pedindex = quad+17;
220 else pedindex = (det-1)/3+22;
223 ped2SubHg = meanPed[pedindex];
224 ped2SubLg = meanPed[pedindex+kNch];
226 else if(fPedSubMode==1){
227 ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
228 ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
231 if(quad != 5){ // ZDC (not reference PTMs!)
232 if(det == 1){ // *** ZNC
233 tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
234 tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
235 if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
236 if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
238 //printf("\t pedindex %d tZN1Corr[%d] = %1.0f tZN1Corr[%d] = %1.0f",
239 // pedindex, quad, tZN1Corr[quad], quad+5, tZN1Corr[quad+5]);
241 else if(det == 2){ // *** ZP1
242 tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
243 tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
244 if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
245 if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
247 //printf("\t pedindex %d tZP1Corr[%d] = %1.0f tZP1Corr[%d] = %1.0f",
248 // pedindex, quad, tZP1Corr[quad], quad+5, tZP1Corr[quad+5]);
251 if(quad == 1){ // *** ZEM1
252 dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
253 dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
254 if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
255 if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
257 //printf("\t pedindex %d tZEM1Corr[%d] = %1.0f tZEM1Corr[%d] = %1.0f",
258 // pedindex, quad, tZEM1Corr[quad], quad+1, tZEM1Corr[quad+1]);
260 else if(quad == 2){ // *** ZEM2
261 dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
262 dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
263 if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
264 if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
266 //printf("\t pedindex %d tZEM2Corr[%d] = %1.0f tZEM2Corr[%d] = %1.0f",
267 // pedindex, quad, tZEM2Corr[quad], quad+1, tZEM2Corr[quad+1]);
270 else if(det == 4){ // *** ZN2
271 tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
272 tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
273 if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
274 if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
276 //printf("\t pedindex %d tZN2Corr[%d] = %1.0f tZN2Corr[%d] = %1.0f\n",
277 // pedindex, quad, tZN2Corr[quad], quad+5, tZN2Corr[quad+5]);
279 else if(det == 5){ // *** ZP2
280 tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
281 tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
282 if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
283 if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
285 //printf("\t pedindex %d tZP2Corr[%d] = %1.0f tZP2Corr[%d] = %1.0f\n",
286 // pedindex, quad, tZP2Corr[quad], quad+5, tZP2Corr[quad+5]);
289 else{ // Reference PMs
291 sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
292 sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
294 if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
295 if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
298 sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
299 sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
301 if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
302 if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
307 /*printf(" - AliZDCReconstructor -> digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
308 iDigit, det, quad, ped2SubHg, ped2SubLg);
309 printf(" HGChain -> RawDig %d DigCorr %1.2f\n", digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg);
310 printf(" LGChain -> RawDig %d DigCorr %1.2f\n", digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);
315 // If CALIBRATION_MB run build the RecoParam object
316 if(fIsCalibrationMB){
317 Float_t ZDCC=0., ZDCA=0., ZEM=0;
318 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
319 for(Int_t jkl=0; jkl<5; jkl++){
320 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
321 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
323 // Using energies in TeV in fRecoParam object
324 BuildRecoParam(ZDCC/1000., ZDCA/1000., ZEM/1000.);
326 Bool_t recFlag1 = kFALSE, recFlag2 = kFALSE, recFlag3 = kFALSE;
327 // reconstruct the event
329 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
330 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, recFlag1, recFlag2, recFlag3);
331 else if(fRecoMode==2)
332 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
333 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, recFlag1, recFlag2, recFlag3);
336 //_____________________________________________________________________________
337 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
339 // *** ZDC raw data reconstruction
340 // Works on the current event
342 // Retrieving calibration data
343 // Parameters for pedestal subtraction
345 Float_t meanPed[2*kNch];
346 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
347 // Parameters pedestal subtraction through correlation with out-of-time signals
348 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
349 for(Int_t jj=0; jj<2*kNch; jj++){
350 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
351 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
354 Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
355 Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
356 Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
357 Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
358 Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
359 Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
360 for(Int_t ich=0; ich<5; ich++){
361 adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
362 adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
363 adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
364 adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
366 adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
367 pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
371 // loop over raw data
372 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
373 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
374 for(Int_t i=0; i<10; i++){
375 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
376 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
379 //fNRun = (Int_t) rawReader->GetRunNumber();
380 Bool_t chOff=kFALSE, isUndflw=kFALSE, isOvflw=kFALSE;
383 AliZDCRawStream rawData(rawReader);
384 while(rawData.Next()){
385 if(rawData.IsCalibration() == kFALSE){ // Reading ADCs
386 //printf(" **** Reading ADC raw data **** \n");
388 if(rawData.GetNChannelsOn() < 48 ) chOff=kTRUE;
389 if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)) isUndflw=kTRUE;
390 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)) isOvflw=kTRUE;
392 if((rawData.IsADCDataWord()) && (isUndflw==kFALSE) && (isOvflw==kFALSE)){
394 Int_t adcMod = rawData.GetADCModule();
395 Int_t det = rawData.GetSector(0);
396 Int_t quad = rawData.GetSector(1);
397 Int_t gain = rawData.GetADCGain();
400 // Mean pedestal value subtraction -------------------------------------------------------
401 if(fPedSubMode == 0){
402 // Not interested in o.o.t. signals (ADC modules 2, 3)
403 if(adcMod == 2 || adcMod == 3) continue;
405 if(quad != 5){ // ZDCs (not reference PTMs)
408 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
409 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
413 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
414 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
419 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
420 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
423 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
424 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
429 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
430 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
434 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
435 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
438 else{ // reference PM
439 pedindex = (det-1)/3 + 22;
441 if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
442 else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
445 if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
446 else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
450 /*printf(" -> AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f\n",
451 det,quad,gain, pedindex, meanPed[pedindex]);
452 printf(" -> AliZDCReconstructor: RawADC %1.0f ADCCorr %1.0f\n",
453 rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);*/
455 }// mean pedestal subtraction
456 // Pedestal subtraction from correlation ------------------------------------------------
457 else if(fPedSubMode == 1){
459 if(adcMod==0 || adcMod==1){
460 if(quad != 5){ // signals from ZDCs
462 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
463 else adcZN1lg[quad] = rawData.GetADCValue();
466 if(gain==0) adcZP1[quad] = rawData.GetADCValue();
467 else adcZP1lg[quad] = rawData.GetADCValue();
470 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
471 else adcZEMlg[quad-1] = rawData.GetADCValue();
474 if(gain==0) adcZN2[quad] = rawData.GetADCValue();
475 else adcZN2lg[quad] = rawData.GetADCValue();
478 if(gain==0) adcZP2[quad] = rawData.GetADCValue();
479 else adcZP2lg[quad] = rawData.GetADCValue();
482 else{ // signals from reference PM
483 if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
484 else pmReflg[quad-1] = rawData.GetADCValue();
487 // Out-of-time pedestals
488 else if(adcMod==2 || adcMod==3){
489 if(quad != 5){ // signals from ZDCs
491 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
492 else adcZN1ootlg[quad] = rawData.GetADCValue();
495 if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
496 else adcZP1ootlg[quad] = rawData.GetADCValue();
499 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
500 else adcZEMootlg[quad-1] = rawData.GetADCValue();
503 if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
504 else adcZN2ootlg[quad] = rawData.GetADCValue();
507 if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
508 else adcZP2ootlg[quad] = rawData.GetADCValue();
511 else{ // signals from reference PM
512 if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
513 else pmRefootlg[quad-1] = rawData.GetADCValue();
516 } // pedestal subtraction from correlation
518 //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
519 // det,quad,gain, pedindex, meanPed[pedindex]);
521 }// Not raw data from calibration run!
527 for(Int_t t=0; t<5; t++){
528 tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
529 tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
531 tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
532 tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
534 tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
535 tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
537 tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
538 tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
539 // 0---------0 Ch. debug 0---------0
540 /* printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
541 printf("\tCorrCoeff0\tCorrCoeff1\n");
542 printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
543 printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
544 printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
545 printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
546 printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
547 printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
548 printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
549 printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
551 printf("ZN1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
552 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
553 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
554 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
556 printf("ZP1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
557 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
558 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
559 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
561 printf("ZN2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
562 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
563 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
564 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
566 printf("ZP2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
567 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
568 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
569 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
572 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
573 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
574 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
575 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
577 sPMRef1[0] = pmRef[0] - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
578 sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
579 sPMRef2[0] = pmRef[0] - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
580 sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
583 // If CALIBRATION_MB run build the RecoParam object
584 if(fIsCalibrationMB){
585 Float_t ZDCC=0., ZDCA=0., ZEM=0;
586 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
587 for(Int_t jkl=0; jkl<5; jkl++){
588 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
589 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
591 BuildRecoParam(ZDCC/100., ZDCA/100., ZEM/100.);
593 // reconstruct the event
595 if(fRecoMode==1) // p-p data
596 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
597 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, chOff, isUndflw, isOvflw);
598 else if(fRecoMode==2) // Pb-Pb data
599 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
600 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, chOff, isUndflw, isOvflw);
604 //_____________________________________________________________________________
605 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* corrADCZN1,
606 Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
607 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
608 Bool_t channelsOff, Bool_t chUnderflow, Bool_t chOverflow) const
610 // ****************** Reconstruct one event ******************
612 // ---- Setting reco flags
615 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
617 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
618 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
619 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
620 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
621 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
622 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
624 if(channelsOff==kTRUE) rFlags[8] = 0x1;
625 if(chUnderflow == kTRUE) rFlags[9] = 0x1;
626 if(chOverflow==kTRUE) rFlags[10] = 0x1;
627 recoFlag = rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
628 rFlags[5] << 5 | rFlags[4] << 4 | rFlags[3] << 3 |
629 rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
631 // ****** Retrieving calibration data
632 // --- Equalization coefficients ---------------------------------------------
633 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
634 for(Int_t ji=0; ji<5; ji++){
635 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
636 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
637 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
638 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
640 // --- Energy calibration factors ------------------------------------
642 // **** Energy calibration coefficient set to 1
643 // **** (no trivial way to calibrate in p-p runs)
644 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
646 // ****** Equalization of detector responses
647 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
648 for(Int_t gi=0; gi<10; gi++){
650 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
651 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
652 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
653 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
656 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
657 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
658 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
659 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
663 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
664 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
665 for(Int_t gi=0; gi<5; gi++){
666 calibSumZN1[0] += equalTowZN1[gi];
667 calibSumZP1[0] += equalTowZP1[gi];
668 calibSumZN2[0] += equalTowZN2[gi];
669 calibSumZP2[0] += equalTowZP2[gi];
671 calibSumZN1[1] += equalTowZN1[gi+5];
672 calibSumZP1[1] += equalTowZP1[gi+5];
673 calibSumZN2[1] += equalTowZN2[gi+5];
674 calibSumZP2[1] += equalTowZP2[gi+5];
677 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
678 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
679 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
680 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
682 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
683 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
684 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
685 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
687 // ****** Energy calibration of detector responses
688 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
689 for(Int_t gi=0; gi<5; gi++){
691 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
692 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
693 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
694 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
696 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
697 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
698 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
699 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
702 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
703 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
704 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
705 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
706 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
707 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
709 // ****** No. of spectator and participants nucleons
710 // Variables calculated to comply with ESD structure
711 // *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
712 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
713 Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
714 Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
715 Double_t impPar=0., impPar1=0., impPar2=0.;
717 // create the output tree
718 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
719 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
720 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
721 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
722 nGenSpec, nGenSpecLeft, nGenSpecRight,
723 nPart, nPartTotLeft, nPartTotRight,
724 impPar, impPar1, impPar2,
727 const Int_t kBufferSize = 4000;
728 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
729 // write the output tree
730 clustersTree->Fill();
733 //_____________________________________________________________________________
734 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
735 Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
736 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
737 Bool_t channelsOff, Bool_t chUnderflow, Bool_t chOverflow) const
739 // ****************** Reconstruct one event ******************
741 // ---- Setting reco flags
744 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
746 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
747 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
748 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
749 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
750 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
751 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
753 if(channelsOff==kTRUE) rFlags[8] = 0x1;
754 if(chUnderflow == kTRUE) rFlags[9] = 0x1;
755 if(chOverflow==kTRUE) rFlags[10] = 0x1;
756 recoFlag = rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
757 rFlags[5] << 5 | rFlags[4] << 4 | rFlags[3] << 3 |
758 rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
760 // ****** Retrieving calibration data
761 // --- Equalization coefficients ---------------------------------------------
762 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
763 for(Int_t ji=0; ji<5; ji++){
764 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
765 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
766 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
767 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
769 // --- Energy calibration factors ------------------------------------
770 Float_t valFromOCDB[6], calibEne[6];
771 for(Int_t ij=0; ij<6; ij++){
772 valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
774 if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
775 else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
777 else calibEne[ij] = valFromOCDB[ij];
780 // ****** Equalization of detector responses
781 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
782 for(Int_t gi=0; gi<10; gi++){
784 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
785 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
786 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
787 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
790 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
791 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
792 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
793 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
797 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
798 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
799 for(Int_t gi=0; gi<5; gi++){
800 calibSumZN1[0] += equalTowZN1[gi];
801 calibSumZP1[0] += equalTowZP1[gi];
802 calibSumZN2[0] += equalTowZN2[gi];
803 calibSumZP2[0] += equalTowZP2[gi];
805 calibSumZN1[1] += equalTowZN1[gi+5];
806 calibSumZP1[1] += equalTowZP1[gi+5];
807 calibSumZN2[1] += equalTowZN2[gi+5];
808 calibSumZP2[1] += equalTowZP2[gi+5];
811 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
812 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
813 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
814 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
816 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
817 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
818 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
819 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
821 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
822 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
823 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
824 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
825 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
826 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
828 // ****** Energy calibration of detector responses
829 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
830 for(Int_t gi=0; gi<5; gi++){
832 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
833 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
834 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
835 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
837 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
838 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
839 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
840 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
843 // ****** Number of detected spectator nucleons
844 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
846 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
847 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
848 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
849 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
851 else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
852 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
853 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
854 nDetSpecNRight, nDetSpecPRight);*/
856 Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
857 Int_t nPart=0, nPartA=0, nPartC=0;
858 Double_t b=0., bA=0., bC=0.;
860 if(fIsCalibrationMB == kFALSE){
861 // ****** Reconstruction parameters ------------------
863 //fRecoParam->Print("");
865 TH2F *hZDCvsZEM = fRecoParam->GethZDCvsZEM();
866 TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
867 TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
868 TH1D *hNpartDist = fRecoParam->GethNpartDist();
869 TH1D *hbDist = fRecoParam->GethbDist();
870 Float_t ClkCenter = fRecoParam->GetClkCenter();
872 Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
873 Double_t origin = xHighEdge*ClkCenter;
875 //printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
877 // ====> Summed ZDC info (sideA+side C)
878 TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
879 Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
880 Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
881 line->SetParameter(0, y/(x-origin));
882 line->SetParameter(1, -origin*y/(x-origin));
884 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
885 //printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
887 Double_t countPerc=0;
888 Double_t xBinCenter=0, yBinCenter=0;
889 for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
890 for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
891 xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
892 yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
894 if(line->GetParameter(0)>0){
895 if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
896 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
898 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
899 xBinCenter, yBinCenter, countPerc);*/
903 if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
904 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
906 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
907 xBinCenter, yBinCenter, countPerc);*/
913 Double_t xSecPerc = 0.;
914 if(hZDCvsZEM->GetEntries()!=0){
915 xSecPerc = countPerc/hZDCvsZEM->GetEntries();
918 AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
921 //printf(" xSecPerc %1.4f \n", xSecPerc);
924 TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
925 Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
926 lineC->SetParameter(0, yC/(x-origin));
927 lineC->SetParameter(1, -origin*yC/(x-origin));
929 //printf(" ***************** Side C \n");
930 //printf(" E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
932 Double_t countPercC=0;
933 Double_t xBinCenterC=0, yBinCenterC=0;
934 for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
935 for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
936 xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
937 yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
938 if(lineC->GetParameter(0)>0){
939 if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
940 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
944 if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
945 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
951 Double_t xSecPercC = 0.;
952 if(hZDCCvsZEM->GetEntries()!=0){
953 xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
956 AliWarning(" Histogram hZDCCvsZEM from OCDB has no entries!!!");
959 //printf(" xSecPercC %1.4f \n", xSecPercC);
962 TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
963 Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
964 lineA->SetParameter(0, yA/(x-origin));
965 lineA->SetParameter(1, -origin*yA/(x-origin));
968 //printf(" ***************** Side A \n");
969 //printf(" E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
971 Double_t countPercA=0;
972 Double_t xBinCenterA=0, yBinCenterA=0;
973 for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
974 for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
975 xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
976 yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
977 if(lineA->GetParameter(0)>0){
978 if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
979 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
983 if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
984 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
990 Double_t xSecPercA = 0.;
991 if(hZDCAvsZEM->GetEntries()!=0){
992 xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
995 AliWarning(" Histogram hZDCAvsZEM from OCDB has no entries!!!");
998 //printf(" xSecPercA %1.4f \n", xSecPercA);
1000 // ****** Number of participants (from E_ZDC vs. E_ZEM correlation)
1001 Int_t nPart=0, nPartC=0, nPartA=0;
1002 Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
1003 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1004 nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1005 if((1.-nPartFrac) < xSecPerc){
1006 nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1008 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
1009 //printf(" nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1013 if(nPart<0) nPart=0;
1015 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1016 nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1017 if((1.-nPartFracC) < xSecPercC){
1018 nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1020 //printf(" ***************** Side C \n");
1021 //printf(" nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1025 if(nPartC<0) nPartC=0;
1027 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1028 nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1029 if((1.-nPartFracA) < xSecPercA){
1030 nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1032 //printf(" ***************** Side A \n");
1033 //printf(" nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1037 if(nPartA<0) nPartA=0;
1039 // ****** Impact parameter (from E_ZDC vs. E_ZEM correlation)
1040 Float_t b=0, bC=0, bA=0;
1041 Double_t bFrac=0., bFracC=0., bFracA=0.;
1042 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1043 bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1044 if((1.-bFrac) < xSecPerc){
1045 b = hbDist->GetBinLowEdge(ibbin);
1050 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1051 bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1052 if((1.-bFracC) < xSecPercC){
1053 bC = hbDist->GetBinLowEdge(ibbin);
1058 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1059 bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1060 if((1.-bFracA) < xSecPercA){
1061 bA = hbDist->GetBinLowEdge(ibbin);
1066 // ****** Number of spectator nucleons
1067 Int_t nGenSpec=0, nGenSpecC=0, nGenSpecA=0;
1069 nGenSpec = 416 - nPart;
1070 nGenSpecC = 416 - nPartC;
1071 nGenSpecA = 416 - nPartA;
1072 if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1073 if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1074 if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1077 /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1078 " calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
1079 calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1080 printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n",
1081 nGenSpecLeft, nGenSpecRight);
1082 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1085 delete lineC; delete lineA;
1087 } // ONLY IF fIsCalibrationMB==kFALSE
1089 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1090 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1091 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1092 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1093 nGenSpec, nGenSpecA, nGenSpecC,
1094 nPart, nPartA, nPartC, b, bA, bC,
1097 const Int_t kBufferSize = 4000;
1098 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1100 // write the output tree
1101 clustersTree->Fill();
1104 //_____________________________________________________________________________
1105 void AliZDCReconstructor::BuildRecoParam(Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1107 // Calculate RecoParam object for Pb-Pb data
1108 (fRecoParam->GethZDCvsZEM())->Fill(ZDCC+ZDCA, ZEM);
1109 (fRecoParam->GethZDCCvsZEM())->Fill(ZDCC, ZEM);
1110 (fRecoParam->GethZDCAvsZEM())->Fill(ZDCA, ZEM);
1114 //_____________________________________________________________________________
1115 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1117 // fill energies and number of participants to the ESD
1119 if(fIsCalibrationMB==kTRUE) WritePbPbRecoParamInOCDB();
1122 AliZDCReco* preco = &reco;
1123 clustersTree->SetBranchAddress("ZDC", &preco);
1125 clustersTree->GetEntry(0);
1127 AliESDZDC * esdzdc = esd->GetESDZDC();
1128 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1129 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1130 for(Int_t i=0; i<5; i++){
1131 tZN1Ene[i] = reco.GetZN1HREnTow(i);
1132 tZN2Ene[i] = reco.GetZN2HREnTow(i);
1133 tZP1Ene[i] = reco.GetZP1HREnTow(i);
1134 tZP2Ene[i] = reco.GetZP2HREnTow(i);
1136 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1137 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1138 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1139 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1142 esdzdc->SetZN1TowerEnergy(tZN1Ene);
1143 esdzdc->SetZN2TowerEnergy(tZN2Ene);
1144 esdzdc->SetZP1TowerEnergy(tZP1Ene);
1145 esdzdc->SetZP2TowerEnergy(tZP2Ene);
1147 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1148 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1149 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1150 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1152 Int_t nPart = reco.GetNParticipants();
1153 Int_t nPartA = reco.GetNPartSideA();
1154 Int_t nPartC = reco.GetNPartSideC();
1155 Double_t b = reco.GetImpParameter();
1156 Double_t bA = reco.GetImpParSideA();
1157 Double_t bC = reco.GetImpParSideC();
1158 UInt_t recoFlag = reco.GetRecoFlag();
1160 esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
1161 reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1162 nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1166 //_____________________________________________________________________________
1167 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1169 // Setting the storage
1171 Bool_t deleteManager = kFALSE;
1173 AliCDBManager *manager = AliCDBManager::Instance();
1174 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1176 if(!defstorage || !(defstorage->Contains("ZDC"))){
1177 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1178 manager->SetDefaultStorage(uri);
1179 deleteManager = kTRUE;
1182 AliCDBStorage *storage = manager->GetDefaultStorage();
1185 AliCDBManager::Instance()->UnsetDefaultStorage();
1186 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1192 //_____________________________________________________________________________
1193 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1196 // Getting pedestal calibration object for ZDC set
1198 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1199 if(!entry) AliFatal("No calibration data loaded!");
1201 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
1202 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1207 //_____________________________________________________________________________
1208 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1211 // Getting energy and equalization calibration object for ZDC set
1213 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1214 if(!entry) AliFatal("No calibration data loaded!");
1216 AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1217 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1222 //_____________________________________________________________________________
1223 AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1226 // Getting energy and equalization calibration object for ZDC set
1228 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1229 if(!entry) AliFatal("No calibration data loaded!");
1231 AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1232 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1237 //_____________________________________________________________________________
1238 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1241 // Getting reconstruction parameters from OCDB
1243 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1244 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1246 AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1247 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1253 //_____________________________________________________________________________
1254 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1257 // Getting reconstruction parameters from OCDB
1259 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1260 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1262 AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1263 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1269 //_____________________________________________________________________________
1270 void AliZDCReconstructor::WritePbPbRecoParamInOCDB() const
1273 // Writing Pb-Pb reconstruction parameters from OCDB
1275 AliCDBManager *man = AliCDBManager::Instance();
1276 AliCDBMetaData *md= new AliCDBMetaData();
1277 md->SetResponsible("Chiara Oppedisano");
1278 md->SetComment("ZDC Pb-Pb reconstruction parameters");
1279 md->SetObjectClassName("AliZDCRecoParamPbPb");
1280 AliCDBId id("ZDC/Calib/RecoParamPbPb",fNRun,AliCDBRunRange::Infinity());
1281 man->Put(fRecoParam, id, md);