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(GetPedData()),
57 fEnCalibData(GetEnCalibData()),
58 fTowCalibData(GetTowCalibData()),
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;
114 fRecoParam = new AliZDCRecoParamPbPb();
116 TH2F* hZDCvsZEM = new TH2F("hZDCvsZEM","hZDCvsZEM",100,0.,10.,100,0.,1000.);
117 hZDCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCvsZEM->SetYTitle("E_{ZDC} (TeV)");
118 fRecoParam->SetZDCvsZEM(hZDCvsZEM);
120 TH2F* hZDCCvsZEM = new TH2F("hZDCCvsZEM","hZDCCvsZEM",100,0.,10.,100,0.,500.);
121 hZDCCvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCCvsZEM->SetYTitle("E_{ZDCC} (TeV)");
122 fRecoParam->SetZDCCvsZEM(hZDCCvsZEM);
124 TH2F* hZDCAvsZEM = new TH2F("hZDCAvsZEM","hZDCAvsZEM",100,0.,10.,100,0.,500.);
125 hZDCAvsZEM->SetXTitle("E_{ZEM} (TeV)"); hZDCAvsZEM->SetYTitle("E_{ZDCA} (TeV)");
126 fRecoParam->SetZDCAvsZEM(hZDCAvsZEM);
129 TString beamType = grpData->GetBeamType();
130 if(beamType==AliGRPObject::GetInvalidString()){
131 AliWarning("GRP/GRP/Data entry: missing value for the beam energy !");
132 AliError("\t ZDC does not reconstruct event 4 UNKNOWN beam type\n");
135 if((beamType.CompareTo("p-p")) == 0){
137 fRecoParam = (AliZDCRecoParampp*) GetppRecoParamFromOCDB();
139 else if((beamType.CompareTo("A-A")) == 0){
141 if(fIsCalibrationMB == kFALSE)
142 fRecoParam = (AliZDCRecoParamPbPb*) GetPbPbRecoParamFromOCDB();
145 fBeamEnergy = grpData->GetBeamEnergy();
146 if(fBeamEnergy==AliGRPObject::GetInvalidFloat()){
147 AliWarning("GRP/GRP/Data entry: missing value for the beam energy ! Using 0.");
151 if(fIsCalibrationMB==kTRUE){
152 AliInfo("\n ***** CALIBRATION_MB data -> building AliZDCRecoParamPbPb object *****");
155 AliInfo(Form("\n\n ***** ZDC reconstruction initialized for %s @ %1.3f GeV *****\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);
314 // If CALIBRATION_MB run build the RecoParam object
315 if(fIsCalibrationMB){
316 Float_t ZDCC=0., ZDCA=0., ZEM=0;
317 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
318 for(Int_t jkl=0; jkl<5; jkl++){
319 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
320 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
322 // Using energies in TeV in fRecoParam object
323 BuildRecoParam(fRecoParam->GethZDCvsZEM(), fRecoParam->GethZDCCvsZEM(),
324 fRecoParam->GethZDCAvsZEM(), ZDCC/1000., ZDCA/1000., ZEM/1000.);
326 // reconstruct the event
328 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
329 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2);
330 else if(fRecoMode==2)
331 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
332 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2);
335 //_____________________________________________________________________________
336 void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree)
338 // *** ZDC raw data reconstruction
339 // Works on the current event
341 // Retrieving calibration data
342 // Parameters for pedestal subtraction
344 Float_t meanPed[2*kNch];
345 for(Int_t jj=0; jj<2*kNch; jj++) meanPed[jj] = fPedData->GetMeanPed(jj);
346 // Parameters pedestal subtraction through correlation with out-of-time signals
347 Float_t corrCoeff0[2*kNch], corrCoeff1[2*kNch];
348 for(Int_t jj=0; jj<2*kNch; jj++){
349 corrCoeff0[jj] = fPedData->GetPedCorrCoeff0(jj);
350 corrCoeff1[jj] = fPedData->GetPedCorrCoeff1(jj);
353 Int_t adcZN1[5], adcZN1oot[5], adcZN1lg[5], adcZN1ootlg[5];
354 Int_t adcZP1[5], adcZP1oot[5], adcZP1lg[5], adcZP1ootlg[5];
355 Int_t adcZN2[5], adcZN2oot[5], adcZN2lg[5], adcZN2ootlg[5];
356 Int_t adcZP2[5], adcZP2oot[5], adcZP2lg[5], adcZP2ootlg[5];
357 Int_t adcZEM[2], adcZEMoot[2], adcZEMlg[2], adcZEMootlg[2];
358 Int_t pmRef[2], pmRefoot[2], pmReflg[2], pmRefootlg[2];
359 for(Int_t ich=0; ich<5; ich++){
360 adcZN1[ich] = adcZN1oot[ich] = adcZN1lg[ich] = adcZN1ootlg[ich] = 0;
361 adcZP1[ich] = adcZP1oot[ich] = adcZP1lg[ich] = adcZP1ootlg[ich] = 0;
362 adcZN2[ich] = adcZN2oot[ich] = adcZN2lg[ich] = adcZN2ootlg[ich] = 0;
363 adcZP2[ich] = adcZP2oot[ich] = adcZP2lg[ich] = adcZP2ootlg[ich] = 0;
365 adcZEM[ich] = adcZEMoot[ich] = adcZEMlg[ich] = adcZEMootlg[ich] = 0;
366 pmRef[ich] = pmRefoot[ich] = pmReflg[ich] = pmRefootlg[ich] = 0;
370 // loop over raw data
371 Float_t tZN1Corr[10], tZP1Corr[10], tZN2Corr[10], tZP2Corr[10];
372 Float_t dZEM1Corr[2], dZEM2Corr[2], sPMRef1[2], sPMRef2[2];
373 for(Int_t i=0; i<10; i++){
374 tZN1Corr[i] = tZP1Corr[i] = tZN2Corr[i] = tZP2Corr[i] = 0.;
375 if(i<2) dZEM1Corr[i] = dZEM2Corr[i] = sPMRef1[i] = sPMRef2[i] = 0.;
379 fNRun = (Int_t) rawReader->GetRunNumber();
380 AliZDCRawStream rawData(rawReader);
381 while(rawData.Next()){
383 Bool_t ch2process = kTRUE;
385 // Setting reco flags (part I)
386 if((rawData.IsADCDataWord()) && (rawData.IsUnderflow() == kTRUE)){
390 if((rawData.IsADCDataWord()) && (rawData.IsOverflow() == kTRUE)){
391 fRecoFlag = 0x1 << 7;
394 if(rawData.GetNChannelsOn() < 48 ) fRecoFlag = 0x1 << 6;
396 if((rawData.IsADCDataWord()) && (ch2process == kTRUE)){
398 Int_t adcMod = rawData.GetADCModule();
399 Int_t det = rawData.GetSector(0);
400 Int_t quad = rawData.GetSector(1);
401 Int_t gain = rawData.GetADCGain();
404 // Mean pedestal value subtraction -------------------------------------------------------
405 if(fPedSubMode == 0){
406 // Not interested in o.o.t. signals (ADC modules 2, 3)
407 if(adcMod == 2 || adcMod == 3) return;
409 if(quad != 5){ // ZDCs (not reference PTMs)
412 if(gain == 0) tZN1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
413 else tZN1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
417 if(gain == 0) tZP1Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
418 else tZP1Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
423 if(gain == 0) dZEM1Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
424 else dZEM1Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
427 if(gain == 0) dZEM2Corr[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
428 else dZEM2Corr[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
433 if(gain == 0) tZN2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
434 else tZN2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
438 if(gain == 0) tZP2Corr[quad] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
439 else tZP2Corr[quad+5] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex+kNch]);
442 else{ // reference PM
443 pedindex = (det-1)/3 + 22;
445 if(gain==0) sPMRef1[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
446 else sPMRef1[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
449 if(gain==0) sPMRef2[0] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
450 else sPMRef2[1] += (Float_t) (rawData.GetADCValue()-meanPed[pedindex]);
454 /*printf(" -> AliZDCReconstructor: det %d quad %d res %d -> Pedestal[%d] %1.0f\n",
455 det,quad,gain, pedindex, meanPed[pedindex]);
456 printf(" -> AliZDCReconstructor: RawADC %1.0f ADCCorr %1.0f\n",
457 rawData.GetADCValue(), rawData.GetADCValue()-meanPed[pedindex]);*/
459 }// mean pedestal subtraction
460 // Pedestal subtraction from correlation ------------------------------------------------
461 else if(fPedSubMode == 1){
463 if(adcMod==0 || adcMod==1){
464 if(quad != 5){ // signals from ZDCs
466 if(gain==0) adcZN1[quad] = rawData.GetADCValue();
467 else adcZN1lg[quad] = rawData.GetADCValue();
470 if(gain==0) adcZP1[quad] = rawData.GetADCValue();
471 else adcZP1lg[quad] = rawData.GetADCValue();
474 if(gain==0) adcZEM[quad-1] = rawData.GetADCValue();
475 else adcZEMlg[quad-1] = rawData.GetADCValue();
478 if(gain==0) adcZN2[quad] = rawData.GetADCValue();
479 else adcZN2lg[quad] = rawData.GetADCValue();
482 if(gain==0) adcZP2[quad] = rawData.GetADCValue();
483 else adcZP2lg[quad] = rawData.GetADCValue();
486 else{ // signals from reference PM
487 if(gain==0) pmRef[quad-1] = rawData.GetADCValue();
488 else pmReflg[quad-1] = rawData.GetADCValue();
491 // Out-of-time pedestals
492 else if(adcMod==2 || adcMod==3){
493 if(quad != 5){ // signals from ZDCs
495 if(gain==0) adcZN1oot[quad] = rawData.GetADCValue();
496 else adcZN1ootlg[quad] = rawData.GetADCValue();
499 if(gain==0) adcZP1oot[quad] = rawData.GetADCValue();
500 else adcZP1ootlg[quad] = rawData.GetADCValue();
503 if(gain==0) adcZEMoot[quad-1] = rawData.GetADCValue();
504 else adcZEMootlg[quad-1] = rawData.GetADCValue();
507 if(gain==0) adcZN2oot[quad] = rawData.GetADCValue();
508 else adcZN2ootlg[quad] = rawData.GetADCValue();
511 if(gain==0) adcZP2oot[quad] = rawData.GetADCValue();
512 else adcZP2ootlg[quad] = rawData.GetADCValue();
515 else{ // signals from reference PM
516 if(gain==0) pmRefoot[quad-1] = rawData.GetADCValue();
517 else pmRefootlg[quad-1] = rawData.GetADCValue();
520 } // pedestal subtraction from correlation
522 //printf("\t AliZDCReconstructor - det %d quad %d res %d -> Ped[%d] = %1.0f\n",
523 // det,quad,gain, pedindex, meanPed[pedindex]);
528 for(Int_t t=0; t<5; t++){
529 tZN1Corr[t] = adcZN1[t] - (corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]);
530 tZN1Corr[t+5] = adcZN1lg[t] - (corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]);
532 tZP1Corr[t] = adcZP1[t] - (corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]);
533 tZP1Corr[t+5] = adcZP1lg[t] - (corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]);
535 tZN2Corr[t] = adcZN2[t] - (corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]);
536 tZN2Corr[t+5] = adcZN2lg[t] - (corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]);
538 tZP2Corr[t] = adcZP2[t] - (corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]);
539 tZP2Corr[t+5] = adcZP2lg[t] - (corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]);
540 // 0---------0 Ch. debug 0---------0
541 /* printf("\n\n ---------- Debug of pedestal subtraction from correlation ----------\n");
542 printf("\tCorrCoeff0\tCorrCoeff1\n");
543 printf(" ZN1 %d\t%1.0f\t%1.0f\n",t,corrCoeff0[t],corrCoeff1[t]);
544 printf(" ZN1lg %d\t%1.0f\t%1.0f\n",t+kNch,corrCoeff0[t+kNch],corrCoeff1[t+kNch]);
545 printf(" ZP1 %d\t%1.0f\t%1.0f\n",t+5,corrCoeff0[t+5],corrCoeff1[t+5]);
546 printf(" ZP1lg %d\t%1.0f\t%1.0f\n",t+5+kNch,corrCoeff0[t+5+kNch],corrCoeff1[t+5+kNch]);
547 printf(" ZN2 %d\t%1.0f\t%1.0f\n",t+12,corrCoeff0[t+12],corrCoeff1[t+12]);
548 printf(" ZN2lg %d\t%1.0f\t%1.0f\n",t+12+kNch,corrCoeff0[t+12+kNch],corrCoeff1[t+12+kNch]);
549 printf(" ZP2 %d\t%1.0f\t%1.0f\n",t+17,corrCoeff0[t+17],corrCoeff1[t+17]);
550 printf(" ZP2lg %d\t%1.0f\t%1.0f\n",t+17+kNch,corrCoeff0[t+17+kNch],corrCoeff1[t+17+kNch]);
552 printf("ZN1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
553 adcZN1[t],(corrCoeff1[t]*adcZN1oot[t]+corrCoeff0[t]),tZN1Corr[t]);
554 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
555 adcZN1lg[t],(corrCoeff1[t+kNch]*adcZN1ootlg[t]+corrCoeff0[t+kNch]),tZN1Corr[t+5]);
557 printf("ZP1 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
558 adcZP1[t],(corrCoeff1[t+5]*adcZP1oot[t]+corrCoeff0[t+5]),tZP1Corr[t]);
559 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
560 adcZP1lg[t],(corrCoeff1[t+5+kNch]*adcZP1ootlg[t]+corrCoeff0[t+5+kNch]),tZP1Corr[t+5]);
562 printf("ZN2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
563 adcZN2[t],(corrCoeff1[t+12]*adcZN2oot[t]+corrCoeff0[t+12]),tZN2Corr[t]);
564 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
565 adcZN2lg[t],(corrCoeff1[t+12+kNch]*adcZN2ootlg[t]+corrCoeff0[t+12+kNch]),tZN2Corr[t+5]);
567 printf("ZP2 -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
568 adcZP2[t],(corrCoeff1[t+17]*adcZP2oot[t]+corrCoeff0[t+17]),tZP2Corr[t]);
569 printf(" lg -> rawADC %d\tpedestal%1.2f\tcorrADC%1.2f\n",
570 adcZP2lg[t],(corrCoeff1[t+17+kNch]*adcZP2ootlg[t]+corrCoeff0[t+17+kNch]),tZP2Corr[t+5]);
573 dZEM1Corr[0] = adcZEM[0] - (corrCoeff1[9]*adcZEMoot[0]+corrCoeff0[9]);
574 dZEM1Corr[1] = adcZEMlg[0] - (corrCoeff1[9+kNch]*adcZEMootlg[0]+corrCoeff0[9+kNch]);
575 dZEM2Corr[0] = adcZEM[1] - (corrCoeff1[10]*adcZEMoot[1]+corrCoeff0[10]);
576 dZEM2Corr[1] = adcZEMlg[1] - (corrCoeff1[10+kNch]*adcZEMootlg[1]+corrCoeff0[10+kNch]);
578 sPMRef1[0] = pmRef[0] - (corrCoeff1[22]*pmRefoot[0]+corrCoeff0[22]);
579 sPMRef1[1] = pmReflg[0] - (corrCoeff1[22+kNch]*pmRefootlg[0]+corrCoeff0[22+kNch]);
580 sPMRef2[0] = pmRef[0] - (corrCoeff1[23]*pmRefoot[1]+corrCoeff0[23]);
581 sPMRef2[1] = pmReflg[0] - (corrCoeff1[23+kNch]*pmRefootlg[1]+corrCoeff0[23+kNch]);
583 // Setting reco flags (part II)
584 Float_t sumZNAhg=0, sumZPAhg=0, sumZNChg=0, sumZPChg=0;
585 for(Int_t jj=0; jj<5; jj++){
586 sumZNAhg += tZN2Corr[jj];
587 sumZPAhg += tZP2Corr[jj];
588 sumZNChg += tZN1Corr[jj];
589 sumZPChg += tZP1Corr[jj];
591 if(sumZNAhg>0.) fRecoFlag = 0x1;
592 if(sumZPAhg>0.) fRecoFlag = 0x1 << 1;
593 if(dZEM1Corr[0]>0.) fRecoFlag = 0x1 << 2;
594 if(dZEM2Corr[0]>0.) fRecoFlag = 0x1 << 3;
595 if(sumZNChg>0.) fRecoFlag = 0x1 << 4;
596 if(sumZPChg>0.) fRecoFlag = 0x1 << 5;
598 // If CALIBRATION_MB run build the RecoParam object
599 if(fIsCalibrationMB){
600 Float_t ZDCC=0., ZDCA=0., ZEM=0;
601 ZEM += dZEM1Corr[0] + dZEM2Corr[0];
602 for(Int_t jkl=0; jkl<5; jkl++){
603 ZDCC += tZN1Corr[jkl] + tZP1Corr[jkl];
604 ZDCA += tZN2Corr[jkl] + tZP2Corr[jkl];
606 BuildRecoParam(fRecoParam->GethZDCvsZEM(), fRecoParam->GethZDCCvsZEM(),
607 fRecoParam->GethZDCAvsZEM(), ZDCC/100., ZDCA/100., ZEM/100.);
609 // reconstruct the event
611 if(fRecoMode==1) // p-p data
612 ReconstructEventpp(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
613 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2);
614 else if(fRecoMode==2) // Pb-Pb data
615 ReconstructEventPbPb(clustersTree, tZN1Corr, tZP1Corr, tZN2Corr, tZP2Corr,
616 dZEM1Corr, dZEM2Corr, sPMRef1, sPMRef2);
620 //_____________________________________________________________________________
621 void AliZDCReconstructor::ReconstructEventpp(TTree *clustersTree, Float_t* corrADCZN1,
622 Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
623 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2) const
625 // ****************** Reconstruct one event ******************
627 // ****** Retrieving calibration data
628 // --- Equalization coefficients ---------------------------------------------
629 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
630 for(Int_t ji=0; ji<5; ji++){
631 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
632 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
633 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
634 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
636 // --- Energy calibration factors ------------------------------------
638 // **** Energy calibration coefficient set to 1
639 // **** (no trivial way to calibrate in p-p runs)
640 for(Int_t ij=0; ij<4; ij++) calibEne[ij] = fEnCalibData->GetEnCalib(ij);
642 // ****** Equalization of detector responses
643 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
644 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 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
652 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
653 for(Int_t gi=0; gi<5; gi++){
654 calibSumZN1[0] += equalTowZN1[gi];
655 calibSumZP1[0] += equalTowZP1[gi];
656 calibSumZN2[0] += equalTowZN2[gi];
657 calibSumZP2[0] += equalTowZP2[gi];
659 calibSumZN1[1] += equalTowZN1[gi+5];
660 calibSumZP1[1] += equalTowZP1[gi+5];
661 calibSumZN2[1] += equalTowZN2[gi+5];
662 calibSumZP2[1] += equalTowZP2[gi+5];
665 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
666 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
667 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
668 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
670 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
671 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
672 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
673 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
675 // ****** Energy calibration of detector responses
676 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
677 for(Int_t gi=0; gi<5; gi++){
679 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
680 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
681 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
682 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
684 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
685 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
686 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
687 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
690 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
691 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
692 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
693 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
694 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
695 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
697 // ****** No. of spectator and participants nucleons
698 // Variables calculated to comply with ESD structure
699 // *** N.B. -> They have a meaning only in Pb-Pb!!!!!!!!!!!!
700 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
701 Int_t nGenSpec=0, nGenSpecLeft=0, nGenSpecRight=0;
702 Int_t nPart=0, nPartTotLeft=0, nPartTotRight=0;
703 Double_t impPar=0., impPar1=0., impPar2=0.;
705 // create the output tree
706 AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
707 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
708 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
709 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
710 nGenSpec, nGenSpecLeft, nGenSpecRight,
711 nPart, nPartTotLeft, nPartTotRight,
712 impPar, impPar1, impPar2);
714 AliZDCReco* preco = &reco;
715 const Int_t kBufferSize = 4000;
716 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
718 // write the output tree
719 clustersTree->Fill();
722 //_____________________________________________________________________________
723 void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
724 Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
725 Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2) const
727 // ****************** Reconstruct one event ******************
729 // ****** Retrieving calibration data
730 // --- Equalization coefficients ---------------------------------------------
731 Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
732 for(Int_t ji=0; ji<5; ji++){
733 equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
734 equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
735 equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
736 equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
738 // --- Energy calibration factors ------------------------------------
739 Float_t valFromOCDB[6], calibEne[6];
740 for(Int_t ij=0; ij<6; ij++){
741 valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
743 if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
744 else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
746 else calibEne[ij] = valFromOCDB[ij];
749 // ****** Equalization of detector responses
750 Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
751 for(Int_t gi=0; gi<10; gi++){
752 equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
753 equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
754 equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
755 equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
758 // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
759 Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
760 for(Int_t gi=0; gi<5; gi++){
761 calibSumZN1[0] += equalTowZN1[gi];
762 calibSumZP1[0] += equalTowZP1[gi];
763 calibSumZN2[0] += equalTowZN2[gi];
764 calibSumZP2[0] += equalTowZP2[gi];
766 calibSumZN1[1] += equalTowZN1[gi+5];
767 calibSumZP1[1] += equalTowZP1[gi+5];
768 calibSumZN2[1] += equalTowZN2[gi+5];
769 calibSumZP2[1] += equalTowZP2[gi+5];
772 calibSumZN1[0] = calibSumZN1[0]*calibEne[0]/8.;
773 calibSumZP1[0] = calibSumZP1[0]*calibEne[1]/8.;
774 calibSumZN2[0] = calibSumZN2[0]*calibEne[2]/8.;
775 calibSumZP2[0] = calibSumZP2[0]*calibEne[3]/8.;
777 calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
778 calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
779 calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
780 calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
782 Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
783 calibZEM1[0] = corrADCZEM1[0]*calibEne[5]/8.;
784 calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
785 calibZEM2[0] = corrADCZEM2[0]*calibEne[5]/8.;
786 calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
787 for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
789 // ****** Energy calibration of detector responses
790 Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
791 for(Int_t gi=0; gi<5; gi++){
793 calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]/8.;
794 calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]/8.;
795 calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]/8.;
796 calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]/8.;
798 calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
799 calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
800 calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
801 calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
804 // ****** Number of detected spectator nucleons
805 Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
807 nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
808 nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
809 nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
810 nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
812 else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
813 /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
814 " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
815 nDetSpecNRight, nDetSpecPRight);*/
817 if(fIsCalibrationMB == kFALSE){
818 // ****** Reconstruction parameters ------------------
820 //fRecoParam->Print("");
822 TH2F *hZDCvsZEM = fRecoParam->GethZDCvsZEM();
823 TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
824 TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
825 TH1D *hNpartDist = fRecoParam->GethNpartDist();
826 TH1D *hbDist = fRecoParam->GethbDist();
827 Float_t ClkCenter = fRecoParam->GetClkCenter();
829 Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
830 Double_t origin = xHighEdge*ClkCenter;
832 printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
834 // ====> Summed ZDC info (sideA+side C)
835 TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
836 Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
837 Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
838 line->SetParameter(0, y/(x-origin));
839 line->SetParameter(1, -origin*y/(x-origin));
841 printf(" ***************** Summed ZDC info (sideA+side C) \n");
842 printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
844 Double_t countPerc=0;
845 Double_t xBinCenter=0, yBinCenter=0;
846 for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
847 for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
848 xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
849 yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
851 if(line->GetParameter(0)>0){
852 if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
853 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
855 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
856 xBinCenter, yBinCenter, countPerc);*/
860 if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
861 countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
863 /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
864 xBinCenter, yBinCenter, countPerc);*/
870 Double_t xSecPerc = 0.;
871 if(hZDCvsZEM->GetEntries()!=0){
872 xSecPerc = countPerc/hZDCvsZEM->GetEntries();
875 AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
878 //printf(" xSecPerc %1.4f \n", xSecPerc);
881 TF1 *lineC = new TF1("lineC","[0]*x+[1]",0.,xHighEdge);
882 Float_t yC = (calibSumZN1[0]+calibSumZP1[0])/1000.;
883 lineC->SetParameter(0, yC/(x-origin));
884 lineC->SetParameter(1, -origin*yC/(x-origin));
886 //printf(" ***************** Side C \n");
887 //printf(" E_{ZEM} %1.4f, E_{ZDCC} %1.2f, TF1: %1.2f*x + %1.2f ", x, yC,yC/(x-origin),-origin*yC/(x-origin));
889 Double_t countPercC=0;
890 Double_t xBinCenterC=0, yBinCenterC=0;
891 for(Int_t nbinx=1; nbinx<=hZDCCvsZEM->GetNbinsX(); nbinx++){
892 for(Int_t nbiny=1; nbiny<=hZDCCvsZEM->GetNbinsY(); nbiny++){
893 xBinCenterC = hZDCCvsZEM->GetXaxis()->GetBinCenter(nbinx);
894 yBinCenterC = hZDCCvsZEM->GetYaxis()->GetBinCenter(nbiny);
895 if(lineC->GetParameter(0)>0){
896 if(yBinCenterC < (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
897 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
901 if(yBinCenterC > (lineC->GetParameter(0)*xBinCenterC + lineC->GetParameter(1))){
902 countPercC += hZDCCvsZEM->GetBinContent(nbinx,nbiny);
908 Double_t xSecPercC = 0.;
909 if(hZDCCvsZEM->GetEntries()!=0){
910 xSecPercC = countPercC/hZDCCvsZEM->GetEntries();
913 AliWarning(" Histogram hZDCCvsZEM from OCDB has no entries!!!");
916 //printf(" xSecPercC %1.4f \n", xSecPercC);
919 TF1 *lineA = new TF1("lineA","[0]*x+[1]",0.,xHighEdge);
920 Float_t yA = (calibSumZN2[0]+calibSumZP2[0])/1000.;
921 lineA->SetParameter(0, yA/(x-origin));
922 lineA->SetParameter(1, -origin*yA/(x-origin));
925 //printf(" ***************** Side A \n");
926 //printf(" E_{ZEM} %1.4f, E_{ZDCA} %1.2f, TF1: %1.2f*x + %1.2f ", x, yA,yA/(x-origin),-origin*yA/(x-origin));
928 Double_t countPercA=0;
929 Double_t xBinCenterA=0, yBinCenterA=0;
930 for(Int_t nbinx=1; nbinx<=hZDCAvsZEM->GetNbinsX(); nbinx++){
931 for(Int_t nbiny=1; nbiny<=hZDCAvsZEM->GetNbinsY(); nbiny++){
932 xBinCenterA = hZDCAvsZEM->GetXaxis()->GetBinCenter(nbinx);
933 yBinCenterA = hZDCAvsZEM->GetYaxis()->GetBinCenter(nbiny);
934 if(lineA->GetParameter(0)>0){
935 if(yBinCenterA < (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
936 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
940 if(yBinCenterA > (lineA->GetParameter(0)*xBinCenterA + lineA->GetParameter(1))){
941 countPercA += hZDCAvsZEM->GetBinContent(nbinx,nbiny);
947 Double_t xSecPercA = 0.;
948 if(hZDCAvsZEM->GetEntries()!=0){
949 xSecPercA = countPercA/hZDCAvsZEM->GetEntries();
952 AliWarning(" Histogram hZDCAvsZEM from OCDB has no entries!!!");
955 //printf(" xSecPercA %1.4f \n", xSecPercA);
957 // ****** Number of participants (from E_ZDC vs. E_ZEM correlation)
958 Int_t nPart=0, nPartC=0, nPartA=0;
959 Double_t nPartFrac=0., nPartFracC=0., nPartFracA=0.;
960 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
961 nPartFrac += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
962 if((1.-nPartFrac) < xSecPerc){
963 nPart = (Int_t) hNpartDist->GetBinLowEdge(npbin);
965 //printf(" ***************** Summed ZDC info (sideA+side C) \n");
966 //printf(" nPartFrac %1.4f, nPart %d\n", nPartFrac, nPart);
972 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
973 nPartFracC += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
974 if((1.-nPartFracC) < xSecPercC){
975 nPartC = (Int_t) hNpartDist->GetBinLowEdge(npbin);
977 //printf(" ***************** Side C \n");
978 //printf(" nPartFracC %1.4f, nPartC %d\n", nPartFracC, nPartC);
982 if(nPartC<0) nPartC=0;
984 for(Int_t npbin=1; npbin<hNpartDist->GetNbinsX(); npbin++){
985 nPartFracA += (hNpartDist->GetBinContent(npbin))/(hNpartDist->GetEntries());
986 if((1.-nPartFracA) < xSecPercA){
987 nPartA = (Int_t) hNpartDist->GetBinLowEdge(npbin);
989 //printf(" ***************** Side A \n");
990 //printf(" nPartFracA %1.4f, nPartA %d\n\n", nPartFracA, nPartA);
994 if(nPartA<0) nPartA=0;
996 // ****** Impact parameter (from E_ZDC vs. E_ZEM correlation)
997 Float_t b=0, bC=0, bA=0;
998 Double_t bFrac=0., bFracC=0., bFracA=0.;
999 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1000 bFrac += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1001 if((1.-bFrac) < xSecPerc){
1002 b = hbDist->GetBinLowEdge(ibbin);
1007 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1008 bFracC += (hbDist->GetBinContent(ibbin))/(hbDist->GetEntries());
1009 if((1.-bFracC) < xSecPercC){
1010 bC = hbDist->GetBinLowEdge(ibbin);
1015 for(Int_t ibbin=1; ibbin<hbDist->GetNbinsX(); ibbin++){
1016 bFracA += (hbDist->GetBinContent(ibbin))/(hNpartDist->GetEntries());
1017 if((1.-bFracA) < xSecPercA){
1018 bA = hbDist->GetBinLowEdge(ibbin);
1023 // ****** Number of spectator nucleons
1024 Int_t nGenSpec=0, nGenSpecC=0, nGenSpecA=0;
1026 nGenSpec = 416 - nPart;
1027 nGenSpecC = 416 - nPartC;
1028 nGenSpecA = 416 - nPartA;
1029 if(nGenSpec>416) nGenSpec=416; if(nGenSpec<0) nGenSpec=0;
1030 if(nGenSpecC>416) nGenSpecC=416; if(nGenSpecC<0) nGenSpecC=0;
1031 if(nGenSpecA>416) nGenSpecA=416; if(nGenSpecA<0) nGenSpecA=0;
1034 /*printf("\n\t AliZDCReconstructor -> calibSumZN1[0] %1.0f, calibSumZP1[0] %1.0f,"
1035 " calibSumZN2[0] %1.0f, calibSumZP2[0] %1.0f, corrADCZEMHG %1.0f\n",
1036 calibSumZN1[0],calibSumZP1[0],calibSumZN2[0],calibSumZP2[0],corrADCZEMHG);
1037 printf("\t AliZDCReconstructor -> nGenSpecLeft %d nGenSpecRight %d\n",
1038 nGenSpecLeft, nGenSpecRight);
1039 printf("\t AliZDCReconstructor -> NpartL %d, NpartR %d, b %1.2f fm\n\n",nPartTotLeft, nPartTotRight, impPar);
1042 // create the output tree
1043 AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1044 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1045 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1046 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1047 nGenSpec, nGenSpecA, nGenSpecC,
1048 nPart, nPartA, nPartC, b, bA, bC);
1050 AliZDCReco* preco = &reco;
1051 const Int_t kBufferSize = 4000;
1052 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
1053 // write the output tree
1054 clustersTree->Fill();
1056 delete lineC; delete lineA;
1057 } // ONLY IF fIsCalibrationMB==kFALSE
1059 // create the output tree
1060 AliZDCReco reco(calibSumZN1, calibSumZP1, calibSumZN2, calibSumZP2,
1061 calibTowZN1, calibTowZP1, calibTowZN2, calibTowZP2,
1062 calibZEM1, calibZEM2, sPMRef1, sPMRef2,
1063 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
1065 0, 0, 0, 0., 0., 0.);
1067 AliZDCReco* preco = &reco;
1068 const Int_t kBufferSize = 4000;
1069 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
1070 // write the output tree
1071 clustersTree->Fill();
1075 //_____________________________________________________________________________
1076 void AliZDCReconstructor::BuildRecoParam(TH2F* hCorr, TH2F* hCorrC, TH2F* hCorrA,
1077 Float_t ZDCC, Float_t ZDCA, Float_t ZEM) const
1079 // Calculate RecoParam object for Pb-Pb data
1080 hCorr->Fill(ZDCC+ZDCA, ZEM);
1081 hCorrC->Fill(ZDCC, ZEM);
1082 hCorrA->Fill(ZDCA, ZEM);
1086 Float_t clkCenter;*/
1090 //_____________________________________________________________________________
1091 void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
1093 // fill energies and number of participants to the ESD
1095 if(fIsCalibrationMB==kTRUE) WritePbPbRecoParamInOCDB();
1098 AliZDCReco* preco = &reco;
1099 clustersTree->SetBranchAddress("ZDC", &preco);
1101 clustersTree->GetEntry(0);
1103 AliESDZDC * esdzdc = esd->GetESDZDC();
1104 Float_t tZN1Ene[5], tZN2Ene[5], tZP1Ene[5], tZP2Ene[5];
1105 Float_t tZN1EneLR[5], tZN2EneLR[5], tZP1EneLR[5], tZP2EneLR[5];
1106 for(Int_t i=0; i<5; i++){
1107 tZN1Ene[i] = reco.GetZN1HREnTow(i);
1108 tZN2Ene[i] = reco.GetZN2HREnTow(i);
1109 tZP1Ene[i] = reco.GetZP1HREnTow(i);
1110 tZP2Ene[i] = reco.GetZP2HREnTow(i);
1112 tZN1EneLR[i] = reco.GetZN1LREnTow(i);
1113 tZN2EneLR[i] = reco.GetZN2LREnTow(i);
1114 tZP1EneLR[i] = reco.GetZP1LREnTow(i);
1115 tZP2EneLR[i] = reco.GetZP2LREnTow(i);
1118 esdzdc->SetZN1TowerEnergy(tZN1Ene);
1119 esdzdc->SetZN2TowerEnergy(tZN2Ene);
1120 esdzdc->SetZP1TowerEnergy(tZP1Ene);
1121 esdzdc->SetZP2TowerEnergy(tZP2Ene);
1123 esdzdc->SetZN1TowerEnergyLR(tZN1EneLR);
1124 esdzdc->SetZN2TowerEnergyLR(tZN2EneLR);
1125 esdzdc->SetZP1TowerEnergyLR(tZP1EneLR);
1126 esdzdc->SetZP2TowerEnergyLR(tZP2EneLR);
1128 Int_t nPart = reco.GetNParticipants();
1129 Int_t nPartA = reco.GetNPartSideA();
1130 Int_t nPartC = reco.GetNPartSideC();
1131 Double_t b = reco.GetImpParameter();
1132 Double_t bA = reco.GetImpParSideA();
1133 Double_t bC = reco.GetImpParSideC();
1135 esd->SetZDC(reco.GetZN1HREnergy(), reco.GetZP1HREnergy(), reco.GetZEM1HRsignal(),
1136 reco.GetZEM2HRsignal(), reco.GetZN2HREnergy(), reco.GetZP2HREnergy(),
1137 nPart, nPartA, nPartC, b, bA, bC, fRecoFlag);
1141 //_____________________________________________________________________________
1142 AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
1144 // Setting the storage
1146 Bool_t deleteManager = kFALSE;
1148 AliCDBManager *manager = AliCDBManager::Instance();
1149 AliCDBStorage *defstorage = manager->GetDefaultStorage();
1151 if(!defstorage || !(defstorage->Contains("ZDC"))){
1152 AliWarning("No default storage set or default storage doesn't contain ZDC!");
1153 manager->SetDefaultStorage(uri);
1154 deleteManager = kTRUE;
1157 AliCDBStorage *storage = manager->GetDefaultStorage();
1160 AliCDBManager::Instance()->UnsetDefaultStorage();
1161 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
1167 //_____________________________________________________________________________
1168 AliZDCPedestals* AliZDCReconstructor::GetPedData() const
1171 // Getting pedestal calibration object for ZDC set
1173 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
1174 if(!entry) AliFatal("No calibration data loaded!");
1176 AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*> (entry->GetObject());
1177 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1182 //_____________________________________________________________________________
1183 AliZDCEnCalib* AliZDCReconstructor::GetEnCalibData() const
1186 // Getting energy and equalization calibration object for ZDC set
1188 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/EnCalib");
1189 if(!entry) AliFatal("No calibration data loaded!");
1191 AliZDCEnCalib *calibdata = dynamic_cast<AliZDCEnCalib*> (entry->GetObject());
1192 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1197 //_____________________________________________________________________________
1198 AliZDCTowerCalib* AliZDCReconstructor::GetTowCalibData() const
1201 // Getting energy and equalization calibration object for ZDC set
1203 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/TowCalib");
1204 if(!entry) AliFatal("No calibration data loaded!");
1206 AliZDCTowerCalib *calibdata = dynamic_cast<AliZDCTowerCalib*> (entry->GetObject());
1207 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
1212 //_____________________________________________________________________________
1213 AliZDCRecoParampp* AliZDCReconstructor::GetppRecoParamFromOCDB() const
1216 // Getting reconstruction parameters from OCDB
1218 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParampp");
1219 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1221 AliZDCRecoParampp *param = dynamic_cast<AliZDCRecoParampp*> (entry->GetObject());
1222 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1228 //_____________________________________________________________________________
1229 AliZDCRecoParamPbPb* AliZDCReconstructor::GetPbPbRecoParamFromOCDB() const
1232 // Getting reconstruction parameters from OCDB
1234 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecoParamPbPb");
1235 if(!entry) AliFatal("No RecoParam data found in OCDB!");
1237 AliZDCRecoParamPbPb *param = dynamic_cast<AliZDCRecoParamPbPb*> (entry->GetObject());
1238 if(!param) AliFatal("No RecoParam object in OCDB entry!");
1244 //_____________________________________________________________________________
1245 void AliZDCReconstructor::WritePbPbRecoParamInOCDB() const
1248 // Writing Pb-Pb reconstruction parameters from OCDB
1250 AliCDBManager *man = AliCDBManager::Instance();
1251 AliCDBMetaData *md= new AliCDBMetaData();
1252 md->SetResponsible("Chiara Oppedisano");
1253 md->SetComment("ZDC Pb-Pb reconstruction parameters");
1254 md->SetObjectClassName("AliZDCRecoParamPbPb");
1255 AliCDBId id("ZDC/Calib/RecoParamPbPb",fNRun,AliCDBRunRange::Infinity());
1256 man->Put(fRecoParam, id, md);