Cleanup of collisions geometries and headers.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSDApi0mip.cxx
CommitLineData
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"
30ClassImp(AliPHOSDApi0mip)
31//----------------------------------------------------------------
32AliPHOSDApi0mip::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//----------------------------------------------------------------
58AliPHOSDApi0mip::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 104AliPHOSDApi0mip& AliPHOSDApi0mip::operator= (const AliPHOSDApi0mip&)
105{
106 // Operator= is not implemented yet
107 AliFatal("Operator= is not implemented");
108 return *this;
109}
110//-------------------------------------------------------------------
12f6dd57 111AliPHOSDApi0mip::~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//-------------------------------------------------------------------
139void 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//-------------------------------------------------------------------
151void AliPHOSDApi0mip::NewEvent(){
152 // Initiator must be called in the beginning of event
153 //
154
155 if(fEvent) fEvent->Reset();
156 fEventClustered = false;
157}
158//----------------------------------------------------------------
159bool 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//----------------------------------------------------------------
173bool 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//-------------------------------------------------------------------
217void 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//-------------------------------------------------------------------
241void 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 335void 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//-------------------------------------------------------------------