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