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