e9605be085380cf291c7b9607411ac05b299f218
[u/mrichter/AliRoot.git] / ZDC / ZDCMAPPINGda.cxx
1 /*\r
2 \r
3 This program reads the DAQ data files passed as argument using the monitoring library.\r
4 \r
5 The program reports about its processing progress.\r
6 \r
7 Messages on stdout are exported to DAQ log system.\r
8 \r
9 DA to write mapping for ADC modules and VME scaler\r
10 \r
11 Contact: Chiara.Oppedisano@to.infn.it\r
12 Link: \r
13 Run Type: PHYSICS\r
14 DA Type: MON\r
15 Number of events needed: at least 50% \r
16 Input Files:  none\r
17 Output Files: ZDCChMapping.dat\r
18 Trigger 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
52 int main(int argc, char **argv) {\r
53 \r
54   // needed for streamer application\r
55   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",\r
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
66 \r
67   //TMinuitMinimizer m; \r
68   TVirtualFitter::SetDefaultFitter("Minuit");\r
69 \r
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
78   int const kNModules = 10;\r
79   int const kNChannels = 24;\r
80   int const kNScChannels = 32;  \r
81   int const kZDCTDCGeo = 4;\r
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
114   TH1F * hTDC[6]={0x0,0x0,0x0,0x0,0x0,0x0};\r
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
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
129   if(argc<2) {\r
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
137   for(int n=1;n<argc;n++){\r
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
161 \r
162       struct eventHeaderStruct *event;\r
163       eventTypeType eventT;\r
164  \r
165       /* check shutdown condition */\r
166       if (daqDA_checkShutdown()) {break;}\r
167 \r
168       /* get next event */\r
169       status=monitorGetEventDynamic((void **)&event);\r
170       if(status==MON_ERR_EOF){\r
171         printf ("End of File detected\n");\r
172         break; /* end of monitoring file has been reached */\r
173       }\r
174       if(status!=0) {\r
175         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));\r
176         return -1;\r
177       }\r
178 \r
179       /* retry if got no event */\r
180       if(event==NULL) continue;\r
181       \r
182       // Initalize raw-data reading and decoding\r
183       AliRawReader *reader = new AliRawReaderDate((void*)event);\r
184       reader->Reset();\r
185       reader->Select("ZDC");\r
186       // --- Reading event header\r
187       //UInt_t evtype = reader->GetType();\r
188       //printf("\t ZDCMAPPINGda -> run # %d\n",reader->GetRunNumber());\r
189       //\r
190       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);\r
191         \r
192         \r
193       /* use event - here, just write event id to result file */\r
194       eventT=event->eventType;\r
195       \r
196       if(eventT==START_OF_DATA){\r
197                 \r
198         rawStreamZDC->SetSODReading(kTRUE);\r
199         \r
200         // --------------------------------------------------------\r
201         // --- Writing ascii data file for the Shuttle preprocessor\r
202         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");\r
203         //if(rawStreamZDC->Next()){\r
204           while((rawStreamZDC->Next())){\r
205             if(rawStreamZDC->IsHeaderMapping()){ // mapping header\r
206                iMod++;\r
207                modGeo[iMod]  = rawStreamZDC->GetADCModule();\r
208                modType[iMod] = rawStreamZDC->GetModType();\r
209                modNCh[iMod]  = rawStreamZDC->GetADCNChannels();\r
210             }\r
211             if(rawStreamZDC->IsChMapping()){ \r
212               if(modType[iMod]==1){ // ADC mapping ----------------------\r
213                 adcMod[ich]  = rawStreamZDC->GetADCModFromMap(ich);\r
214                 adcCh[ich]   = rawStreamZDC->GetADCChFromMap(ich);\r
215                 sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);\r
216                 det[ich]     = rawStreamZDC->GetDetectorFromMap(ich);\r
217                 sec[ich]     = rawStreamZDC->GetTowerFromMap(ich);\r
218                 ich++;\r
219               }\r
220               else if(modType[iMod]==2){ // VME scaler mapping -------------\r
221                 scMod[iScCh]     = rawStreamZDC->GetScalerModFromMap(iScCh);\r
222                 scCh[iScCh]      = rawStreamZDC->GetScalerChFromMap(iScCh);\r
223                 scSigCode[iScCh] = rawStreamZDC->GetScalerSignFromMap(iScCh);\r
224                 scDet[iScCh]     = rawStreamZDC->GetScDetectorFromMap(iScCh);\r
225                 scSec[iScCh]     = rawStreamZDC->GetScTowerFromMap(iScCh);\r
226                 iScCh++;\r
227               }\r
228               else if(modType[iMod]==6 && modGeo[iMod]==4){ // ZDC TDC mapping --------------------\r
229                 tdcMod[itdcCh]     = rawStreamZDC->GetTDCModFromMap(itdcCh);\r
230                 tdcCh[itdcCh]      = rawStreamZDC->GetTDCChFromMap(itdcCh);\r
231                 tdcSigCode[itdcCh] = rawStreamZDC->GetTDCSignFromMap(itdcCh);\r
232                 itdcCh++;\r
233               }\r
234             }\r
235           }\r
236           // Writing data on output FXS file\r
237           for(Int_t is=0; is<2*kNChannels; is++){\r
238              fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",\r
239                is,adcMod[is],adcCh[is],sigCode[is],det[is],sec[is]);\r
240              //printf("  Mapping DA -> %d ADC: mod %d ch %d, code %d det %d, sec %d\n",\r
241              //  is,adcMod[is],adcCh[is],sigCode[is],det[is],sec[is]);\r
242           }\r
243           for(Int_t is=0; is<kNScChannels; is++){\r
244              fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",\r
245                is,scMod[is],scCh[is],scSigCode[is],scDet[is],scSec[is]);\r
246              //if(scMod[is]!=-1) printf("  Mapping DA -> %d Scaler: mod %d ch %d, code %d det %d, sec %d\n",\r
247              //  is,scMod[is],scCh[is],scSigCode[is],scDet[is],scSec[is]);\r
248           }\r
249           for(Int_t is=0; is<kNScChannels; is++){\r
250              fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\n",\r
251                is,tdcMod[is],tdcCh[is],tdcSigCode[is]);\r
252              //if(tdcMod[is]!=-1) printf("  Mapping DA -> %d TDC: mod %d ch %d, code %d\n",\r
253              //  is,tdcMod[is],tdcCh[is],tdcSigCode[is]);\r
254           }\r
255           for(Int_t is=0; is<kNModules; is++){\r
256              fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\n",\r
257              modGeo[is],modType[is],modNCh[is]);\r
258              //printf("  Mapping DA -> Module mapping: geo %d type %d #ch %d\n",\r
259              //  modGeo[is],modType[is],modNCh[is]);\r
260           }\r
261           \r
262         //} //if (rawstream)\r
263         fclose(mapFile4Shuttle);\r
264       }// SOD event\r
265       else if(eventT==PHYSICS_EVENT){ \r
266 \r
267         for(int ij=0; ij<6; ij++) tdcData[ij]=-999.;\r
268         rawStreamZDC->SetSODReading(kTRUE);\r
269 \r
270         // ----- Setting ch. mapping -----\r
271         for(int jk=0; jk<2*kNChannels; jk++){\r
272           rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);\r
273           rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);\r
274           rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);\r
275           rawStreamZDC->SetMapDet(jk, det[jk]);\r
276           rawStreamZDC->SetMapTow(jk, sec[jk]);\r
277         }\r
278         \r
279         while(rawStreamZDC->Next()){\r
280          if(rawStreamZDC->GetADCModule()!=kZDCTDCGeo) continue; //skipping ADCs, scalers and trigger cards\r
281 \r
282           if(rawStreamZDC->GetADCModule()==kZDCTDCGeo && rawStreamZDC->IsZDCTDCDatum()==kTRUE){\r
283              //\r
284              itdc = rawStreamZDC->GetChannel(); \r
285              if(itdc==iprevtdc) ihittdc++;\r
286              else ihittdc=0;\r
287              iprevtdc=itdc;\r
288              //\r
289              if(((itdc>=8 && itdc<=13) || itdc==15) && ihittdc==0){\r
290                if(itdc!=15){\r
291                  tdcData[itdc-8] = 0.025*rawStreamZDC->GetZDCTDCDatum();\r
292                  //\r
293                  //printf("   ev.%d **** TDC%d %1.0f ns  \n",nphys,itdc, tdcData[itdc-8]);\r
294                }\r
295                //\r
296                else if(itdc==15){\r
297                   tdcL0 = 0.025*rawStreamZDC->GetZDCTDCDatum();\r
298                   //\r
299                   //printf("   ev.%d ++++ TDCL0 %1.0f ns  \n",nphys,tdcL0);\r
300                   //\r
301                   for(int ic=0; ic<6; ic++){\r
302                     if(tdcData[ic]!=-999. && tdcL0!=-999.){\r
303                       hTDC[ic]->Fill(tdcData[ic]-tdcL0);\r
304                       //printf(" ev.%d -> Filling histo%d: %f ns\n",nphys,ic, tdcData[ic]-tdcL0);\r
305                     }\r
306                   }\r
307                }\r
308              }\r
309           }\r
310         }\r
311         \r
312         nphys++;\r
313       \r
314         delete rawStreamZDC;\r
315         rawStreamZDC = 0x0;     \r
316 \r
317       }//(if PHYSICS_EVENT) \r
318       else if(eventT==END_OF_RUN){\r
319         printf("End Of Run detected\n");\r
320         break;\r
321       }\r
322       \r
323       iev++; \r
324 \r
325       /* free resources */\r
326       free(event);\r
327  \r
328     } // event loop    \r
329       \r
330   }\r
331   printf(" \n # of physics events analyzed = %d\n\n", nphys);\r
332 \r
333   /* Analysis of the histograms */\r
334   //\r
335   FILE *fileShuttle;\r
336   fileShuttle = fopen(TDCDATA_FILE,"w");\r
337   //\r
338   Float_t xUp=0., xLow=0., deltaX=0;\r
339   Int_t binMax=0, nBinsx=0;\r
340   Float_t mean[6]={0.,0.,0.,0.,0.,0.}, sigma[6]={0.,0.,0.,0.,0.,0.};\r
341   TF1 *fitfun[6]={0x0,0x0,0x0,0x0,0x0,0x0};\r
342   for(int k=0; k<6; k++){\r
343     if(hTDC[k]->GetEntries()!=0){\r
344        binMax = hTDC[k]->GetMaximumBin();\r
345        if(binMax<=1){\r
346          printf("\n WARNING! Something wrong with det %d histo \n\n", k);\r
347          continue;\r
348        }\r
349        // \r
350        xUp = (hTDC[k]->GetXaxis())->GetXmax();\r
351        xLow = (hTDC[k]->GetXaxis())->GetXmin();\r
352        deltaX = xUp-xLow;\r
353        nBinsx = (hTDC[k]->GetXaxis())->GetNbins();\r
354        //printf(" xMax = %f\n", xLow+binMax*deltaX/nBinsx);\r
355        hTDC[k]->Fit("gaus","Q","",xLow+binMax*deltaX/nBinsx*0.75,xLow+binMax*deltaX/nBinsx*1.25);\r
356        fitfun[k] = hTDC[k]->GetFunction("gaus");\r
357        mean[k] = (Float_t) (fitfun[k]->GetParameter(1));\r
358        sigma[k] = (Float_t) (fitfun[k]->GetParameter(2));\r
359        //printf("\t Mean value from fit = %1.2f\n", mean[k]);\r
360        //\r
361        fprintf(fileShuttle,"\t%f\t%f\n",mean[k], sigma[k]);\r
362      }\r
363      else printf("  WARNING! TDC histo %d has no entries and can't be processed!\n",k);\r
364   }\r
365   //                                                   \r
366   fclose(fileShuttle);\r
367   \r
368   TFile *histofile = new TFile(TDCHISTO_FILE,"RECREATE");\r
369   histofile->cd();\r
370   for(int k=0; k<6; k++) hTDC[k]->Write();                                                     \r
371   histofile->Close();\r
372   //\r
373   for(Int_t j=0; j<6; j++) delete hTDC[j];\r
374   \r
375   /* store the result files on FES */\r
376   // [1] File with mapping\r
377   status = daqDA_FES_storeFile(MAPDATA_FILE, "MAPPING");\r
378   if(status){\r
379     printf("Failed to export file : %d\n",status);\r
380     return -1;\r
381   }\r
382   // [2] File with TDC data\r
383   status = daqDA_FES_storeFile(TDCDATA_FILE, "TDCDATA");\r
384   if(status){\r
385     printf("Failed to export TDC data file to DAQ FES\n");\r
386     return -1;\r
387   }\r
388   // [3] File with TDC histos\r
389   status = daqDA_FES_storeFile(TDCHISTO_FILE, "TDCHISTOS");\r
390   if(status){\r
391     printf("Failed to export TDC histos file to DAQ FES\n");\r
392     return -1;\r
393   }\r
394 \r
395   return status;\r
396 }\r