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 if(beamType==AliGRPObject::GetInvalidString()){
117 AliWarning("GRP/GRP/Data entry: missing value for the beam energy !");
118 AliError("\t ZDC does not reconstruct event 4 UNKNOWN beam type\n");
121 if(((beamType.CompareTo("pp"))==0) || ((beamType.CompareTo("p-p"))==0)
122 ||((beamType.CompareTo("PP"))==0) || ((beamType.CompareTo("P-P"))==0 ) fRecoMode=1;
123 else if((beamType.CompareTo("A-A")) == 0){
125 if(fIsCalibrationMB == kTRUE){
126 fRecoParam = new AliZDCRecoParamPbPb();
128 TH2F* hZDCvsZEM = new TH2F("hZDCvsZEM","hZDCvsZEM",100,0.,10.,100,0.,1000.);
129 hZDCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCvsZEM->SetYTitle("E_{ZDC} (TeV)");
130 fRecoParam->SetZDCvsZEM(hZDCvsZEM);
132 TH2F* hZDCCvsZEM = new TH2F("hZDCCvsZEM","hZDCCvsZEM",100,0.,10.,100,0.,500.);
133 hZDCCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCCvsZEM->SetYTitle("E_{ZDCC} (TeV)");
134 fRecoParam->SetZDCCvsZEM(hZDCCvsZEM);
136 TH2F* hZDCAvsZEM = new TH2F("hZDCAvsZEM","hZDCAvsZEM",100,0.,10.,100,0.,500.);
137 hZDCAvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCAvsZEM->SetYTitle("E_{ZDCA} (TeV)");
138 fRecoParam->SetZDCAvsZEM(hZDCAvsZEM);
140 AliInfo("\n ***** CALIBRATION_MB data -> building AliZDCRecoParamPbPb object *****");
144 fBeamEnergy = grpData->GetBeamEnergy();
145 if(fBeamEnergy==AliGRPObject::GetInvalidFloat()){
146 AliWarning("GRP/GRP/Data entry: missing value for the beam energy ! Using 0.");
150 if(fIsCalibrationMB==kFALSE)
151 printf("\n\n ***** ZDC reconstruction initialized for %s @ %1.0f GeV *****\n\n",beamType.Data(), fBeamEnergy);
154 AliError(" ATTENTION!!!!!! No beam type nor beam energy has been set!!!!!!\n");
159 //_____________________________________________________________________________
160 void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
162 // *** Local ZDC reconstruction for digits
163 // Works on the current event
165 // Retrieving calibration data
166 // Parameters for mean value pedestal subtraction
168 Float_t meanPed[2*kNch];
169 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
170 // Parameters pedestal subtraction through correlation with out-of-time signals
171 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
172 for(Int_t jj=0; jj<2*kNch; jj++){
173 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
174 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
179 AliZDCDigit* pdigit = &digit;
180 digitsTree->SetBranchAddress("ZDC", &pdigit);
181 //printf("\n\t # of digits in tree: %d\n",(Int_t) digitsTree->GetEntries());
184 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
185 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
186 for(Int_t i=0; i<10; i++){
187 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
188 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
191 Int_t digNentries = digitsTree->GetEntries();
192 Float_t ootDigi[kNch];
193 // -- Reading out-of-time signals (last kNch entries) for current event
195 for(Int_t iDigit=kNch; iDigit<digNentries; iDigit++){
196 ootDigi[iDigit] = digitsTree->GetEntry(iDigit);
200 for(Int_t iDigit=0; iDigit<(digNentries/2); iDigit++) {
201 digitsTree->GetEntry(iDigit);
202 if (!pdigit) continue;
204 Int_t det = digit.GetSector(0);
205 Int_t quad = digit.GetSector(1);
207 Float_t ped2SubHg=0., ped2SubLg=0.;
209 if(det==1) pedindex = quad;
210 else if(det==2) pedindex = quad+5;
211 else if(det==3) pedindex = quad+9;
212 else if(det==4) pedindex = quad+12;
213 else if(det==5) pedindex = quad+17;
215 else pedindex = (det-1)/3+22;
218 ped2SubHg = meanPed[pedindex];
219 ped2SubLg = meanPed[pedindex+kNch];
221 else if(fPedSubMode==1){
222 ped2SubHg = corrCoeff1[pedindex]*ootDigi[pedindex]+corrCoeff0[pedindex];
223 ped2SubLg = corrCoeff1[pedindex+kNch]*ootDigi[pedindex+kNch]+corrCoeff0[pedindex+kNch];
226 if(quad != 5){ // ZDC (not reference PTMs!)
227 if(det == 1){ // *** ZNC
228 tZN1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
229 tZN1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
230 if(tZN1Corr[quad]<0.) tZN1Corr[quad] = 0.;
231 if(tZN1Corr[quad+5]<0.) tZN1Corr[quad+5] = 0.;
233 //printf("\t pedindex %d tZN1Corr[%d] = %1.0f tZN1Corr[%d] = %1.0f",
234 // pedindex, quad, tZN1Corr[quad], quad+5, tZN1Corr[quad+5]);
236 else if(det == 2){ // *** ZP1
237 tZP1Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
238 tZP1Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
239 if(tZP1Corr[quad]<0.) tZP1Corr[quad] = 0.;
240 if(tZP1Corr[quad+5]<0.) tZP1Corr[quad+5] = 0.;
242 //printf("\t pedindex %d tZP1Corr[%d] = %1.0f tZP1Corr[%d] = %1.0f",
243 // pedindex, quad, tZP1Corr[quad], quad+5, tZP1Corr[quad+5]);
246 if(quad == 1){ // *** ZEM1
247 dZEM1Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
248 dZEM1Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
249 if(dZEM1Corr[0]<0.) dZEM1Corr[0] = 0.;
250 if(dZEM1Corr[1]<0.) dZEM1Corr[1] = 0.;
252 //printf("\t pedindex %d tZEM1Corr[%d] = %1.0f tZEM1Corr[%d] = %1.0f",
253 // pedindex, quad, tZEM1Corr[quad], quad+1, tZEM1Corr[quad+1]);
255 else if(quad == 2){ // *** ZEM2
256 dZEM2Corr[0] += (Float_t) (digit.GetADCValue(0)-ped2SubHg);
257 dZEM2Corr[1] += (Float_t) (digit.GetADCValue(1)-ped2SubLg);
258 if(dZEM2Corr[0]<0.) dZEM2Corr[0] = 0.;
259 if(dZEM2Corr[1]<0.) dZEM2Corr[1] = 0.;
261 //printf("\t pedindex %d tZEM2Corr[%d] = %1.0f tZEM2Corr[%d] = %1.0f",
262 // pedindex, quad, tZEM2Corr[quad], quad+1, tZEM2Corr[quad+1]);
265 else if(det == 4){ // *** ZN2
266 tZN2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
267 tZN2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
268 if(tZN2Corr[quad]<0.) tZN2Corr[quad] = 0.;
269 if(tZN2Corr[quad+5]<0.) tZN2Corr[quad+5] = 0.;
271 //printf("\t pedindex %d tZN2Corr[%d] = %1.0f tZN2Corr[%d] = %1.0f\n",
272 // pedindex, quad, tZN2Corr[quad], quad+5, tZN2Corr[quad+5]);
274 else if(det == 5){ // *** ZP2
275 tZP2Corr[quad] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
276 tZP2Corr[quad+5] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
277 if(tZP2Corr[quad]<0.) tZP2Corr[quad] = 0.;
278 if(tZP2Corr[quad+5]<0.) tZP2Corr[quad+5] = 0.;
280 //printf("\t pedindex %d tZP2Corr[%d] = %1.0f tZP2Corr[%d] = %1.0f\n",
281 // pedindex, quad, tZP2Corr[quad], quad+5, tZP2Corr[quad+5]);
284 else{ // Reference PMs
286 sPMRef1[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
287 sPMRef1[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
289 if(sPMRef1[0]<0.) sPMRef1[0] = 0.;
290 if(sPMRef2[1]<0.) sPMRef1[1] = 0.;
293 sPMRef2[0] = (Float_t) (digit.GetADCValue(0)-ped2SubHg);
294 sPMRef2[1] = (Float_t) (digit.GetADCValue(1)-ped2SubLg);
296 if(sPMRef2[0]<0.) sPMRef2[0] = 0.;
297 if(sPMRef2[1]<0.) sPMRef2[1] = 0.;
302 /*printf(" - AliZDCReconstructor -> digit #%d det %d quad %d pedHG %1.0f pedLG %1.0f\n",
303 iDigit, det, quad, ped2SubHg, ped2SubLg);
304 printf(" HGChain -> RawDig %d DigCorr %1.2f\n", digit.GetADCValue(0), digit.GetADCValue(0)-ped2SubHg);
305 printf(" LGChain -> RawDig %d DigCorr %1.2f\n", digit.GetADCValue(1), digit.GetADCValue(1)-ped2SubLg);
310 // If CALIBRATION_MB run build the RecoParam object
311 if(fIsCalibrationMB){
312 Float_t ZDCC=0., ZDCA=0., ZEM=0;
313 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
314 for(Int_t jkl=0; jkl<5; jkl++){
315 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
316 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
318 // Using energies in TeV in fRecoParam object
319 BuildRecoParam(ZDCC/1000., ZDCA/1000., ZEM/1000.);
321 Bool_t recFlag1 = kFALSE, recFlag2 = kFALSE, recFlag3 = kFALSE;
322 // reconstruct the event
324 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
325 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, recFlag1, recFlag2, recFlag3);
326 else if(fRecoMode==2)
327 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
328 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, recFlag1, recFlag2, recFlag3);
331 //_____________________________________________________________________________
332 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
334 // *** ZDC raw data reconstruction
335 // Works on the current event
337 // Retrieving calibration data
338 // Parameters for pedestal subtraction
340 Float_t meanPed[2*kNch];
341 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
342 // Parameters pedestal subtraction through correlation with out-of-time signals
343 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
344 for(Int_t jj=0; jj<2*kNch; jj++){
345 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
346 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
349 Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
350 Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
351 Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
352 Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
353 Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
354 Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
355 for(Int_t ich=0; ich<5; ich++){
356 adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
357 adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
358 adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
359 adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
361 adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
362 pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
366 // loop over raw data
367 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
368 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
369 for(Int_t i=0; i<10; i++){
370 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
371 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
374 //fNRun = (Int_t) rawReader->GetRunNumber();
375 Bool_t chOff=kFALSE, isUndflw=kFALSE, isOvflw=kFALSE;
378 AliZDCRawStream rawData(rawReader);
379 while(rawData.Next()){
380 if(rawData.IsCalibration() == kFALSE){ // Reading ADCs
381 //printf(" **** Reading ADC raw data **** \n");
383 if(rawData.GetNChannelsOn() < 48 ) chOff=kTRUE;
384 if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)) isUndflw=kTRUE;
385 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)) isOvflw=kTRUE;
387 if((rawData.IsADCDataWord()) && (isUndflw==kFALSE) && (isOvflw==kFALSE)){
389 Int_t adcMod = rawData.GetADCModule();
390 Int_t det = rawData.GetSector(0);
391 Int_t quad = rawData.GetSector(1);
392 Int_t gain = rawData.GetADCGain();
395 // Mean pedestal value subtraction -------------------------------------------------------
396 if(fPedSubMode == 0){
397 // Not interested in o.o.t. signals (ADC modules 2, 3)
398 if(adcMod == 2 || adcMod == 3) continue;
400 if(quad != 5){ // ZDCs (not reference PTMs)
403 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
404 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
408 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
409 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
414 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
415 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
418 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
419 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
424 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
425 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
429 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
430 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
433 else{ // reference PM
434 pedindex = (det-1)/3 + 22;
436 if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
437 else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
440 if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
441 else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
445 /*printf(" -> AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f\n",
446 det,quad,gain, pedindex, meanPed[pedindex]);
447 printf(" -> AliZDCReconstructor: RawADC %1.0f ADCCorr %1.0f\n",
448 rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);*/
450 }// mean pedestal subtraction
451 // Pedestal subtraction from correlation ------------------------------------------------
452 else if(fPedSubMode == 1){
454 if(adcMod==0 || adcMod==1){
455 if(quad != 5){ // signals from ZDCs
457 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
458 else adcZN1lg[quad] = rawData.GetADCValue();
461 if(gain==0) adcZP1[quad] = rawData.GetADCValue();
462 else adcZP1lg[quad] = rawData.GetADCValue();
465 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
466 else adcZEMlg[quad-1] = rawData.GetADCValue();
469 if(gain==0) adcZN2[quad] = rawData.GetADCValue();
470 else adcZN2lg[quad] = rawData.GetADCValue();
473 if(gain==0) adcZP2[quad] = rawData.GetADCValue();
474 else adcZP2lg[quad] = rawData.GetADCValue();
477 else{ // signals from reference PM
478 if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
479 else pmReflg[quad-1] = rawData.GetADCValue();
482 // Out-of-time pedestals
483 else if(adcMod==2 || adcMod==3){
484 if(quad != 5){ // signals from ZDCs
486 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
487 else adcZN1ootlg[quad] = rawData.GetADCValue();
490 if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
491 else adcZP1ootlg[quad] = rawData.GetADCValue();
494 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
495 else adcZEMootlg[quad-1] = rawData.GetADCValue();
498 if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
499 else adcZN2ootlg[quad] = rawData.GetADCValue();
502 if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
503 else adcZP2ootlg[quad] = rawData.GetADCValue();
506 else{ // signals from reference PM
507 if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
508 else pmRefootlg[quad-1] = rawData.GetADCValue();
511 } // pedestal subtraction from correlation
513 //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
514 // det,quad,gain, pedindex, meanPed[pedindex]);
516 }// Not raw data from calibration run!
522 for(Int_t t=0; t<5; t++){
523 tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
524 tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
526 tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
527 tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
529 tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
530 tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
532 tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
533 tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
534 // 0---------0 Ch. debug 0---------0
535 /* printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
536 printf("\tCorrCoeff0\tCorrCoeff1\n");
537 printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
538 printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
539 printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
540 printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
541 printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
542 printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
543 printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
544 printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
546 printf("ZN1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
547 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
548 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
549 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
551 printf("ZP1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
552 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
553 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
554 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
556 printf("ZN2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
557 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
558 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
559 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
561 printf("ZP2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
562 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
563 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
564 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
567 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
568 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
569 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
570 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
572 sPMRef1[0] = pmRef[0] - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
573 sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
574 sPMRef2[0] = pmRef[0] - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
575 sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
578 // If CALIBRATION_MB run build the RecoParam object
579 if(fIsCalibrationMB){
580 Float_t ZDCC=0., ZDCA=0., ZEM=0;
581 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
582 for(Int_t jkl=0; jkl<5; jkl++){
583 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
584 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
586 BuildRecoParam(ZDCC/100., ZDCA/100., ZEM/100.);
588 // reconstruct the event
590 if(fRecoMode==1) // p-p data
591 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
592 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, chOff, isUndflw, isOvflw);
593 else if(fRecoMode==2) // Pb-Pb data
594 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
595 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2, chOff, isUndflw, isOvflw);
599 //_____________________________________________________________________________
600 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* corrADCZN1,
601 Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
602 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
603 Bool_t channelsOff, Bool_t chUnderflow, Bool_t chOverflow) const
605 // ****************** Reconstruct one event ******************
607 // ---- Setting reco flags
610 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
612 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
613 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
614 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
615 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
616 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
617 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
619 if(channelsOff==kTRUE) rFlags[8] = 0x1;
620 if(chUnderflow == kTRUE) rFlags[9] = 0x1;
621 if(chOverflow==kTRUE) rFlags[10] = 0x1;
622 recoFlag = rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
623 rFlags[5] << 5 | rFlags[4] << 4 | rFlags[3] << 3 |
624 rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
626 // ****** Retrieving calibration data
627 // --- Equalization coefficients ---------------------------------------------
628 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
629 for(Int_t ji=0; ji<5; ji++){
630 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
631 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
632 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
633 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
635 // --- Energy calibration factors ------------------------------------
637 // **** Energy calibration coefficient set to 1
638 // **** (no trivial way to calibrate in p-p runs)
639 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
641 // ****** Equalization of detector responses
642 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
643 for(Int_t gi=0; gi<10; gi++){
645 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
646 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
647 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
648 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
651 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
652 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
653 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
654 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
658 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
659 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
660 for(Int_t gi=0; gi<5; gi++){
661 calibSumZN1[0] += equalTowZN1[gi];
662 calibSumZP1[0] += equalTowZP1[gi];
663 calibSumZN2[0] += equalTowZN2[gi];
664 calibSumZP2[0] += equalTowZP2[gi];
666 calibSumZN1[1] += equalTowZN1[gi+5];
667 calibSumZP1[1] += equalTowZP1[gi+5];
668 calibSumZN2[1] += equalTowZN2[gi+5];
669 calibSumZP2[1] += equalTowZP2[gi+5];
672 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
673 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
674 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
675 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
677 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
678 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
679 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
680 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
682 // ****** Energy calibration of detector responses
683 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
684 for(Int_t gi=0; gi<5; gi++){
686 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
687 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
688 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
689 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
691 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
692 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
693 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
694 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
697 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
698 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
699 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
700 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
701 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
702 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
704 // ****** No. of spectator and participants nucleons
705 // Variables calculated to comply with ESD structure
706 // *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
707 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
708 Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
709 Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
710 Double_t impPar=0., impPar1=0., impPar2=0.;
712 // create the output tree
713 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
714 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
715 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
716 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
717 nGenSpec, nGenSpecLeft, nGenSpecRight,
718 nPart, nPartTotLeft, nPartTotRight,
719 impPar, impPar1, impPar2,
722 const Int_t kBufferSize = 4000;
723 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
724 // write the output tree
725 clustersTree->Fill();
728 //_____________________________________________________________________________
729 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
730 Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
731 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
732 Bool_t channelsOff, Bool_t chUnderflow, Bool_t chOverflow) const
734 // ****************** Reconstruct one event ******************
736 // ---- Setting reco flags
739 for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
741 if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
742 if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
743 if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
744 if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
745 if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
746 if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
748 if(channelsOff==kTRUE) rFlags[8] = 0x1;
749 if(chUnderflow == kTRUE) rFlags[9] = 0x1;
750 if(chOverflow==kTRUE) rFlags[10] = 0x1;
751 recoFlag = rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
752 rFlags[5] << 5 | rFlags[4] << 4 | rFlags[3] << 3 |
753 rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
755 // ****** Retrieving calibration data
756 // --- Equalization coefficients ---------------------------------------------
757 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
758 for(Int_t ji=0; ji<5; ji++){
759 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
760 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
761 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
762 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
764 // --- Energy calibration factors ------------------------------------
765 Float_t valFromOCDB[6], calibEne[6];
766 for(Int_t ij=0; ij<6; ij++){
767 valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
769 if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
770 else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
772 else calibEne[ij] = valFromOCDB[ij];
775 // ****** Equalization of detector responses
776 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
777 for(Int_t gi=0; gi<10; gi++){
779 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
780 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
781 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
782 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
785 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
786 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
787 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
788 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
792 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
793 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
794 for(Int_t gi=0; gi<5; gi++){
795 calibSumZN1[0] += equalTowZN1[gi];
796 calibSumZP1[0] += equalTowZP1[gi];
797 calibSumZN2[0] += equalTowZN2[gi];
798 calibSumZP2[0] += equalTowZP2[gi];
800 calibSumZN1[1] += equalTowZN1[gi+5];
801 calibSumZP1[1] += equalTowZP1[gi+5];
802 calibSumZN2[1] += equalTowZN2[gi+5];
803 calibSumZP2[1] += equalTowZP2[gi+5];
806 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
807 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
808 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
809 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
811 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
812 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
813 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
814 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
816 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
817 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
818 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
819 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
820 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
821 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
823 // ****** Energy calibration of detector responses
824 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
825 for(Int_t gi=0; gi<5; gi++){
827 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
828 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
829 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
830 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
832 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
833 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
834 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
835 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
838 // ****** Number of detected spectator nucleons
839 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
841 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
842 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
843 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
844 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
846 else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
847 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
848 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
849 nDetSpecNRight, nDetSpecPRight);*/
851 Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
852 Int_t nPart=0, nPartA=0, nPartC=0;
853 Double_t b=0., bA=0., bC=0.;
855 if(fIsCalibrationMB == kFALSE){
856 // ****** Reconstruction parameters ------------------
858 //fRecoParam->Print("");
860 TH2F *hZDCvsZEM = fRecoParam->GethZDCvsZEM();
861 TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
862 TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
863 TH1D *hNpartDist = fRecoParam->GethNpartDist();
864 TH1D *hbDist = fRecoParam->GethbDist();
865 Float_t ClkCenter = fRecoParam->GetClkCenter();
867 Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
868 Double_t origin = xHighEdge*ClkCenter;
870 //printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
872 // ====> Summed ZDC info (sideA+side C)
873 TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
874 Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
875 Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
876 line->SetParameter(0, y/(x-origin));
877 line->SetParameter(1, -origin*y/(x-origin));
879 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
880 //printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
882 Double_t countPerc=0;
883 Double_t xBinCenter=0, yBinCenter=0;
884 for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
885 for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
886 xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
887 yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
889 if(line->GetParameter(0)>0){
890 if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
891 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
893 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
894 xBinCenter, yBinCenter, countPerc);*/
898 if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
899 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
901 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
902 xBinCenter, yBinCenter, countPerc);*/
908 Double_t xSecPerc = 0.;
909 if(hZDCvsZEM->GetEntries()!=0){
910 xSecPerc = countPerc/hZDCvsZEM->GetEntries();
913 AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
916 //printf(" xSecPerc %1.4f \n", xSecPerc);
919 TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
920 Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
921 lineC->SetParameter(0, yC/(x-origin));
922 lineC->SetParameter(1, -origin*yC/(x-origin));
924 //printf(" ***************** Side C \n");
925 //printf(" E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
927 Double_t countPercC=0;
928 Double_t xBinCenterC=0, yBinCenterC=0;
929 for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
930 for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
931 xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
932 yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
933 if(lineC->GetParameter(0)>0){
934 if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
935 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
939 if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
940 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
946 Double_t xSecPercC = 0.;
947 if(hZDCCvsZEM->GetEntries()!=0){
948 xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
951 AliWarning(" Histogram hZDCCvsZEM from OCDB has no entries!!!");
954 //printf(" xSecPercC %1.4f \n", xSecPercC);
957 TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
958 Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
959 lineA->SetParameter(0, yA/(x-origin));
960 lineA->SetParameter(1, -origin*yA/(x-origin));
963 //printf(" ***************** Side A \n");
964 //printf(" E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
966 Double_t countPercA=0;
967 Double_t xBinCenterA=0, yBinCenterA=0;
968 for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
969 for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
970 xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
971 yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
972 if(lineA->GetParameter(0)>0){
973 if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
974 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
978 if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
979 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
985 Double_t xSecPercA = 0.;
986 if(hZDCAvsZEM->GetEntries()!=0){
987 xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
990 AliWarning(" Histogram hZDCAvsZEM from OCDB has no entries!!!");
993 //printf(" xSecPercA %1.4f \n", xSecPercA);
995 // ****** Number of participants (from E_ZDC vs. E_ZEM correlation)
996 Int_t nPart=0, nPartC=0, nPartA=0;
997 Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
998 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
999 nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1000 if((1.-nPartFrac) < xSecPerc){
1001 nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1003 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
1004 //printf(" nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
1008 if(nPart<0) nPart=0;
1010 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1011 nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1012 if((1.-nPartFracC) < xSecPercC){
1013 nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1015 //printf(" ***************** Side C \n");
1016 //printf(" nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
1020 if(nPartC<0) nPartC=0;
1022 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
1023 nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
1024 if((1.-nPartFracA) < xSecPercA){
1025 nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
1027 //printf(" ***************** Side A \n");
1028 //printf(" nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
1032 if(nPartA<0) nPartA=0;
1034 // ****** Impact parameter (from E_ZDC vs. E_ZEM correlation)
1035 Float_t b=0, bC=0, bA=0;
1036 Double_t bFrac=0., bFracC=0., bFracA=0.;
1037 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1038 bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1039 if((1.-bFrac) < xSecPerc){
1040 b = hbDist->GetBinLowEdge(ibbin);
1045 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1046 bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1047 if((1.-bFracC) < xSecPercC){
1048 bC = hbDist->GetBinLowEdge(ibbin);
1053 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1054 bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1055 if((1.-bFracA) < xSecPercA){
1056 bA = hbDist->GetBinLowEdge(ibbin);
1061 // ****** Number of spectator nucleons
1062 Int_t nGenSpec=0, nGenSpecC=0, nGenSpecA=0;
1064 nGenSpec = 416 - nPart;
1065 nGenSpecC = 416 - nPartC;
1066 nGenSpecA = 416 - nPartA;
1067 if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1068 if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1069 if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1072 /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1073 " calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
1074 calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1075 printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n",
1076 nGenSpecLeft, nGenSpecRight);
1077 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1080 delete lineC; delete lineA;
1082 } // ONLY IF fIsCalibrationMB==kFALSE
1084 AliZDCReco* reco = new AliZDCReco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1085 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1086 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1087 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1088 nGenSpec, nGenSpecA, nGenSpecC,
1089 nPart, nPartA, nPartC, b, bA, bC,
1092 const Int_t kBufferSize = 4000;
1093 clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
1095 // write the output tree
1096 clustersTree->Fill();
1099 //_____________________________________________________________________________
1100 void AliZDCReconstructor::BuildRecoParam(Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1102 // Calculate RecoParam object for Pb-Pb data
1103 (fRecoParam->GethZDCvsZEM())->Fill(ZDCC+ZDCA, ZEM);
1104 (fRecoParam->GethZDCCvsZEM())->Fill(ZDCC, ZEM);
1105 (fRecoParam->GethZDCAvsZEM())->Fill(ZDCA, ZEM);
1109 //_____________________________________________________________________________
1110 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1112 // fill energies and number of participants to the ESD
1114 if(fIsCalibrationMB==kTRUE) WritePbPbRecoParamInOCDB();
1117 AliZDCReco* preco = &reco;
1118 clustersTree->SetBranchAddress("ZDC", &preco);
1120 clustersTree->GetEntry(0);
1122 AliESDZDC * esdzdc = esd->GetESDZDC();
1123 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1124 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1125 for(Int_t i=0; i<5; i++){
1126 tZN1Ene[i] = reco.GetZN1HREnTow(i);
1127 tZN2Ene[i] = reco.GetZN2HREnTow(i);
1128 tZP1Ene[i] = reco.GetZP1HREnTow(i);
1129 tZP2Ene[i] = reco.GetZP2HREnTow(i);
1131 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1132 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1133 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1134 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1137 esdzdc->SetZN1TowerEnergy(tZN1Ene);
1138 esdzdc->SetZN2TowerEnergy(tZN2Ene);
1139 esdzdc->SetZP1TowerEnergy(tZP1Ene);
1140 esdzdc->SetZP2TowerEnergy(tZP2Ene);
1142 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1143 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1144 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1145 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1147 Int_t nPart = reco.GetNParticipants();
1148 Int_t nPartA = reco.GetNPartSideA();
1149 Int_t nPartC = reco.GetNPartSideC();
1150 Double_t b = reco.GetImpParameter();
1151 Double_t bA = reco.GetImpParSideA();
1152 Double_t bC = reco.GetImpParSideC();
1153 UInt_t recoFlag = reco.GetRecoFlag();
1155 esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
1156 reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1157 nPart, nPartA, nPartC, b, bA, bC, recoFlag);
1161 //_____________________________________________________________________________
1162 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1164 // Setting the storage
1166 Bool_t deleteManager = kFALSE;
1168 AliCDBManager *manager = AliCDBManager::Instance();
1169 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1171 if(!defstorage || !(defstorage->Contains("ZDC"))){
1172 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1173 manager->SetDefaultStorage(uri);
1174 deleteManager = kTRUE;
1177 AliCDBStorage *storage = manager->GetDefaultStorage();
1180 AliCDBManager::Instance()->UnsetDefaultStorage();
1181 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1187 //_____________________________________________________________________________
1188 AliZDCPedestals* AliZDCReconstructor::GetPedestalData() const
1191 // Getting pedestal calibration object for ZDC set
1193 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1194 if(!entry) AliFatal("No calibration data loaded!");
1196 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
1197 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1202 //_____________________________________________________________________________
1203 AliZDCEnCalib* AliZDCReconstructor::GetEnergyCalibData() const
1206 // Getting energy and equalization calibration object for ZDC set
1208 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnergyCalib");
1209 if(!entry) AliFatal("No calibration data loaded!");
1211 AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1212 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1217 //_____________________________________________________________________________
1218 AliZDCTowerCalib* AliZDCReconstructor::GetTowerCalibData() const
1221 // Getting energy and equalization calibration object for ZDC set
1223 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowerCalib");
1224 if(!entry) AliFatal("No calibration data loaded!");
1226 AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1227 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1232 //_____________________________________________________________________________
1233 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1236 // Getting reconstruction parameters from OCDB
1238 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1239 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1241 AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1242 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1248 //_____________________________________________________________________________
1249 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1252 // Getting reconstruction parameters from OCDB
1254 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1255 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1257 AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1258 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1264 //_____________________________________________________________________________
1265 void AliZDCReconstructor::WritePbPbRecoParamInOCDB() const
1268 // Writing Pb-Pb reconstruction parameters from OCDB
1270 AliCDBManager *man = AliCDBManager::Instance();
1271 AliCDBMetaData *md= new AliCDBMetaData();
1272 md->SetResponsible("Chiara Oppedisano");
1273 md->SetComment("ZDC Pb-Pb reconstruction parameters");
1274 md->SetObjectClassName("AliZDCRecoParamPbPb");
1275 AliCDBId id("ZDC/Calib/RecoParamPbPb",fNRun,AliCDBRunRange::Infinity());
1276 man->Put(fRecoParam, id, md);