]>
Commit | Line | Data |
---|---|---|
c335b3c5 | 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 | **************************************************************************/ | |
e7386040 | 15 | /* $Id$ */ |
c335b3c5 | 16 | |
12f6dd57 | 17 | // -- |
18 | // -- | |
19 | // Implementation for TTree output in PHOS DA | |
20 | // for calibrating energy by pi0 and MIP. | |
21 | // -- | |
22 | // -- Author: Hisayuki Torii (Hiroshima Univ.) | |
23 | // -- | |
24 | ||
25 | #include <stdio.h> | |
26 | #include <iostream> | |
27 | #include <math.h> | |
c335b3c5 | 28 | #include "AliLog.h" |
12f6dd57 | 29 | #include "AliPHOSDApi0mip.h" |
30 | ClassImp(AliPHOSDApi0mip) | |
31 | //---------------------------------------------------------------- | |
32 | AliPHOSDApi0mip::AliPHOSDApi0mip(int module,int iterid,char* fopt) : | |
33 | TNamed(), fCreateTree(false), fCreateHist(false), fMod(0), fIterId(0), | |
34 | fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0), | |
35 | fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0), | |
36 | fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) { | |
37 | // Constructor | |
38 | ||
39 | char hname[1024], htitle[1024]; | |
40 | sprintf(hname,"AliPHOSDApi0mip_mod%d_iter%d",module,iterid); | |
41 | SetName(hname); | |
42 | sprintf(htitle,"PHOS MIP/pi0 Calibration DA for Module:%d Iteration:%d",module,iterid); | |
43 | SetTitle(htitle); | |
44 | ||
45 | fMod = module; | |
46 | fIterId = iterid; | |
47 | fEvent = 0; | |
48 | fEventClustered = false; | |
49 | fTime = 0; | |
50 | fCreateTree = false; | |
51 | fCreateHist = false; | |
52 | ||
53 | char fname[1024]; | |
54 | sprintf(fname,"AliPHOSDApi0mip_mod%d.root",module); | |
55 | fTFile = TFile::Open(fname,fopt); | |
56 | } | |
57 | //---------------------------------------------------------------- | |
58 | AliPHOSDApi0mip::AliPHOSDApi0mip(const AliPHOSDApi0mip& da): | |
59 | TNamed(da), fCreateTree(false), fCreateHist(false), fMod(da.fMod), fIterId(da.fIterId), | |
60 | fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0), | |
61 | fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0), | |
62 | fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) { | |
63 | // Copy Constructor | |
64 | ||
65 | char fname[1024], hname[1024], htitle[1024];; | |
66 | sprintf(fname,"%s.root",GetName()); | |
67 | fTFile = TFile::Open(fname,"RECREATE"); | |
68 | fEvent = 0; | |
69 | fEventClustered = false; | |
70 | fTime = 0; | |
71 | ||
72 | if( da.fCreateHist ){ | |
73 | fCreateHist = true; | |
74 | fH1Time = (TH1I*) da.fH1Time->Clone(); | |
75 | fH1DigitNum = (TH1F*) da.fH1DigitNum->Clone(); | |
76 | fH1ClusterNum = (TH1F*) da.fH1ClusterNum->Clone(); | |
77 | fH2EneDigitId = (TH2F*) da.fH2EneDigitId->Clone(); | |
78 | fH2MipDigitId = (TH2F*) da.fH2MipDigitId->Clone(); | |
79 | fH2Pi0DigitId = (TH2F*) da.fH2Pi0DigitId->Clone(); | |
80 | fH3Pi0AsymPt = (TH3F*) da.fH3Pi0AsymPt->Clone(); | |
81 | } else { | |
82 | fCreateHist = false; | |
83 | fH1Time = 0; | |
84 | fH1DigitNum = 0; | |
85 | fH1ClusterNum = 0; | |
86 | fH2EneDigitId = 0; | |
87 | fH2MipDigitId = 0; | |
88 | fH2Pi0DigitId = 0; | |
89 | fH3Pi0AsymPt = 0; | |
90 | } | |
91 | if( da.fCreateTree ){ | |
92 | // Create new ttree instead of copy. | |
93 | fCreateTree = true; | |
94 | sprintf(hname,"tevt_mod%d_iter%d",fMod,fIterId); | |
95 | sprintf(htitle,"Calibration for Module:%d Iteration:%d",fMod,fIterId); | |
96 | fTTree = new TTree(hname,htitle); | |
97 | fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent); | |
98 | } else { | |
99 | fCreateTree = false; | |
100 | fTTree = 0; | |
101 | } | |
102 | } | |
103 | //------------------------------------------------------------------- | |
c335b3c5 | 104 | AliPHOSDApi0mip& AliPHOSDApi0mip::operator= (const AliPHOSDApi0mip&) |
105 | { | |
106 | // Operator= is not implemented yet | |
107 | AliFatal("Operator= is not implemented"); | |
108 | return *this; | |
109 | } | |
110 | //------------------------------------------------------------------- | |
12f6dd57 | 111 | AliPHOSDApi0mip::~AliPHOSDApi0mip(){ |
112 | // Destructor | |
113 | ||
114 | //std::cout<<" DEBUGDEBUG :: AliPHOSDApi0mip:: deconstructor "<<std::endl; | |
115 | if( fCreateHist ){ | |
116 | fH1Time->Write(0,TObject::kOverwrite); | |
117 | fH1DigitNum->Write(0,TObject::kOverwrite); | |
118 | fH1ClusterNum->Write(0,TObject::kOverwrite); | |
119 | fH2EneDigitId->Write(0,TObject::kOverwrite); | |
120 | fH2MipDigitId->Write(0,TObject::kOverwrite); | |
121 | fH2Pi0DigitId->Write(0,TObject::kOverwrite); | |
122 | fH3Pi0AsymPt->Write(0,TObject::kOverwrite); | |
123 | delete fH1Time; | |
124 | delete fH1DigitNum; | |
125 | delete fH1ClusterNum; | |
126 | delete fH2EneDigitId; | |
127 | delete fH2MipDigitId; | |
128 | delete fH2Pi0DigitId; | |
129 | delete fH3Pi0AsymPt; | |
130 | } | |
131 | if( fCreateTree ){ | |
132 | fTTree->Write(); | |
133 | } | |
134 | if( fEvent ) delete fEvent; | |
135 | fTFile->Close(); | |
136 | delete fTFile; | |
137 | } | |
138 | //------------------------------------------------------------------- | |
139 | void AliPHOSDApi0mip::FillDigit(float adc,int row,int col){ | |
140 | // Put new digit information in event | |
141 | // | |
142 | ||
143 | if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent; | |
144 | float energy = adc * 0.005; // 5MeV/ch is fixed. | |
145 | if( energy > 0.020 ){ // Threshold at 20MeV. | |
146 | //std::cout<<" AliPHOSDApi0mip::AppendDigit() DEBUG:: (row,col)=("<<row<<","<<col<<") adc="<<adc<<" energy="<<energy<<std::endl; | |
147 | fEvent->Fill(energy,row,col); | |
148 | } | |
149 | } | |
150 | //------------------------------------------------------------------- | |
151 | void AliPHOSDApi0mip::NewEvent(){ | |
152 | // Initiator must be called in the beginning of event | |
153 | // | |
154 | ||
155 | if(fEvent) fEvent->Reset(); | |
156 | fEventClustered = false; | |
157 | } | |
158 | //---------------------------------------------------------------- | |
159 | bool AliPHOSDApi0mip::CreateTree(){ | |
160 | if(! fTFile ) { | |
161 | return false; | |
162 | } | |
163 | fTFile->cd(); | |
164 | char hname[1024], htitle[1024]; | |
165 | sprintf(hname,"trevt_mod%d_iter%d",fMod,fIterId); | |
166 | sprintf(htitle,"Calibration for Module:%d Iteration:%d",fMod,fIterId); | |
167 | fTTree = new TTree(hname,htitle); | |
168 | fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent); | |
169 | fCreateTree = true; | |
170 | return true; | |
171 | } | |
172 | //---------------------------------------------------------------- | |
173 | bool AliPHOSDApi0mip::CreateHist(){ | |
174 | // Create histogram routine called by the constructor only. | |
175 | ||
176 | char hname[1024], htitle[1024]; | |
177 | ||
178 | sprintf(hname,"h1_time_mod%d_iter%d",fMod,fIterId); | |
179 | sprintf(htitle,"Time : Mod:%d Iter:%d",fMod,fIterId); | |
180 | fH1Time = (TH1I*) gDirectory->Get(hname); | |
181 | if( fH1Time>0 ){ | |
182 | std::cout<<" AliPHOSDApi0mip:Warning!! Output object already exist : "<<fH1Time->GetName()<<std::endl; | |
183 | return false; | |
184 | } | |
185 | fH1Time = new TH1I(hname,htitle,2,0,2); | |
186 | ||
187 | sprintf(hname,"h1_digitnum_mod%d_iter%d",fMod,fIterId); | |
188 | sprintf(htitle,"Number of Digits : Mod:%d Iter:%d",fMod,fIterId); | |
189 | fH1DigitNum = new TH1F(hname,htitle,100,0,100); | |
190 | ||
191 | sprintf(hname,"h1_clusternum_mod%d_iter%d",fMod,fIterId); | |
192 | sprintf(htitle,"Number of Clusters : Mod:%d Iter:%d",fMod,fIterId); | |
193 | fH1ClusterNum = new TH1F(hname,htitle,100,0,100); | |
194 | ||
195 | sprintf(hname,"h2_pi0digitid_mod%d_iter%d",fMod,fIterId); | |
196 | sprintf(htitle,"PHOS pi0 mass vs Digit Id : Mod:%d Iter:%d",fMod,fIterId); | |
197 | fH2Pi0DigitId = new TH2F(hname,htitle,17920,0,17920,120,0,0.3); | |
198 | ||
199 | sprintf(hname,"h2_mipdigitid_mod%d_iter%d",fMod,fIterId); | |
200 | sprintf(htitle,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId); | |
201 | fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025); | |
202 | //fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,50,0.5,1.5); | |
203 | ||
204 | sprintf(hname,"h2_enedigitid_mod%d_iter%d",fMod,fIterId); | |
205 | sprintf(htitle,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId); | |
206 | //fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025); | |
207 | fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,50,0,5.0); | |
208 | ||
209 | sprintf(hname,"h3_pi0asympt_mod%d_iter%d",fMod,fIterId); | |
210 | sprintf(htitle,"PHOS pi0 mass vs pT vs Asym : Mod:%d Iter:%d",fMod,fIterId); | |
211 | fH3Pi0AsymPt = new TH3F(hname,htitle,20,0,1,20,0,10,200,0,1); | |
212 | ||
213 | fCreateHist = true; | |
214 | return true; | |
215 | } | |
216 | //------------------------------------------------------------------- | |
217 | void AliPHOSDApi0mip::FillTree(AliPHOSDATreeEvent* event){ | |
218 | // Fill all information (AliPHOSDATreeEvent) into TTree | |
219 | // if event==NULL, the obtained information (=fEvent) is used. | |
220 | ||
221 | if(! fCreateTree ) CreateTree(); | |
222 | if( event == 0 ){ | |
223 | if( fEvent == 0 ){ | |
224 | fEvent = new AliPHOSDATreeEvent(); | |
225 | } | |
226 | // No need to make clustering for tree saving. | |
227 | //if(! fEventClustered ){ | |
228 | //fEvent->ExecuteClustering(); | |
229 | //fEventClustered = true; | |
230 | //} | |
231 | fEvent->SetTime(fTime); | |
232 | fTTree->Fill(); | |
233 | } else { | |
234 | if( fEvent ) delete fEvent; | |
235 | fEvent = event; | |
236 | fTTree->Fill(); | |
237 | fEvent = 0; | |
238 | } | |
239 | } | |
240 | //------------------------------------------------------------------- | |
241 | void AliPHOSDApi0mip::FillHist(AliPHOSDATreeEvent* event){ | |
242 | // Fill all information (AliPHOSDATreeEvent) into histograms | |
243 | // if event==NULL, the obtained information (=fEvent) is used. | |
244 | ||
245 | if(! fCreateHist ) CreateHist(); | |
246 | if( event == 0 ){ | |
247 | if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent; | |
248 | if(! fEventClustered ){ | |
249 | fEvent->ExecuteClustering(); | |
250 | fEventClustered = true; | |
251 | } | |
252 | fEvent->SetTime(fTime); | |
253 | event = fEvent; | |
254 | } | |
255 | ||
256 | // -- Constant | |
257 | int nclt1, nclt2; | |
258 | int ndigits; | |
259 | int ndigitsall = 0; | |
260 | AliPHOSDATreeDigit digit1, digit2; | |
261 | const float kDistanceFromVertex = 460.; | |
262 | float dist1, dist2, px, py; | |
263 | float mass, cosine, pt, asym; | |
264 | ||
265 | // -- Filling into Time Histogram | |
266 | double start = fH1Time->GetBinContent(1); | |
267 | double stop = fH1Time->GetBinContent(2); | |
268 | double currenttime = (double) (event->GetTime()); | |
269 | if( start == 0 && stop == 0 ) start = stop = currenttime; | |
270 | if( currenttime < start ) start = currenttime; | |
271 | if( currenttime > stop ) stop = currenttime; | |
272 | fH1Time->SetBinContent(1,start); | |
273 | fH1Time->SetBinContent(2,stop); | |
274 | ||
275 | // -- Start looping clusters | |
276 | nclt1 = event->GetNClusters(); | |
277 | fH1ClusterNum->Fill(nclt1); | |
278 | //std::cout<<" AliPHOSDApi0mip::FillHisto() DEBUG:: Number of clusters: "<<nclt1<<std::endl; | |
279 | while( nclt1-- ){ | |
280 | AliPHOSDATreeCluster& clt1 = event->GetCluster(nclt1); | |
281 | ndigits = clt1.GetNDigits(); | |
282 | while( ndigits-- ){ | |
283 | AliPHOSDATreeDigit& digit = clt1.GetDigit(ndigits); | |
284 | fH2EneDigitId->Fill(digit.GetDigitId(),digit.GetEnergy()); | |
285 | ndigitsall++; | |
286 | } | |
287 | digit1 = clt1.GetMaxDigit(); | |
288 | if( clt1.GetEnergy() > 0 && clt1.GetNDigits() >= 2 ){ | |
289 | fH2MipDigitId->Fill(digit1.GetDigitId(),clt1.GetEnergy()); | |
290 | } | |
291 | if( clt1.GetEnergy() > 0.2 ){ | |
292 | nclt2 = nclt1; | |
293 | while( nclt2-- ){ | |
294 | AliPHOSDATreeCluster& clt2 = event->GetCluster(nclt2); | |
295 | digit2 = clt2.GetMaxDigit(); | |
296 | if( clt2.GetEnergy() > 0.2 ){ | |
297 | // | |
298 | // Following approximation is applied in order to avoid many calling "sqrt" function. | |
299 | // sqrt( 1 + x^2 ) =~ 1 + 0.5x^2 (for x<<1) | |
300 | // | |
301 | //cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex) | |
302 | // / ( kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*2.42/kDistanceFromVertex ) | |
303 | // / ( kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*2.42/kDistanceFromVertex ); | |
304 | // | |
305 | ||
306 | ||
307 | dist1 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*4.84 ); | |
308 | dist2 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 ); | |
309 | cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex) | |
310 | / dist1 / dist2; | |
311 | mass = sqrt( 2.0 * clt1.GetEnergy() * clt2.GetEnergy() * (1 - cosine) ); | |
312 | //std::cout<<" DEBUG: dist1:"<<dist1<<" dist2:"<<dist2<<" cosine:"<<cosine<<" c1:"<<clt1.GetEnergy()<<" c2:"<<clt2.GetEnergy()<<" mass:"<<mass<<std::endl; | |
313 | //std::cout<<" : clt1(row,col)=("<<clt1.GetRow()<<","<<clt1.GetCol()<<") clt2(row,col)=("<<clt2.GetRow()<<","<<clt2.GetCol()<<") "<<std::endl; | |
314 | asym = fabs(clt1.GetEnergy() - clt2.GetEnergy())/(clt1.GetEnergy() + clt2.GetEnergy()); | |
315 | py = clt1.GetEnergy() * kDistanceFromVertex / dist1 + clt2.GetEnergy() * kDistanceFromVertex / dist2; | |
316 | px = clt1.GetEnergy() * (clt1.GetCol()-27.5)*2.42/dist1 + clt2.GetEnergy() * (clt2.GetCol()-27.5)*2.42/ dist2; | |
317 | pt = sqrt(px*px+py*py); | |
318 | fH3Pi0AsymPt->Fill(asym,pt,mass); | |
319 | if( digit1.IsValid() && digit1.IsValid() && | |
320 | digit1.GetDigitId() != digit2.GetDigitId() && | |
321 | //pt > 1.5 && clt1.GetEnergy() > 0.3 && clt2.GetEnergy() > 0.3 ){ | |
322 | //pt > 1.3 && clt1.GetEnergy() > 0.25 && clt2.GetEnergy() > 0.25 ){ | |
323 | pt > 1.0 && clt1.GetEnergy() > 0.2 && clt2.GetEnergy() > 0.2 ){ | |
324 | fH2Pi0DigitId->Fill(digit1.GetDigitId(),mass); | |
325 | fH2Pi0DigitId->Fill(digit2.GetDigitId(),mass); | |
326 | } | |
327 | } // End of energy cut | |
328 | } // End of second cluster loop (clt2) | |
329 | } // End of energy cut | |
330 | } // End of first cluster loop (clt1) | |
331 | fH1DigitNum->Fill(ndigitsall); | |
332 | // | |
333 | } | |
334 | //------------------------------------------------------------------- | |
c335b3c5 | 335 | void AliPHOSDApi0mip::Print(Option_t *option) const |
336 | { | |
12f6dd57 | 337 | // Print Out |
338 | ||
339 | //fTFile->ls(); | |
340 | //fTTree->Print(); | |
341 | if( fEvent ){ | |
342 | //fEvent->ExecuteClustering(); | |
c335b3c5 | 343 | fEvent->Print(option); |
12f6dd57 | 344 | } |
345 | } | |
346 | //------------------------------------------------------------------- |