88dff12d3e4b0aba1dd0a2ef0e5bb1cfa9642e3f
[u/mrichter/AliRoot.git] / ZDC / ZDCLASERda.cxx
1 /*
2
3 This program reads the DAQ data files passed as argument using the monitoring library.
4
5 It computes the average event size and populates local "./result.txt" file with the 
6 result.
7
8 The program reports about its processing progress.
9
10 Messages on stdout are exported to DAQ log system.
11
12 DA for ZDC standalone pedestal runs
13
14 Contact: Chiara.Oppedisano@to.infn.it
15 Link: 
16 Run Type: STANDALONE_LASER_RUN
17 DA Type: LDC
18 Number of events needed: no constraint (tipically ~10^3)
19 Input Files: 
20 Output Files: ZDCLaser.dat
21 Trigger Types Used: Standalone Trigger
22
23 */
24 #define PEDDATA_FILE  "ZDCPedestal.dat"
25 #define MAPDATA_FILE  "ZDCChMapping.dat"
26 #define LASDATA_FILE  "ZDCLaserCalib.dat"
27
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <Riostream.h>
31
32 // DATE
33 #include <event.h>
34 #include <monitor.h>
35 #include <daqDA.h>
36
37 //ROOT
38 #include <TRandom.h>
39 #include <TH1F.h>
40 #include <TF1.h>
41 #include <TFile.h>
42 #include <TFitter.h>
43
44 //AliRoot
45 #include <AliRawReaderDate.h>
46 #include <AliRawEventHeaderBase.h>
47 #include <AliZDCRawStream.h>
48
49
50 /* Main routine
51       Arguments: list of DATE raw data files
52 */
53 int main(int argc, char **argv) {
54   
55   TFitter *minuitFit = new TFitter(4);
56   TVirtualFitter::SetFitter(minuitFit);
57
58   int status = 0;
59
60   /* log start of process */
61   printf("\nZDC LASER program started\n");  
62
63   /* check that we got some arguments = list of files */
64   if (argc<2) {
65     printf("Wrong number of arguments\n");
66     return -1;
67   }
68
69   // --- Histograms for LASER runs 
70   //     20 signal channels + 2 reference PTMs
71   //
72   TH1F::AddDirectory(0);
73   //
74   TH1F *hPMRefChg = new TH1F("hPMRefChg","hPMRefChg", 100,0.,1000.);
75   TH1F *hPMRefAhg = new TH1F("hPMRefAhg","hPMRefAhg", 100,0.,1000.);
76   //
77   TH1F *hPMRefClg = new TH1F("hPMRefClg","hPMRefClg", 100,0.,4000.);
78   TH1F *hPMRefAlg = new TH1F("hPMRefAlg","hPMRefAlg", 100,0.,4000.);
79
80
81   /* open result file */
82   FILE *fp=NULL;
83   fp=fopen("./result.txt","a");
84   if (fp==NULL) {
85     printf("Failed to open file\n");
86     return -1;
87   }
88   
89   FILE *mapFile4Shuttle;
90         
91   // *** To analyze LASER events you MUST have a pedestal data file!!!
92   // *** -> check if a pedestal run has been analyzied
93   int read = 0;
94   read = daqDA_DB_getFile(PEDDATA_FILE,PEDDATA_FILE);
95   if(read){
96     printf("\t ERROR!!! ZDCPedestal.dat file NOT FOUND in DAQ db!!!\n");
97     return -1;
98   }
99   else printf("\t ZDCPedestal.dat file retrieved from DAQ db\n");
100   
101   FILE *filePed = fopen(PEDDATA_FILE,"r");
102   if (filePed==NULL) {
103     printf("\t ERROR!!! Can't open ZDCPedestal.dat file!!!\n");
104     return -1;
105   }
106
107   // 144 = 48 in-time + 48 out-of-time + 48 correlations
108   Float_t readValues[2][144], MeanPed[44], MeanPedWidth[44], 
109         MeanPedOOT[44], MeanPedWidthOOT[44];
110   // ***************************************************
111   //   Unless we have a narrow correlation to fit we
112   //    don't fit and store in-time vs. out-of-time
113   //    histograms -> mean pedstal subtracted!!!!!!
114   // ***************************************************
115   //Float_t CorrCoeff0[44], CorrCoeff1[44];
116   //
117   for(int jj=0; jj<144; jj++){
118     for(int ii=0; ii<2; ii++){
119        fscanf(filePed,"%f",&readValues[ii][jj]);
120     }
121     if(jj<48){
122       MeanPed[jj] = readValues[0][jj];
123       MeanPedWidth[jj] = readValues[1][jj];
124       //printf("\t MeanPed[%d] = %1.1f\n",jj, MeanPed[jj]);
125     }
126     else if(jj>48 && jj<96){
127       MeanPedOOT[jj-48] = readValues[0][jj];
128       MeanPedWidthOOT[jj-48] = readValues[1][jj];
129     }
130     /*else if(jj>144){
131       CorrCoeff0[jj-96] = readValues[0][jj]; 
132       CorrCoeff1[jj-96] = readValues[1][jj];;
133     }
134     */
135   }
136
137   /* report progress */
138   daqDA_progressReport(10);
139
140
141   /* init some counters */
142   int nevents_physics=0;
143   int nevents_total=0;
144
145   /* read the data files */
146   int n;
147   for (n=1;n<argc;n++) {
148    
149     status=monitorSetDataSource( argv[n] );
150     if (status!=0) {
151       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
152       return -1;
153     }
154
155     /* report progress */
156     /* in this example, indexed on the number of files */
157     daqDA_progressReport(10+80*n/argc);
158
159     /* read the file */
160     for(;;) {
161       struct eventHeaderStruct *event;
162       eventTypeType eventT;
163
164       /* get next event */
165       status=monitorGetEventDynamic((void **)&event);
166       if (status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
167       if (status!=0) {
168         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
169         return -1;
170       }
171
172       /* retry if got no event */
173       if (event==NULL) {
174         break;
175       }
176
177       // Initalize raw-data reading and decoding
178       AliRawReader *reader = new AliRawReaderDate((void*)event);
179       reader->Select("ZDC");
180       // --- Reading event header
181       //UInt_t evtype = reader->GetType();
182       //printf("\n\t ZDCLASERda -> ev. type %d\n",evtype);
183       //printf("\t ZDCLASERda -> run # %d\n",reader->GetRunNumber());
184       //
185       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
186         
187
188       /* use event - here, just write event id to result file */
189       eventT=event->eventType;
190       
191       Int_t ich=0, adcMod[48], adcCh[48], sigCode[48], det[48], sec[48];
192       if(eventT==START_OF_DATA){
193                 
194         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
195         else{
196           while(rawStreamZDC->Next()){
197             if(rawStreamZDC->IsChMapping()){
198               adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);
199               adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);
200               sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
201               det[ich] = rawStreamZDC->GetDetectorFromMap(ich);
202               sec[ich] = rawStreamZDC->GetTowerFromMap(ich);
203               ich++;
204             }
205           }
206         }
207         // --------------------------------------------------------
208         // --- Writing ascii data file for the Shuttle preprocessor
209         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
210         for(Int_t i=0; i<ich; i++){
211            fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",i,
212              adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
213            //
214            //printf("ZDCLASERDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
215            //      i,adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
216         }
217         fclose(mapFile4Shuttle);
218       }
219
220       /* use event - here, just write event id to result file */
221       eventT=event->eventType;
222     
223       if(eventT==PHYSICS_EVENT){
224         //
225         // --- Reading data header
226         reader->ReadHeader();
227         const AliRawDataHeader* header = reader->GetDataHeader();
228         if(header) {
229          UChar_t message = header->GetAttributes();
230          if(message & 0x20){ // DEDICATED LASER RUN
231             //printf("\t STANDALONE_LASER_RUN raw data found\n");
232             continue;
233          }
234          else{
235             printf("\t NO STANDALONE_LASER_RUN raw data found\n");
236             return -1;
237          }
238         }
239         else{
240            printf("\t ATTENTION! No Raw Data Header found!!!\n");
241            return -1;
242         }
243
244         if (!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
245         //
246         // ----- Setting ch. mapping -----
247         for(Int_t jk=0; jk<48; jk++){
248           rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
249           rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
250           rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
251           rawStreamZDC->SetMapDet(jk, det[jk]);
252           rawStreamZDC->SetMapTow(jk, sec[jk]);
253         }
254         //
255         while(rawStreamZDC->Next()){
256           Int_t index=-1;
257           // Getting data only for reference PMTs (sector[1]=5)
258           if((rawStreamZDC->IsADCDataWord()) && (rawStreamZDC->GetSector(1)==5)){
259             index = rawStreamZDC->GetADCChannel();
260             Float_t Pedestal = MeanPed[index];
261             Float_t CorrADC = rawStreamZDC->GetADCValue() - Pedestal;
262             
263             // ==== HIGH GAIN CHAIN
264             if(rawStreamZDC->GetADCGain() == 0){
265               // %%%%% PMRef chain side C
266               if(rawStreamZDC->GetSector(0)==1) hPMRefChg->Fill(CorrADC);
267               // %%%%% PMRef side A
268               else if(rawStreamZDC->GetSector(0)==4) hPMRefAhg->Fill(CorrADC);
269             }
270             // ==== LOW GAIN CHAIN
271             else{
272               // %%%%% PMRef chain side C
273               if(rawStreamZDC->GetSector(0)==1) hPMRefClg->Fill(CorrADC);
274               // %%%%% PMRef side A
275               else if(rawStreamZDC->GetSector(0)==4) hPMRefAlg->Fill(CorrADC);
276             }
277           }//IsADCDataWord()
278           //
279          }
280          //
281          nevents_physics++;
282          //
283          delete reader;
284          delete rawStreamZDC;
285
286       }//(if PHYSICS_EVENT) 
287       nevents_total++;
288
289       /* free resources */
290       free(event);
291     
292     }
293   }  
294   
295   /* Analysis of the histograms */
296   //
297   Int_t maxBinRef[4], nBinRef[4];
298   Float_t xMaxRef[4], maxXvalRef[4], xlowRef[4]; 
299   Float_t meanRef[2], sigmaRef[2];
300   TF1 *funRef[4];
301   
302   // ~~~~~~~~ PM Ref side C high gain chain ~~~~~~~~
303   maxBinRef[0] = hPMRefChg->GetMaximumBin();
304   nBinRef[0] = (hPMRefChg->GetXaxis())->GetNbins();
305   xMaxRef[0] = (hPMRefChg->GetXaxis())->GetXmax();
306   maxXvalRef[0] = maxBinRef[0]*xMaxRef[0]/nBinRef[0];
307   // 
308   if(maxXvalRef[0]-100.<0.) {xlowRef[0]=0.;}
309   else xlowRef[0] = maxXvalRef[0];
310   hPMRefChg->Fit("gaus","Q","",xlowRef[0],maxXvalRef[0]+100.);
311   funRef[0] = hPMRefChg->GetFunction("gaus");
312   meanRef[0] = (Float_t) (funRef[0]->GetParameter(1));
313   sigmaRef[0] = (Float_t) (funRef[0]->GetParameter(2));
314   
315   // ~~~~~~~~ PM Ref side A high gain chain ~~~~~~~~
316   maxBinRef[1] = hPMRefAhg->GetMaximumBin();
317   nBinRef[1] = (hPMRefAhg->GetXaxis())->GetNbins();
318   xMaxRef[1] = (hPMRefAhg->GetXaxis())->GetXmax();
319   maxXvalRef[1] = maxBinRef[1]*xMaxRef[1]/nBinRef[1];
320   //
321   if(maxXvalRef[1]-100.<0.) {xlowRef[1]=0.;}
322   else xlowRef[1] = maxXvalRef[1];
323   hPMRefAhg->Fit("gaus","Q","",xlowRef[1],maxXvalRef[1]+100.);
324   funRef[1] = hPMRefAhg->GetFunction("gaus");
325   meanRef[1] = (Float_t) (funRef[1]->GetParameter(1));
326   sigmaRef[1] = (Float_t) (funRef[1]->GetParameter(2));
327   
328   // ~~~~~~~~ PM Ref side C low gain chain ~~~~~~~~
329   maxBinRef[2] = hPMRefClg->GetMaximumBin();
330   nBinRef[2] = (hPMRefClg->GetXaxis())->GetNbins();
331   xMaxRef[2] = (hPMRefClg->GetXaxis())->GetXmax();
332   maxXvalRef[2] = maxBinRef[2]*xMaxRef[2]/nBinRef[2];
333   //
334   if(maxXvalRef[2]-100.<0.) {xlowRef[2]=0.;}
335   else xlowRef[2] = maxXvalRef[2];
336   hPMRefClg->Fit("gaus","Q","",xlowRef[2],maxXvalRef[2]+100.);
337   funRef[2] = hPMRefClg->GetFunction("gaus");
338   meanRef[2] = (Float_t) (funRef[2]->GetParameter(1));
339   sigmaRef[2] = (Float_t) (funRef[2]->GetParameter(2));
340   
341   // ~~~~~~~~ PM Ref side A low gain chain ~~~~~~~~
342   maxBinRef[3] = hPMRefAlg->GetMaximumBin();
343   nBinRef[3] = (hPMRefAlg->GetXaxis())->GetNbins();
344   xMaxRef[3] = (hPMRefAlg->GetXaxis())->GetXmax();
345   maxXvalRef[3] = maxBinRef[3]*xMaxRef[3]/nBinRef[3];
346   //
347   if(maxXvalRef[3]-100.<0.) {xlowRef[3]=0.;}
348   else xlowRef[3] = maxXvalRef[3];
349   hPMRefAlg->Fit("gaus","Q","",xlowRef[3],maxXvalRef[3]+100.);
350   funRef[3] = hPMRefAlg->GetFunction("gaus");
351   meanRef[3] = (Float_t) (funRef[3]->GetParameter(1));
352   sigmaRef[3] = (Float_t) (funRef[3]->GetParameter(2));
353   //
354   FILE *fileShuttle;
355   fileShuttle = fopen(LASDATA_FILE,"w");
356   for(Int_t i=0; i<4; i++)  fprintf(fileShuttle,"\t%f\t%f\n",meanRef[i], sigmaRef[i]); 
357   //                                                   
358   fclose(fileShuttle);
359   //
360   delete hPMRefChg;
361   delete hPMRefAhg;
362   delete hPMRefClg;
363   delete hPMRefAlg;
364
365   //delete minuitFit;
366   TVirtualFitter::SetFitter(0);
367
368   /* write report */
369   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
370
371   /* close result file */
372   fclose(fp);
373   
374   /* report progress */
375   daqDA_progressReport(90);
376
377   /* store the result file on FES */
378   status = daqDA_FES_storeFile(MAPDATA_FILE,MAPDATA_FILE);
379   if(status){
380     printf("Failed to export file : %d\n",status);
381     return -1;
382   }
383   //
384   status = daqDA_FES_storeFile(LASDATA_FILE,LASDATA_FILE);
385   if(status){
386     printf("Failed to export file : %d\n",status);
387     return -1;
388   }
389
390   /* report progress */
391   daqDA_progressReport(100);
392
393
394   return status;
395 }