Little changes needed for running with new AliESDEvent
[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx
CommitLineData
8309c1ab 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// class for ZDC reconstruction //
21// //
22///////////////////////////////////////////////////////////////////////////////
23
24
25#include <TF1.h>
26
27#include "AliRunLoader.h"
28#include "AliRawReader.h"
af885e0f 29#include "AliESDEvent.h"
8309c1ab 30#include "AliZDCDigit.h"
31#include "AliZDCRawStream.h"
32#include "AliZDCReco.h"
33#include "AliZDCReconstructor.h"
48642b09 34#include "AliZDCCalibData.h"
8309c1ab 35
36
37ClassImp(AliZDCReconstructor)
38
39
40//_____________________________________________________________________________
cc2abffd 41AliZDCReconstructor:: AliZDCReconstructor() :
42
43 fZNCen (new TF1("fZNCen",
44 "(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.)),
45 fZNPer (new TF1("fZNPer",
46 "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.)),
47 fZPCen (new TF1("fZPCen",
48 "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.)),
49 fZPPer (new TF1("fZPPer",
50 "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.)),
51 fZDCCen (new TF1("fZDCCen",
52 "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.)),
53 fZDCPer (new TF1("fZDCPer",
54 "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.)),
55 fbCen (new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.)),
56 fbPer (new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.)),
57 fZEMn (new TF1("fZEMn","126.2-0.05399*x+0.000005679*x*x",0.,4000.)),
58 fZEMp (new TF1("fZEMp","82.49-0.03611*x+0.00000385*x*x",0.,4000.)),
59 fZEMsp (new TF1("fZEMsp","208.7-0.09006*x+0.000009526*x*x",0.,4000.)),
60 fZEMb (new TF1("fZEMb",
61 "16.06-0.01633*x+1.44e-5*x*x-6.778e-9*x*x*x+1.438e-12*x*x*x*x-1.112e-16*x*x*x*x*x",0.,4000.)),
62 fCalibData(GetCalibData())
8309c1ab 63
8309c1ab 64{
cc2abffd 65 // **** Default constructor
8309c1ab 66
8309c1ab 67}
68
8309c1ab 69
70//_____________________________________________________________________________
71AliZDCReconstructor::~AliZDCReconstructor()
72{
73// destructor
74
75 delete fZNCen;
76 delete fZNPer;
77 delete fZPCen;
78 delete fZPPer;
79 delete fZDCCen;
80 delete fZDCPer;
81 delete fbCen;
82 delete fbPer;
83 delete fZEMn;
84 delete fZEMp;
85 delete fZEMsp;
86 delete fZEMb;
87}
88
89
90//_____________________________________________________________________________
70f04f6d 91void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
8309c1ab 92{
48642b09 93 // *** Local ZDC reconstruction for digits
70f04f6d 94 // Works on the current event
78d18275 95
48642b09 96 Float_t meanPed[47];
78d18275 97 for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
8309c1ab 98
70f04f6d 99 // get digits
8309c1ab 100 AliZDCDigit digit;
101 AliZDCDigit* pdigit = &digit;
70f04f6d 102 digitsTree->SetBranchAddress("ZDC", &pdigit);
103
104 // loop over digits
105 Float_t zn1corr=0, zp1corr=0, zn2corr=0, zp2corr=0, zemcorr=0;
106 for (Int_t iDigit = 0; iDigit < digitsTree->GetEntries(); iDigit++) {
107 digitsTree->GetEntry(iDigit);
108 if (!pdigit) continue;
109
110 if(digit.GetSector(0) == 1)
111 zn1corr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)]); // high gain ZN1 ADCs
112 else if(digit.GetSector(0) == 2)
113 zp1corr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+10]); // high gain ZP1 ADCs
114 else if(digit.GetSector(0) == 3){
115 if(digit.GetSector(1)==1)
116 zemcorr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+20]); // high gain ZEM1 ADCs
117 else if(digit.GetSector(1)==2)
118 zemcorr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+22]); // high gain ZEM2 ADCs
8309c1ab 119 }
70f04f6d 120 else if(digit.GetSector(0) == 4)
121 zn2corr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+24]); // high gain ZN2 ADCs
122 else if(digit.GetSector(0) == 5)
123 zp2corr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+34]); // high gain ZP2 ADCs
8309c1ab 124 }
70f04f6d 125 if(zn1corr<0) zn1corr=0;
126 if(zp1corr<0) zp1corr=0;
127 if(zn2corr<0) zn2corr=0;
128 if(zp2corr<0) zp2corr=0;
129 if(zemcorr<0) zemcorr=0;
130
131 // reconstruct the event
132 //printf("\n \t ZDCReco from digits-> ZN = %.0f, ZP = %.0f, ZEM = %.0f\n",zncorr,zpcorr,zemcorr);
133 ReconstructEvent(clustersTree, zn1corr, zp1corr, zemcorr, zn2corr, zp2corr);
8309c1ab 134
8309c1ab 135}
136
137//_____________________________________________________________________________
70f04f6d 138void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
8309c1ab 139{
70f04f6d 140 // *** ZDC raw data reconstruction
141 // Works on the current event
48642b09 142
48642b09 143 Float_t meanPed[47];
78d18275 144 for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
8309c1ab 145
70f04f6d 146 rawReader->Reset();
147
148 // loop over raw data digits
149 Float_t zn1corr=0, zp1corr=0, zn2corr=0, zp2corr=0,zemcorr=0;
150 AliZDCRawStream digit(rawReader);
151 while (digit.Next()) {
152 if(digit.IsADCDataWord()){
153 if(digit.GetADCGain() == 0){
154 if(digit.GetSector(0) == 1)
155 zn1corr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)]); // high gain ZN1 ADCs;
156 else if(digit.GetSector(0) == 2)
157 zp1corr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+10]); // high gain ZP1 ADCs;
158 else if(digit.GetSector(0) == 3){
159 if(digit.GetSector(1)==1)
160 zemcorr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+20]); // high gain ZEM1 ADCs
161 else if(digit.GetSector(1)==2)
162 zemcorr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+22]); // high gain ZEM2 ADCs
48642b09 163 }
70f04f6d 164 else if(digit.GetSector(0) == 4)
165 zn2corr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+24]); // high gain ZN2 ADCs;
166 else if(digit.GetSector(0) == 5)
167 zp2corr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+34]); // high gain ZP2 ADCs;
8309c1ab 168 }
169 }
8309c1ab 170 }
70f04f6d 171 if(zn1corr<0) zn1corr=0;
172 if(zp1corr<0) zp1corr=0;
173 if(zn2corr<0) zn2corr=0;
174 if(zp2corr<0) zp2corr=0;
175 if(zemcorr<0) zemcorr=0;
176
177 // reconstruct the event
178 //printf("\n\t ZDCReco from raw-> ZN = %.0f, ZP = %.0f, ZEM = %.0f\n",zncorr,zpcorr,zemcorr);
179 ReconstructEvent(clustersTree, zn1corr, zp1corr, zemcorr, zn2corr, zp2corr);
8309c1ab 180
8309c1ab 181}
182
183//_____________________________________________________________________________
70f04f6d 184void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree, Float_t zn1corr,
980685f2 185 Float_t zp1corr, Float_t zemcorr, Float_t zn2corr, Float_t zp2corr) const
8309c1ab 186{
48642b09 187 // ***** Reconstruct one event
8309c1ab 188
189 // --- ADCchannel -> photoelectrons
190 // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
191 // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
980685f2 192 Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
193 zn1phe = zn1corr/convFactor;
194 zp1phe = zp1corr/convFactor;
8309c1ab 195 zemphe = zemcorr/convFactor;
980685f2 196 zn2phe = zn2corr/convFactor;
197 zp2phe = zp2corr/convFactor;
48642b09 198 //if AliDebug(1,Form("\n znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
8309c1ab 199
200 // --- Energy calibration
201 // Conversion factors for hadronic ZDCs goes from phe yield to TRUE
202 // incident energy (conversion from GeV to TeV is included); while for EM
203 // calos conversion is from light yield to detected energy calculated by
204 // GEANT NB -> ZN and ZP conversion factors are constant since incident
205 // spectators have all the same energy, ZEM energy is obtained through a
206 // fit over the whole range of incident particle energies
207 // (obtained with full HIJING simulations)
980685f2 208 Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
209 Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
210 zn1energy = zn1phe/zn1phexTeV;
211 zp1energy = zp1phe/zp1phexTeV;
212 zdc1energy = zn1energy+zp1energy;
213 zn2energy = zn2phe/zn2phexTeV;
214 zp2energy = zp2phe/zp2phexTeV;
215 zdc2energy = zn2energy+zp2energy;
8309c1ab 216 zemenergy = -4.81+0.3238*zemphe;
217 if(zemenergy<0) zemenergy=0;
48642b09 218 // if AliDebug(1,Form(" znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
219 // "\n zemenergy = %f TeV\n", znenergy, zpenergy,
220 // zdcenergy, zemenergy);
48642b09 221 // if(zdcenergy==0)
222 // if AliDebug(1,Form("\n\n ### ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
223 // " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy);
8309c1ab 224
980685f2 225 // --- Number of detected spectator nucleons
226 // *** N.B. -> It works only in Pb-Pb
227 Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
228 nDetSpecNLeft = (Int_t) (zn1energy/2.760);
229 nDetSpecPLeft = (Int_t) (zp1energy/2.760);
230 nDetSpecNRight = (Int_t) (zn2energy/2.760);
231 nDetSpecPRight = (Int_t) (zp2energy/2.760);
232
233 // --- Number of generated spectator nucleons (from HIJING parameterization)
234 // *** N.B. -> Only one side!!!
235 Int_t nGenSpecN=0, nGenSpecP=0, nGenSpec=0;
8309c1ab 236 Double_t impPar=0;
237 // Cut value for Ezem (GeV)
980685f2 238 // ### Results from production -> 0<b<18 fm (Apr 2002)
8309c1ab 239 Float_t eZEMCut = 420.;
240 Float_t deltaEZEMSup = 690.;
241 Float_t deltaEZEMInf = 270.;
242 if(zemenergy > (eZEMCut+deltaEZEMSup)){
980685f2 243 nGenSpecN = (Int_t) (fZNCen->Eval(zn1energy));
244 nGenSpecP = (Int_t) (fZPCen->Eval(zp1energy));
245 nGenSpec = (Int_t) (fZDCCen->Eval(zdc1energy));
246 impPar = fbCen->Eval(zdc1energy);
8309c1ab 247 }
248 else if(zemenergy < (eZEMCut-deltaEZEMInf)){
980685f2 249 nGenSpecN = (Int_t) (fZNPer->Eval(zn1energy));
250 nGenSpecP = (Int_t) (fZPPer->Eval(zp1energy));
251 nGenSpec = (Int_t) (fZDCPer->Eval(zdc1energy));
252 impPar = fbPer->Eval(zdc1energy);
8309c1ab 253 }
254 else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
255 nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy));
256 nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy));
257 nGenSpec = (Int_t)(fZEMsp->Eval(zemenergy));
258 impPar = fZEMb->Eval(zemenergy);
8309c1ab 259 }
980685f2 260 // ### Results from production -> 0<b<18 fm (Apr 2002)
261 if(zn1energy>162.) nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy));
262 if(zp1energy>59.75) nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy));
263 if(zdc1energy>221.5) nGenSpec = (Int_t)(fZEMsp->Eval(zemenergy));
264 if(zdc1energy>220.) impPar = fZEMb->Eval(zemenergy);
8309c1ab 265
266 if(nGenSpecN>125) nGenSpecN=125;
267 else if(nGenSpecN<0) nGenSpecN=0;
268 if(nGenSpecP>82) nGenSpecP=82;
269 else if(nGenSpecP<0) nGenSpecP=0;
270 if(nGenSpec>207) nGenSpec=207;
271 else if(nGenSpec<0) nGenSpec=0;
8309c1ab 272
980685f2 273 // --- Number of generated participants (from HIJING parameterization)
8309c1ab 274 Int_t nPart, nPartTot;
275 nPart = 207-nGenSpecN-nGenSpecP;
276 nPartTot = 207-nGenSpec;
78d18275 277 //printf("\t ZDCeventReco-> ZNEn = %.0f GeV, ZPEn = %.0f GeV, ZEMEn = %.0f GeV\n",
278 // znenergy, zpenergy, zemenergy);
8309c1ab 279
70f04f6d 280 // get the output tree and store the ZDC reco object there
980685f2 281 AliZDCReco reco(zn1energy, zp1energy, zdc1energy, zemenergy,
282 zn2energy, zp2energy, zdc2energy,
283 nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
284 nGenSpecN, nGenSpecP, nGenSpec,nPartTot, impPar);
8309c1ab 285 AliZDCReco* preco = &reco;
286 const Int_t kBufferSize = 4000;
70f04f6d 287 clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
8309c1ab 288
289 // write the output tree
70f04f6d 290 clustersTree->Fill();
8309c1ab 291}
292
293//_____________________________________________________________________________
70f04f6d 294void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
8309c1ab 295{
70f04f6d 296 // fill energies and number of participants to the ESD
8309c1ab 297
8309c1ab 298 AliZDCReco reco;
299 AliZDCReco* preco = &reco;
70f04f6d 300 clustersTree->SetBranchAddress("ZDC", &preco);
8309c1ab 301
70f04f6d 302 clustersTree->GetEntry(0);
980685f2 303 esd->SetZDC(reco.GetZN1energy(), reco.GetZP1energy(), reco.GetZEMenergy(),
304 reco.GetZN2energy(), reco.GetZP2energy(), reco.GetNPart());
8309c1ab 305}
48642b09 306
307//_____________________________________________________________________________
78d18275 308AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
48642b09 309{
cc2abffd 310 // Setting the storage
48642b09 311
78d18275 312 Bool_t deleteManager = kFALSE;
48642b09 313
78d18275 314 AliCDBManager *manager = AliCDBManager::Instance();
315 AliCDBStorage *defstorage = manager->GetDefaultStorage();
48642b09 316
78d18275 317 if(!defstorage || !(defstorage->Contains("ZDC"))){
318 AliWarning("No default storage set or default storage doesn't contain ZDC!");
319 manager->SetDefaultStorage(uri);
320 deleteManager = kTRUE;
321 }
322
323 AliCDBStorage *storage = manager->GetDefaultStorage();
324
325 if(deleteManager){
326 AliCDBManager::Instance()->UnsetDefaultStorage();
327 defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
328 }
329
330 return storage;
331}
48642b09 332
78d18275 333//_____________________________________________________________________________
4fda3ba1 334AliZDCCalibData* AliZDCReconstructor::GetCalibData() const
78d18275 335{
48642b09 336
4fda3ba1 337 // Getting calibration object for ZDC set
338
339 AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Data");
457a440d 340 if(!entry) AliFatal("No calibration data loaded!");
4fda3ba1 341
457a440d 342 AliZDCCalibData *calibdata = dynamic_cast<AliZDCCalibData*> (entry->GetObject());
343 if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
48642b09 344
78d18275 345 return calibdata;
48642b09 346}