Fast implementation of gamma-distributed random numbers (Andreas)
[u/mrichter/AliRoot.git] / ZDC / ZDCMAPPINGda.cxx
CommitLineData
76493b40 1/*\r
2\r
3This program reads the DAQ data files passed as argument using the monitoring library.\r
4\r
5The program reports about its processing progress.\r
6\r
7Messages on stdout are exported to DAQ log system.\r
8\r
9DA to write mapping for ADC modules and VME scaler\r
10\r
11Contact: Chiara.Oppedisano@to.infn.it\r
12Link: \r
13Run Type: PHYSICS\r
14DA Type: MON\r
15Number of events needed: at least 50% \r
16Input Files: none\r
17Output Files: ZDCChMapping.dat\r
18Trigger Types Used: different trigger types are used\r
19\r
20*/\r
21\r
22#define MAPDATA_FILE "ZDCChMapping.dat"\r
23#define TDCDATA_FILE "ZDCTDCCalib.dat"\r
24#define TDCHISTO_FILE "ZDCTDCHisto.root"\r
25\r
26#include <stdio.h>\r
27#include <stdlib.h>\r
28#include <Riostream.h>\r
29#include <signal.h>\r
30\r
31// DATE\r
32#include <event.h>\r
33#include <monitor.h>\r
34#include <daqDA.h>\r
35\r
36//ROOT\r
37#include <TROOT.h>\r
38#include <TPluginManager.h>\r
39#include <TH1F.h>\r
40#include <TH2F.h>\r
41#include <TProfile.h>\r
42#include <TF1.h>\r
43#include <TFile.h>\r
44#include <TFitter.h>\r
45#include "TMinuitMinimizer.h"\r
46\r
47//AliRoot\r
48#include <AliRawReaderDate.h>\r
49#include <AliRawEventHeaderBase.h>\r
50#include <AliZDCRawStream.h>\r
51\r
52int main(int argc, char **argv) {\r
53\r
45b889f8 54 // needed for streamer application\r
76493b40 55 gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",\r
45b889f8 56 "*",\r
57 "TStreamerInfo",\r
58 "RIO",\r
59 "TStreamerInfo()"); \r
60 // needed for Minuit plugin\r
61 gROOT->GetPluginManager()->AddHandler("ROOT::Math::Minimizer",\r
62 "Minuit",\r
63 "TMinuitMinimizer",\r
64 "Minuit",\r
65 "TMinuitMinimizer(const char*)");\r
76493b40 66\r
45b889f8 67 //TMinuitMinimizer m; \r
76493b40 68 TVirtualFitter::SetDefaultFitter("Minuit");\r
69\r
76493b40 70 char *monitor_table[] = { "ALL", "no", "PHY", "yes", "SOD", "all", NULL };\r
71 int err = monitorDeclareTable(monitor_table);\r
72 if(err){\r
73 printf("monitorDeclareTable() failed: %s\n", monitorDecodeError(err));\r
74 return -1;\r
75 } \r
76 \r
77 int status=0, nphys=0;\r
ed8a15fd 78 int const kNModules = 9;\r
76493b40 79 int const kNChannels = 24;\r
80 int const kNScChannels = 32; \r
616b80bc 81 int const kZDCTDCGeo = 4;\r
76493b40 82 \r
83 int itdc=0, iprevtdc=-1, ihittdc=0;\r
84 float tdcData[6], tdcL0=-999.; \r
85 for(int ij=0; ij<6; ij++) tdcData[ij]=-999.;\r
86 \r
87 Int_t iMod=-1;\r
88 Int_t modGeo[kNModules], modType[kNModules],modNCh[kNModules];\r
89 for(Int_t kl=0; kl<kNModules; kl++){\r
90 modGeo[kl]=modType[kl]=modNCh[kl]=0;\r
91 }\r
92 \r
93 Int_t ich=0;\r
94 Int_t adcMod[2*kNChannels], adcCh[2*kNChannels], sigCode[2*kNChannels];\r
95 Int_t det[2*kNChannels], sec[2*kNChannels];\r
96 for(Int_t y=0; y<2*kNChannels; y++){\r
97 adcMod[y]=adcCh[y]=sigCode[y]=det[y]=sec[y]=0;\r
98 }\r
99 \r
100 Int_t iScCh=0;\r
101 Int_t scMod[kNScChannels], scCh[kNScChannels], scSigCode[kNScChannels];\r
102 Int_t scDet[kNScChannels], scSec[kNScChannels];\r
103 for(Int_t y=0; y<kNScChannels; y++){\r
104 scMod[y]=scCh[y]=scSigCode[y]=scDet[y]=scSec[y]=0;\r
105 }\r
106 \r
107 Int_t itdcCh=0;\r
108 Int_t tdcMod[kNScChannels], tdcCh[kNScChannels], tdcSigCode[kNScChannels];\r
109 Int_t tdcDet[kNScChannels], tdcSec[kNScChannels];\r
110 for(Int_t y=0; y<kNScChannels; y++){\r
111 tdcMod[y]=tdcCh[y]=tdcSigCode[y]=tdcDet[y]=tdcSec[y]=-1;\r
112 }\r
113 \r
616b80bc 114 TH1F * hTDC[6]={0x0,0x0,0x0,0x0,0x0,0x0};\r
45b889f8 115 for(int it=0; it<6; it++){\r
116 if(it==0) hTDC[it] = new TH1F("TDCZEM1","TDC ZEM1",200, -200., 200.);\r
117 else if(it==1) hTDC[it] = new TH1F("TDCZEM2","TDC ZEM2",200, -200., 200.);\r
118 else if(it==2) hTDC[it] = new TH1F("TDCZNC", "TDC ZNC", 200, -200., 200.);\r
119 else if(it==3) hTDC[it] = new TH1F("TDCZPC", "TDC ZPC", 200, -200., 200.);\r
120 else if(it==4) hTDC[it] = new TH1F("TDCZNA", "TDC ZNA", 200, -200., 200.);\r
121 else if(it==5) hTDC[it] = new TH1F("TDCZPA", "TDC ZPA", 200, -200., 200.);\r
76493b40 122 }\r
123 \r
124 /* log start of process */\r
125 printf("\n ZDC MAPPING program started\n"); \r
126 signal(SIGSEGV, SIG_DFL);\r
127\r
128 /* check that we got some arguments = list of files */\r
45b889f8 129 if(argc<2) {\r
76493b40 130 printf("Wrong number of arguments\n");\r
131 return -1;\r
132 }\r
133 \r
134 FILE *mapFile4Shuttle;\r
135\r
136 /* read the data files */\r
45b889f8 137 for(int n=1;n<argc;n++){\r
76493b40 138 \r
139 status=monitorSetDataSource( argv[n] );\r
140 if (status!=0) {\r
141 printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));\r
142 return -1;\r
143 }\r
144 \r
145 /* declare monitoring program */\r
146 status=monitorDeclareMp( __FILE__ );\r
147 if (status!=0) {\r
148 printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));\r
149 return -1;\r
150 }\r
151 monitorSetNowait();\r
152 monitorSetNoWaitNetworkTimeout(1000);\r
153\r
154\r
155 int iev = 0;\r
156 //Bool_t sodRead = kFALSE;\r
157 //while(!sodRead){\r
158 \r
159 /* loop on events (infinite) */\r
160 for(;;) {\r
222fabe9 161 \r
162 if(nphys > 50000) break;\r
163 \r
76493b40 164 struct eventHeaderStruct *event;\r
165 eventTypeType eventT;\r
166 \r
167 /* check shutdown condition */\r
168 if (daqDA_checkShutdown()) {break;}\r
169\r
170 /* get next event */\r
171 status=monitorGetEventDynamic((void **)&event);\r
172 if(status==MON_ERR_EOF){\r
173 printf ("End of File detected\n");\r
174 break; /* end of monitoring file has been reached */\r
175 }\r
176 if(status!=0) {\r
177 printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));\r
178 return -1;\r
179 }\r
180\r
181 /* retry if got no event */\r
182 if(event==NULL) continue;\r
183 \r
184 // Initalize raw-data reading and decoding\r
185 AliRawReader *reader = new AliRawReaderDate((void*)event);\r
186 reader->Reset();\r
187 reader->Select("ZDC");\r
188 // --- Reading event header\r
189 //UInt_t evtype = reader->GetType();\r
190 //printf("\t ZDCMAPPINGda -> run # %d\n",reader->GetRunNumber());\r
191 //\r
192 AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);\r
193 \r
194 \r
195 /* use event - here, just write event id to result file */\r
196 eventT=event->eventType;\r
197 \r
198 if(eventT==START_OF_DATA){\r
199 \r
200 rawStreamZDC->SetSODReading(kTRUE);\r
201 \r
202 // --------------------------------------------------------\r
203 // --- Writing ascii data file for the Shuttle preprocessor\r
204 mapFile4Shuttle = fopen(MAPDATA_FILE,"w");\r
205 //if(rawStreamZDC->Next()){\r
206 while((rawStreamZDC->Next())){\r
207 if(rawStreamZDC->IsHeaderMapping()){ // mapping header\r
208 iMod++;\r
209 modGeo[iMod] = rawStreamZDC->GetADCModule();\r
210 modType[iMod] = rawStreamZDC->GetModType();\r
211 modNCh[iMod] = rawStreamZDC->GetADCNChannels();\r
212 }\r
213 if(rawStreamZDC->IsChMapping()){ \r
214 if(modType[iMod]==1){ // ADC mapping ----------------------\r
215 adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);\r
216 adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);\r
217 sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);\r
218 det[ich] = rawStreamZDC->GetDetectorFromMap(ich);\r
219 sec[ich] = rawStreamZDC->GetTowerFromMap(ich);\r
220 ich++;\r
221 }\r
222 else if(modType[iMod]==2){ // VME scaler mapping -------------\r
223 scMod[iScCh] = rawStreamZDC->GetScalerModFromMap(iScCh);\r
224 scCh[iScCh] = rawStreamZDC->GetScalerChFromMap(iScCh);\r
225 scSigCode[iScCh] = rawStreamZDC->GetScalerSignFromMap(iScCh);\r
226 scDet[iScCh] = rawStreamZDC->GetScDetectorFromMap(iScCh);\r
227 scSec[iScCh] = rawStreamZDC->GetScTowerFromMap(iScCh);\r
228 iScCh++;\r
229 }\r
230 else if(modType[iMod]==6 && modGeo[iMod]==4){ // ZDC TDC mapping --------------------\r
231 tdcMod[itdcCh] = rawStreamZDC->GetTDCModFromMap(itdcCh);\r
232 tdcCh[itdcCh] = rawStreamZDC->GetTDCChFromMap(itdcCh);\r
233 tdcSigCode[itdcCh] = rawStreamZDC->GetTDCSignFromMap(itdcCh);\r
234 itdcCh++;\r
235 }\r
236 }\r
237 }\r
238 // Writing data on output FXS file\r
239 for(Int_t is=0; is<2*kNChannels; is++){\r
240 fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",\r
241 is,adcMod[is],adcCh[is],sigCode[is],det[is],sec[is]);\r
242 //printf(" Mapping DA -> %d ADC: mod %d ch %d, code %d det %d, sec %d\n",\r
ed8a15fd 243 //is,adcMod[is],adcCh[is],sigCode[is],det[is],sec[is]);\r
76493b40 244 }\r
245 for(Int_t is=0; is<kNScChannels; is++){\r
246 fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",\r
247 is,scMod[is],scCh[is],scSigCode[is],scDet[is],scSec[is]);\r
248 //if(scMod[is]!=-1) printf(" Mapping DA -> %d Scaler: mod %d ch %d, code %d det %d, sec %d\n",\r
ed8a15fd 249 //is,scMod[is],scCh[is],scSigCode[is],scDet[is],scSec[is]);\r
76493b40 250 }\r
251 for(Int_t is=0; is<kNScChannels; is++){\r
252 fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\n",\r
253 is,tdcMod[is],tdcCh[is],tdcSigCode[is]);\r
254 //if(tdcMod[is]!=-1) printf(" Mapping DA -> %d TDC: mod %d ch %d, code %d\n",\r
ed8a15fd 255 //is,tdcMod[is],tdcCh[is],tdcSigCode[is]);\r
76493b40 256 }\r
257 for(Int_t is=0; is<kNModules; is++){\r
258 fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\n",\r
259 modGeo[is],modType[is],modNCh[is]);\r
260 //printf(" Mapping DA -> Module mapping: geo %d type %d #ch %d\n",\r
ed8a15fd 261 //modGeo[is],modType[is],modNCh[is]);\r
76493b40 262 }\r
263 \r
264 //} //if (rawstream)\r
265 fclose(mapFile4Shuttle);\r
266 }// SOD event\r
267 else if(eventT==PHYSICS_EVENT){ \r
268\r
269 for(int ij=0; ij<6; ij++) tdcData[ij]=-999.;\r
270 rawStreamZDC->SetSODReading(kTRUE);\r
271\r
272 // ----- Setting ch. mapping -----\r
45b889f8 273 for(int jk=0; jk<2*kNChannels; jk++){\r
76493b40 274 rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);\r
275 rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);\r
276 rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);\r
277 rawStreamZDC->SetMapDet(jk, det[jk]);\r
278 rawStreamZDC->SetMapTow(jk, sec[jk]);\r
279 }\r
280 \r
281 while(rawStreamZDC->Next()){\r
616b80bc 282 if(rawStreamZDC->GetADCModule()!=kZDCTDCGeo) continue; //skipping ADCs, scalers and trigger cards\r
283\r
284 if(rawStreamZDC->GetADCModule()==kZDCTDCGeo && rawStreamZDC->IsZDCTDCDatum()==kTRUE){\r
76493b40 285 //\r
286 itdc = rawStreamZDC->GetChannel(); \r
616b80bc 287 if(itdc==iprevtdc) ihittdc++;\r
288 else ihittdc=0;\r
289 iprevtdc=itdc;\r
290 //\r
45b889f8 291 if(((itdc>=8 && itdc<=13) || itdc==15) && ihittdc==0){\r
616b80bc 292 if(itdc!=15){\r
293 tdcData[itdc-8] = 0.025*rawStreamZDC->GetZDCTDCDatum();\r
294 //\r
45b889f8 295 //printf(" ev.%d **** TDC%d %1.0f ns \n",nphys,itdc, tdcData[itdc-8]);\r
616b80bc 296 }\r
297 //\r
298 else if(itdc==15){\r
299 tdcL0 = 0.025*rawStreamZDC->GetZDCTDCDatum();\r
300 //\r
45b889f8 301 //printf(" ev.%d ++++ TDCL0 %1.0f ns \n",nphys,tdcL0);\r
616b80bc 302 //\r
303 for(int ic=0; ic<6; ic++){\r
304 if(tdcData[ic]!=-999. && tdcL0!=-999.){\r
305 hTDC[ic]->Fill(tdcData[ic]-tdcL0);\r
45b889f8 306 //printf(" ev.%d -> Filling histo%d: %f ns\n",nphys,ic, tdcData[ic]-tdcL0);\r
616b80bc 307 }\r
308 }\r
309 }\r
76493b40 310 }\r
311 }\r
312 }\r
313 \r
45b889f8 314 nphys++;\r
315 \r
316 delete rawStreamZDC;\r
317 rawStreamZDC = 0x0; \r
d5be89ef 318 delete reader;\r
45b889f8 319\r
76493b40 320 }//(if PHYSICS_EVENT) \r
321 else if(eventT==END_OF_RUN){\r
322 printf("End Of Run detected\n");\r
323 break;\r
324 }\r
325 \r
76493b40 326 iev++; \r
327\r
328 /* free resources */\r
329 free(event);\r
330 \r
331 } // event loop \r
332 \r
333 }\r
334 printf(" \n # of physics events analyzed = %d\n\n", nphys);\r
335\r
336 /* Analysis of the histograms */\r
337 //\r
338 FILE *fileShuttle;\r
339 fileShuttle = fopen(TDCDATA_FILE,"w");\r
340 //\r
341 Float_t xUp=0., xLow=0., deltaX=0;\r
342 Int_t binMax=0, nBinsx=0;\r
616b80bc 343 Float_t mean[6]={0.,0.,0.,0.,0.,0.}, sigma[6]={0.,0.,0.,0.,0.,0.};\r
344 TF1 *fitfun[6]={0x0,0x0,0x0,0x0,0x0,0x0};\r
45b889f8 345 for(int k=0; k<6; k++){\r
76493b40 346 if(hTDC[k]->GetEntries()!=0){\r
347 binMax = hTDC[k]->GetMaximumBin();\r
348 if(binMax<=1){\r
349 printf("\n WARNING! Something wrong with det %d histo \n\n", k);\r
350 continue;\r
351 }\r
352 // \r
353 xUp = (hTDC[k]->GetXaxis())->GetXmax();\r
354 xLow = (hTDC[k]->GetXaxis())->GetXmin();\r
355 deltaX = xUp-xLow;\r
356 nBinsx = (hTDC[k]->GetXaxis())->GetNbins();\r
357 //printf(" xMax = %f\n", xLow+binMax*deltaX/nBinsx);\r
358 hTDC[k]->Fit("gaus","Q","",xLow+binMax*deltaX/nBinsx*0.75,xLow+binMax*deltaX/nBinsx*1.25);\r
359 fitfun[k] = hTDC[k]->GetFunction("gaus");\r
360 mean[k] = (Float_t) (fitfun[k]->GetParameter(1));\r
361 sigma[k] = (Float_t) (fitfun[k]->GetParameter(2));\r
ed8a15fd 362 //printf("\t Mean value from fit = %1.2f ns\n", mean[k]);\r
76493b40 363 //\r
364 fprintf(fileShuttle,"\t%f\t%f\n",mean[k], sigma[k]);\r
365 }\r
366 else printf(" WARNING! TDC histo %d has no entries and can't be processed!\n",k);\r
367 }\r
368 // \r
369 fclose(fileShuttle);\r
370 \r
371 TFile *histofile = new TFile(TDCHISTO_FILE,"RECREATE");\r
372 histofile->cd();\r
373 for(int k=0; k<6; k++) hTDC[k]->Write(); \r
374 histofile->Close();\r
375 //\r
376 for(Int_t j=0; j<6; j++) delete hTDC[j];\r
377 \r
378 /* store the result files on FES */\r
379 // [1] File with mapping\r
380 status = daqDA_FES_storeFile(MAPDATA_FILE, "MAPPING");\r
381 if(status){\r
382 printf("Failed to export file : %d\n",status);\r
383 return -1;\r
384 }\r
385 // [2] File with TDC data\r
386 status = daqDA_FES_storeFile(TDCDATA_FILE, "TDCDATA");\r
387 if(status){\r
388 printf("Failed to export TDC data file to DAQ FES\n");\r
389 return -1;\r
390 }\r
391 // [3] File with TDC histos\r
392 status = daqDA_FES_storeFile(TDCHISTO_FILE, "TDCHISTOS");\r
393 if(status){\r
394 printf("Failed to export TDC histos file to DAQ FES\n");\r
395 return -1;\r
396 }\r
397\r
398 return status;\r
399}\r