typo fixed: HMP_DET/HMP_ENV/HMP_ENV_TENV.actual.value
[u/mrichter/AliRoot.git] / T0 / T0da.cxx
CommitLineData
310b89b2 1extern "C" {
2#include <daqDA.h>
3}
4
5#include "event.h"
6#include "monitor.h"
7
8#include <Riostream.h>
9#include <stdio.h>
10#include <stdlib.h>
11
12//AliRoot
13#include "AliRawReaderDate.h"
14#include "AliRawReader.h"
15#include "AliT0digit.h"
16#include "AliT0RawReader.h"
17#include "AliT0Dqclass.h"
18
19//ROOT
20#include "TFile.h"
21#include "TKey.h"
22#include "TH2S.h"
23#include "TObject.h"
24#include "TBenchmark.h"
25#include "TRandom.h"
26
27#include "TCanvas.h"
28#include "TString.h"
29#include "TH1.h"
30#include "TF1.h"
31#include "TSpectrum.h"
32#include "TVirtualFitter.h"
33
34void fitv2DA(char *filename = "t0histdate.root");
35
36/* Main routine
37 Arguments:
38 1- monitoring data source
39*/
40int main(int argc, char **argv) {
41 int status;
42
43 if (argc!=2) {
44 printf("Wrong number of arguments\n");
45 return -1;
46 }
47
48
49 /* define data source : this is argument 1 */
50 status=monitorSetDataSource( argv[1] );
51 if (status!=0) {
52 printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
53 return -1;
54 }
55
56
57 /* declare monitoring program */
58 status=monitorDeclareMp( __FILE__ );
59 if (status!=0) {
60 printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
61 return -1;
62 }
63
64
65 /* define wait event timeout - 1s max */
66 monitorSetNowait();
67 monitorSetNoWaitNetworkTimeout(1000);
68
69
70 /* log start of process */
71 printf("T0 monitoring program started\n");
72
73
74 // Allocation of histograms - start
75 TH1F * hCFD[24]; TH1F *hLED[24]; TH1F*hQT01[4];TH1F*hQT02[4];
76 TH1F*hQTD[4]; TH1F*hADC[24]; TH2*hQTCFD[4];
77 Char_t buf1[10], buf2[10], buf3[10], buf4[10], buf5[10], buf6[10],buf7[10];
78 for (Int_t ic=0; ic<24; ic++)
79 {
80 sprintf(buf1,"CFD%i",ic+1);
81 hCFD[ic]= new TH1F(buf1,"CFD",6000,0,6000);
82 sprintf(buf2,"LED%i",ic+1);
83 hLED[ic]= new TH1F(buf2,"LED",2000,0,2000);
84 //LED-CFD
85 sprintf(buf6,"ADC%i",ic+1);
86 hADC[ic]= new TH1F(buf6,"ADC",6000,0,6000);
87 }
88
89 for (Int_t iq=0; iq<4; iq++)
90 {
91 //QT01 - QT04
92 sprintf(buf3,"QT0%i",iq+1);
93 hQT01[iq]= new TH1F(buf3,"QT01",6000,0,6000);
94 sprintf(buf4,"QT1_%i",iq+1);
95 //QT11 - QT14
96 hQT02[iq]= new TH1F(buf4,"QT02",6000,0,6000);
97 sprintf(buf5,"QTD_%i",iq+1);
98 //QT11-QT01 ....
99 hQTD[iq]= new TH1F(buf5,"QTdiff",4500,1500,6000);
100
101 sprintf(buf7,"QTCFD_%i",iq+1);
102 hQTCFD[iq]= new TH2F(buf7,"QT vs CFD",500,0,6000,500,0,5000);
103 }
104
105 TH1F *hORA= new TH1F("hORA"," T0 A ",1000, 1000,2000);
106 TH1F *hORC= new TH1F("hORC"," T0 C ",1000, 1000,2000);
107 TH1F*hEffCFD= new TH1F("hEffCFD","Effeciency",8,0.25,4.25);
108
109 // TH2F*hQTCFD= new TH2F("hQTCFD","QT vs CFD",500,0.5,6000.5,500,0.5,5000.5);
110 TH2F*hQTLED= new TH2F("hQTLED","QT vs LED",500,0.5,6000.5,500,0.5,5000.5);
111 TH2F*hLEDCFD= new TH2F("hLEDCFD","LEd vs CFD",500,-0.5,10000.5,500,-0.5,10000.5);
112 // Allocation of histograms - end
113
114 Int_t iev=0;
115 /* main loop (infinite) */
116 for(;;) {
117 struct eventHeaderStruct *event;
118 eventTypeType eventT;
119
120 /* check shutdown condition */
121 if (daqDA_checkShutdown()) {break;}
122
123 /* get next event (blocking call until timeout) */
124 status=monitorGetEventDynamic((void **)&event);
125 if (status==(int)MON_ERR_EOF) {
126 printf ("End of File detected\n");
127 break; /* end of monitoring file has been reached */
128 }
129
130 if (status!=0) {
131 printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
132 break;
133 }
134
135 /* retry if got no event */
136 if (event==NULL) {
137 continue;
138 }
139
140 iev++;
141
142 /* use event - here, just write event id to result file */
143 eventT=event->eventType;
144
145 if (eventT==PHYSICS_EVENT) {
146 printf(" event number = %i \n",iev);
147
148 // Initalize raw-data reading and decoding
149 AliRawReader *reader = new AliRawReaderDate((void*)event);
150 // Enable the following two lines in case of real-data
151 // reader->LoadEquipmentIdsMap("T0map.txt");
152 // reader->RequireHeader(kFALSE);
153 AliT0RawReader *start = new AliT0RawReader(reader);
154
155 // Read raw data
156 Int_t allData[110][5];
157 start->Next();
158 for (Int_t i=0; i<110; i++) {
159 allData[i][0]= start->GetData(i,0);
160 if (allData[i][0] != 0) cout<<"event "<<event<<" i "<<iev<<" "<<allData[i][0] - allData[0][0]<<endl;
161 }
162
163 // Fill the histograms
164 for(Int_t ik=0; ik<24; ik++) {
165 for (Int_t iHit=0; iHit<5; iHit++) {
166 if(allData[ik+1][iHit] != 0 ) {
167 // cout<<"event "<<event<<" iHit "<< iHit<< allData[ik+1][0]<<" "<< allData[0][0]<< endl;
168 hCFD[ik]->Fill((allData[ik+1][iHit] - allData[0][0])/1000.);
169 }
170 if(allData[ik+25][iHit] != 0 ) hLED[ik]->Fill(allData[ik+25][iHit] - allData[0][0]);
171 if(allData[ik+25][iHit] != 0 || allData[ik+1][iHit] !=0)
172 hADC[ik]->Fill(allData[ik+25][iHit] - allData[ik+1][iHit]);
173 }
174 }
175 hLEDCFD->Fill(allData[9][0],allData[1][0]);
176
177 hQT01[1]->Fill(allData[18][0] - allData[0][0]);
178 hQT01[2]->Fill(allData[19][0] - allData[0][0]);
179 hQT01[3]->Fill(allData[20][0] - allData[0][0]);
180 hQT02[0]->Fill(allData[21][0] - allData[0][0]);
181 hQT02[0]->Fill(allData[22][0] - allData[0][0]);
182 hQT02[0]->Fill(allData[23][0] - allData[0][0]);
183 hQT02[0]->Fill(allData[24][0] - allData[0][0]);
184
185
186 hQTD[0]->Fill(allData[21][0] - allData[17][0]);
187 hQTLED->Fill(allData[9][0] - allData[0][0], allData[21][0]-allData[17][0]);
188 hQTCFD[0]->Fill(allData[21][0]-allData[17][0], allData[1][0] - allData[0][0]);
189
190 hQTD[1]->Fill(allData[22][0]-allData[18][0]);
191 hQTCFD[1]->Fill(allData[22][0]-allData[18][0], allData[2][0] - allData[0][0]);
192
193 hQTD[2]->Fill(allData[23][0]-allData[19][0]);
194 hQTCFD[2]->Fill(allData[23][0]-allData[19][0], allData[3][0] - allData[0][0]);
195
196 hQTD[3]->Fill(allData[24][0]-allData[20][0]);
197 hQTCFD[3]->Fill(allData[24][0]-allData[20][0], allData[4][0] - allData[0][0]);
198
199
200 hORA->Fill(allData[25][0] - allData[0][0]);
201 hORA->Fill(allData[26][0] - allData[0][0]);
202 // End of fill histograms
203
204 }
205
206
207
208 /* free resources */
209 free(event);
210
211 /* exit when last event received, no need to wait for TERM signal */
212 if (eventT==END_OF_RUN) {
213 printf("EOR event detected\n");
214 break;
215 }
216 }
217
218 // write a file with the histograms
219 Char_t filehist[20];
220 sprintf(filehist,"t0histdate.root");
221 TFile *hist = new TFile(filehist,"RECREATE");
222
223 hLEDCFD->Write();
224
225 for (Int_t i=0; i<24; i++)
226 {
227 hCFD[i]->Write();
228 hLED[i]->Write();
229 hADC[i]->Write();
230 }
231 for (Int_t i=0; i<4; i++)
232 {
233 hQT01[i]->Write();
234 hQT02[i]->Write();
235 hQTD[i]->Write();
236 hQTCFD[i]->Write();
237 }
238
239 hEffCFD->Write();
240 hORA->Write();
241 hORC->Write();
242 hQTLED->Write();
243 hist->Close();
244
245 // Fit the histograms and write the output file
246 fitv2DA();
247
248 return status;
249}
250
251void fitv2DA(char *filename)
252{
253 TFile *file = TFile::Open(filename);
254 FILE *dafile;
255 dafile = fopen("da.txt","w");
256 fprintf(dafile,"Const\t\tMean\t\tSigma\t\tLBorder\tRBorder\n");
257
258// TFile *fithist = new TFile("fitedhist.root","RECREATE");
259
260 char histoname[10];
261 char canvname[10];
262 Float_t p[3] = {0.,0.,0.};
263 Float_t LBordX = 0.;
264 Float_t RBordX = 0.;
265/*
266 TIter next(file.GetListOfKeys());
267 TKey *key;
268 while ((key = (TKey*)nextkey())) {
269 if(key->GetName()==)
270 }
271*/
272 Int_t npeaks = 10;
273 cout<<"npeaks = "<<npeaks<<endl;
274 TH1F *hist[24];
275 TCanvas *c[24];
276 for(Int_t i=0;i<24;i++){///////////#######Cos tu sie dzieje....
277 sprintf(histoname,"CFD%d",i+1);
278 sprintf(canvname,"C%d",i+1);
279 TH1F *histtemp = (TH1F*)file->Get(histoname);
280 hist[i] = histtemp;
281 }
282 cout<<"Wczytany plik z wykresami..."<<endl;
283// TCanvas *c[8];
284 for(Int_t ii=0;ii<24;ii++){
285 c[ii] = new TCanvas(canvname,canvname,800,600);
286 p[0] = 0.;
287 p[1] = 0.;
288 p[2] = 0.;
289 LBordX = 0.;
290 RBordX = 0.;
291
292 Double_t max=0.0;
293// Double_t peak,peakmax=0.0;
294 Double_t tabmax[2] = {0.0, 0.0};
295
296 TSpectrum *s = new TSpectrum(2*npeaks);
297 Int_t nfound = s->Search(hist[ii],2,"goff",0.05);
298 cout<<"Found "<<nfound<<" peaks\n";
299 if(nfound!=0){
300 Float_t *xpeak = s->GetPositionX();
301 for(Int_t k=0;k<nfound;k++)
302 {
303 cout<<"Petla po pikach..."<<endl;
304 Float_t xp = xpeak[k];
305 Int_t xbin = hist[ii]->GetXaxis()->FindBin(xp);
306 Float_t yp = hist[ii]->GetBinContent(xbin);
307 if(yp>max) {
308 max = yp;
309 tabmax[0] = xp;
310 tabmax[1] = yp;
311 }
312 cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[k]<<"\typeak = "<<yp<<endl;
313 }
314
315 Double_t MinBordX, MaxBordX, MaxBordY;
316
317 MinBordX = tabmax[0]-tabmax[0]*0.05;
318 MaxBordX = tabmax[0]+tabmax[0]*0.05;
319 MaxBordY = tabmax[1]+tabmax[1]*0.5;
320
321 Double_t border;
322 border = tabmax[1]*0.1;
323
324 Double_t LBordY = 99999.;
325 Double_t RBordY = 99999.;
326 LBordX = tabmax[0];
327 RBordX = tabmax[0];
328 do{
329 LBordY = hist[ii]->GetBinContent(LBordX);
330 LBordX -=1;
331 }while(LBordY>border);
332 do{
333 RBordX +=1;
334 RBordY = hist[ii]->GetBinContent(RBordX);
335 }while(RBordY>border);
336
337
338 hist[ii]->GetXaxis()->SetRangeUser(MinBordX,MaxBordX);
339 hist[ii]->GetYaxis()->SetRangeUser(0,MaxBordY);
340 TF1 *fit = new TF1("fit","gaus",LBordX,RBordX);
341 hist[ii]->Fit("fit","IR");
342 Double_t maxx=0.;
343// Double_t chi2 = fit->GetChisquare();
344// Float_t p[3];
345
346 for(Int_t j=0;j<3;j++){
347 p[j] = fit->GetParameter(j);
348 }
349
350 max = fit->GetMaximum();
351 maxx = fit->GetMaximumX();
352// hist[ii]->Write();
353 cout<<"\nFitMaxY = "<<max<<"\tFitMaxX "<<maxx<<"\tPeakMax*10% "<<border<<endl;
354 cout<<"LBordY "<<LBordY<<"\tLBordX "<<LBordX+1<<"\tRBordY "<<RBordY<<"\tRBordX "<<RBordX<<endl;
355 cout<<"Parameters: Const - MaxY - "<<p[0]<<"\tMean - "<<p[1]<<"\tSigma - "<<p[2]<<endl;
356 cout<<"Borders: Min - "<<MinBordX<<"\tMax - "<<MaxBordX<<endl;
357 if(p[0]==0.&&p[1]==0.&&p[2]==0.) RBordX = 0;
358 }
359// else{
360// return;
361// }
362 AliT0Dqclass *daqpar = new AliT0Dqclass();
363 cout<<"Alit0dqclass..."<<endl;
364 daqpar->SetTime(0,p[0]);
365 cout<<"First parameter..."<<endl;
366 daqpar->SetTime(1,p[1]);
367 daqpar->SetTime(2,p[2]);
368 cout<<"gauss parameters finished..."<<endl;
369 daqpar->SetTime(3,LBordX+1);
370 cout<<"first border..."<<endl;
371 daqpar->SetTime(4,RBordX);
372 cout<<"second border..."<<endl;
373 TFile *DA = new TFile("DAfile.root","RECREATE");
374 DA->cd();
375 cout<<"file"<<endl;
376 daqpar->Write("Time");
377 cout<<"file written"<<endl;
378 DA->Close();
379 cout<<"file closed"<<endl;
380 delete DA;
381 cout<<"file deleted"<<endl;
382
383 fprintf(dafile,"%lf\t%lf\t%lf\t%d\t%d\n",p[0],p[1],p[2],LBordX+1,RBordX);
384 }
385
386
387 fclose(dafile);
388}