Changes suggested by Sylvain
[u/mrichter/AliRoot.git] / ZDC / ZDCPEDESTALda.cxx
1 /*
2
3 DAcase2.c
4
5 DAcase1.c
6
7 This program reads the DAQ data files passed as argument using the monitoring library.
8
9 It computes the average event size and populates local "./result.txt" file with the 
10 result.
11
12 The program reports about its processing progress.
13
14 Messages on stdout are exported to DAQ log system.
15
16 DA for ZDC standalone pedestal runs
17
18 Contact: Chiara.Oppedisano@to.infn.it
19 Link: /afs/cern.ch/user/c/chiarao/public/RawPed.date
20 Run Type: STANDALONE_PEDESTAL_RUN
21 DA Type: MON
22 Number of events needed: no constraint (tipically ~10^3)
23 Input Files: 
24 Output Files: ZDCPedestal.dat
25 Trigger Types Used: Standalone Trigger
26
27 */
28
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <Riostream.h>
32
33 // DATE
34 #include <event.h>
35 #include <monitor.h>
36 #include <daqDA.h>
37
38 //ROOT
39 #include <TRandom.h>
40 #include <TH1F.h>
41 #include <TH2F.h>
42 #include <TProfile.h>
43 #include <TF1.h>
44 #include <TFile.h>
45
46 //AliRoot
47 #include <AliRawReaderDate.h>
48 #include <AliZDCRawStream.h>
49
50
51 /* Main routine
52       Arguments: list of DATE raw data files
53 */
54 int main(int argc, char **argv) {
55
56   int status = 0;
57
58   /* log start of process */
59   printf("ZDC PEDESTAL monitoring program started\n");  
60
61   /* check that we got some arguments = list of files */
62   if (argc<2) {
63     printf("Wrong number of arguments\n");
64     return -1;
65   }
66
67   // --- Histograms for ADC pedestals 
68   //     [22 signal channels + 2 reference PTMs]  x 2 gain chains
69   //
70   TH1F::AddDirectory(0);
71   int const kNChannels = 24;
72   TH1F *hPedhg[kNChannels], *hPedOutOfTimehg[kNChannels];
73   TH2F *hPedCorrhg[kNChannels];
74   TH1F *hPedlg[kNChannels], *hPedOutOfTimelg[kNChannels];
75   TH2F *hPedCorrlg[kNChannels];
76   //
77   char namhist1hg[50], namhist2hg[50], namhist3hg[50];
78   char namhist1lg[50], namhist2lg[50], namhist3lg[50];
79   for(Int_t j=0; j<kNChannels; j++){
80      if(j<5){ // ZN1
81        sprintf(namhist1hg,"PedZN1hg_%d",j);
82        sprintf(namhist2hg,"PedZN1hgOutOfTime_%d",j);
83        sprintf(namhist3hg,"PedCorrZN1hg_%d",j);
84        //
85        sprintf(namhist1lg,"PedZN1lg_%d",j);
86        sprintf(namhist2lg,"PedZN1lgOutOfTime_%d",j);
87        sprintf(namhist3lg,"PedCorrZN1lg_%d",j);
88      }
89      else if(j>=5 && j<10){ // ZP1
90        sprintf(namhist1hg,"PedZP1hg_%d",j-5);
91        sprintf(namhist2hg,"PedZP1hgOutOfTime_%d",j-5);
92        sprintf(namhist3hg,"PedCorrZP1hg_%d",j-5);
93        //
94        sprintf(namhist1lg,"PedZP1lg_%d",j-5);
95        sprintf(namhist2lg,"PedZP1lgOutOfTime_%d",j-5);
96        sprintf(namhist3lg,"PedCorrZP1lg_%d",j-5);       
97      }
98      else if(j>=10 && j<12){ // ZEM
99        sprintf(namhist1hg,"PedZEMhg_%d",j-10);
100        sprintf(namhist2hg,"PedZEMhgOutOfTime_%d",j-10);
101        sprintf(namhist3hg,"PedCorrZEMhg_%d",j-10);
102        //
103        sprintf(namhist1lg,"PedZEMlg_%d",j-10);
104        sprintf(namhist2lg,"PedZEMlgOutOfTime_%d",j-10);
105        sprintf(namhist3lg,"PedCorrZEMlg_%d",j-10);
106      }
107      else if(j>=12 && j<17){ // ZN2
108        sprintf(namhist1hg,"PedZN2hg_%d",j-12);
109        sprintf(namhist2hg,"PedZN2hgOutOfTime_%d",j-12);
110        sprintf(namhist3hg,"PedCorrZN2hg_%d",j-12);
111        //
112        sprintf(namhist1lg,"PedZN2lg_%d",j-12);
113        sprintf(namhist2lg,"PedZN2lgOutOfTime_%d",j-12);
114        sprintf(namhist3lg,"PedCorrZN2lg_%d",j-12);
115      }
116      else if(j>=17 && j<22){ // ZP2
117        sprintf(namhist1hg,"PedZP2hg_%d",j-17);
118        sprintf(namhist2hg,"PedZP2hgOutOfTime_%d",j-17);
119        sprintf(namhist3hg,"PedCorrZP2hg_%d",j-17);
120        //
121        sprintf(namhist1lg,"PedZP2lg_%d",j-17);
122        sprintf(namhist2lg,"PedZP2lgOutOfTime_%d",j-17);
123        sprintf(namhist3lg,"PedCorrZP2lg_%d",j-17);
124      }
125      else if(j>=22 && j<24){ //Reference PMs
126        sprintf(namhist1hg,"PedRefhg_%d",j-22);
127        sprintf(namhist2hg,"PedRefhgOutOfTime_%d",j-22);
128        sprintf(namhist3hg,"PedCorrRefhg_%d",j-22);
129        //
130        sprintf(namhist1lg,"PedReflg_%d",j-22);
131        sprintf(namhist2lg,"PedReflgOutOfTime_%d",j-22);
132        sprintf(namhist3lg,"PedCorrReflg_%d",j-22);
133      }
134      // --- High gain chain histos
135      hPedhg[j] = new TH1F(namhist1hg, namhist1hg, 100,0., 200.);
136      hPedOutOfTimehg[j] = new TH1F(namhist2hg, namhist2hg, 100,0., 200.);
137      hPedCorrhg[j] = new TH2F(namhist3hg,namhist3hg,100,0.,200.,100,0.,200.);
138      // --- Low gain chain histos
139      hPedlg[j] = new TH1F(namhist1lg, namhist1lg, 100,0., 600.);
140      hPedOutOfTimelg[j] = new TH1F(namhist2lg, namhist2lg, 100,0., 600.);
141      hPedCorrlg[j] = new TH2F(namhist3lg,namhist3lg,100,0.,600.,100,0.,600.);
142   }
143
144
145   /* open result file */
146   FILE *fp=NULL;
147   fp=fopen("./result.txt","a");
148   if (fp==NULL) {
149     printf("Failed to open file\n");
150     return -1;
151   }
152   
153
154   /* report progress */
155   daqDA_progressReport(10);
156
157
158   /* init some counters */
159   int nevents_physics=0;
160   int nevents_total=0;
161
162   /* read the data files */
163   int n;
164   for (n=1;n<argc;n++) {
165    
166     status=monitorSetDataSource( argv[n] );
167     if (status!=0) {
168       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
169       return -1;
170     }
171
172     /* report progress */
173     /* in this example, indexed on the number of files */
174     daqDA_progressReport(10+80*n/argc);
175
176     /* read the file */
177     for(;;) {
178       struct eventHeaderStruct *event;
179       eventTypeType eventT;
180
181       /* get next event */
182       status=monitorGetEventDynamic((void **)&event);
183       if (status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
184       if (status!=0) {
185         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
186         return -1;
187       }
188
189       /* retry if got no event */
190       if (event==NULL) {
191         break;
192       }
193
194
195       /* use event - here, just write event id to result file */
196       eventT=event->eventType;
197     
198       if(eventT==PHYSICS_EVENT){
199         fprintf(fp,"Run #%lu, event size: %lu, BC:%u, Orbit:%u, Period:%u\n",
200           (unsigned long)event->eventRunNb,
201           (unsigned long)event->eventSize,
202           EVENT_ID_GET_BUNCH_CROSSING(event->eventId),
203           EVENT_ID_GET_ORBIT(event->eventId),
204           EVENT_ID_GET_PERIOD(event->eventId));
205         //
206         // Initalize raw-data reading and decoding
207         AliRawReader *reader = new AliRawReaderDate((void*)event);
208         const AliRawDataHeader* header = reader->GetDataHeader();
209         if(header) {
210          UChar_t message = header->GetAttributes();
211          if(message & 0x20){ // DEDICATED PEDESTAL RUN
212             printf("\t STANDALONE_PEDESTAL_RUN raw data found\n");
213             continue;
214          }
215          else{
216             printf("\t NO STANDALONE_PEDESTAL_RUN raw data found\n");
217             return -1;
218          }
219         }
220         //Commented until we won't have true Raw Data Header...
221         //else{
222         //   printf("\t ATTENTION! No Raw Data Header found!!!\n");
223           // return -1;
224         //}
225         //
226         AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
227         //
228         if (!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
229         Int_t counter=0;
230         Int_t RawADChg[kNChannels], RawADCoothg[kNChannels];
231         Int_t RawADClg[kNChannels], RawADCootlg[kNChannels];
232         for(Int_t j=0; j<kNChannels; j++){
233            RawADChg[j]=0; RawADCoothg[j]=0;
234            RawADClg[j]=0; RawADCootlg[j]=0;
235         }
236         while(rawStreamZDC->Next()){
237           Int_t index=-1;
238           if(rawStreamZDC->IsADCDataWord()){
239           if(rawStreamZDC->GetSector(1)!=5){ // Physics signals
240             if(rawStreamZDC->GetSector(0)==1) index = rawStreamZDC->GetSector(1); // *** ZN1
241             else if(rawStreamZDC->GetSector(0)==2) index = rawStreamZDC->GetSector(1)+5; // *** ZP1
242             else if(rawStreamZDC->GetSector(0)==3) index = rawStreamZDC->GetSector(1)+9; // *** ZEM
243             else if(rawStreamZDC->GetSector(0)==4) index = rawStreamZDC->GetSector(1)+12; // *** ZN2
244             else if(rawStreamZDC->GetSector(0)==5) index = rawStreamZDC->GetSector(1)+17; // *** ZP2
245           }
246           else{ // Reference PMs
247             index = (rawStreamZDC->GetSector(0)-1)/3+22;
248           }
249           //
250           /*printf("\t counter %d index %d det %d quad %d res %d ADC %d\n", counter, index,
251                 rawStreamZDC->GetSector(0), rawStreamZDC->GetSector(1), 
252                 rawStreamZDC->GetADCGain(), rawStreamZDC->GetADCValue());
253           */
254           //
255           if(counter<2*kNChannels){ // --- In-time pedestals (1st 48 raw data)
256             if(rawStreamZDC->GetADCGain()==0){ 
257               hPedhg[index]->Fill(rawStreamZDC->GetADCValue()); 
258               RawADChg[index] = rawStreamZDC->GetADCValue();
259             }
260             else{
261               hPedlg[index]->Fill(rawStreamZDC->GetADCValue()); 
262               RawADClg[index] = rawStreamZDC->GetADCValue();
263             }
264             }
265             else{  // --- Out-of-time pedestals
266               if(rawStreamZDC->GetADCGain()==0){
267                 hPedOutOfTimehg[index]->Fill(rawStreamZDC->GetADCValue());
268                 RawADCoothg[index] = rawStreamZDC->GetADCValue();
269               }
270               else{
271                 hPedOutOfTimelg[index]->Fill(rawStreamZDC->GetADCValue());
272                 RawADCootlg[index] = rawStreamZDC->GetADCValue();
273               }
274             }
275             counter++;
276            }//IsADCDataWord()
277            //
278            if(counter == 4*kNChannels){ // Last ADC channel -> Filling correlation histos
279              for(Int_t k=0; k<kNChannels; k++){
280               hPedCorrhg[k]->Fill(RawADCoothg[k], RawADChg[k]);
281               hPedCorrlg[k]->Fill(RawADCootlg[k], RawADClg[k]);
282              }
283            }
284          }
285          //
286          nevents_physics++;
287          //
288          delete reader;
289          delete rawStreamZDC;
290
291       }//(if PHYSICS_EVENT) 
292       nevents_total++;
293
294       /* free resources */
295       free(event);
296     
297     }
298   }  
299   
300   /* Analysis of the histograms */
301   //
302   FILE *fileShuttle;
303   const char *fName = "ZDCPedestal.dat";
304   fileShuttle = fopen(fName,"w");
305   //
306   Float_t MeanPed[2*kNChannels], MeanPedWidth[2*kNChannels], 
307         MeanPedOOT[2*kNChannels], MeanPedWidthOOT[2*kNChannels],
308         CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
309   // --- In-time pedestals
310   TF1 *ADCfunchg[kNChannels];
311   for(Int_t i=0; i<kNChannels; i++){
312      hPedhg[i]->Fit("gaus","Q");
313      ADCfunchg[i] = hPedhg[i]->GetFunction("gaus");
314      MeanPed[i] = ADCfunchg[i]->GetParameter(1);
315      MeanPedWidth[i] = ADCfunchg[i]->GetParameter(2);
316      fprintf(fileShuttle,"\t%f\t%f\n",MeanPed[i],MeanPedWidth[i]);
317      //printf("\t MeanPed[%d] = %f\n",i, MeanPed[i]);
318   }
319   TF1 *ADCfunclg[kNChannels];
320   for(Int_t i=0; i<kNChannels; i++){
321      hPedlg[i]->Fit("gaus","Q");
322      ADCfunclg[i] = hPedlg[i]->GetFunction("gaus");
323      MeanPed[i+kNChannels] = ADCfunclg[i]->GetParameter(1);
324      MeanPedWidth[i+kNChannels] = ADCfunclg[i]->GetParameter(2);
325      fprintf(fileShuttle,"\t%f\t%f\n",MeanPed[i+kNChannels],MeanPedWidth[i+kNChannels]);
326      //printf("\t MeanPed[%d] = %f\n",i+kNChannels, MeanPed[i+kNChannels]);
327   }
328   // --- Out-of-time pedestals
329   TF1 *ADCootfunchg[kNChannels];
330   for(Int_t i=0; i<kNChannels; i++){
331      hPedOutOfTimehg[i]->Fit("gaus","Q");
332      ADCootfunchg[i] = hPedOutOfTimehg[i]->GetFunction("gaus");
333      MeanPedOOT[i] = ADCootfunchg[i]->GetParameter(1);
334      MeanPedWidthOOT[i] = ADCootfunchg[i]->GetParameter(2);
335      fprintf(fileShuttle,"\t%f\t%f\n",MeanPedOOT[i],MeanPedWidthOOT[i]);
336      //printf("\t MeanPedOOT[%d] = %f\n",i, MeanPedOOT[i]);
337   }
338   TF1 *ADCootfunclg[kNChannels];
339   for(Int_t i=0; i<kNChannels; i++){
340      hPedOutOfTimelg[i]->Fit("gaus","Q");
341      ADCootfunclg[i] = hPedOutOfTimelg[i]->GetFunction("gaus");
342      MeanPedOOT[i+kNChannels] = ADCootfunclg[i]->GetParameter(1);
343      MeanPedWidthOOT[i+kNChannels] = ADCootfunclg[i]->GetParameter(2);
344      fprintf(fileShuttle,"\t%f\t%f\n",MeanPedOOT[i+kNChannels],MeanPedWidthOOT[i+kNChannels]);
345      //printf("\t MeanPedOOT[%d] = %f\n",i+kNChannels, MeanPedOOT[i+kNChannels]);
346   }
347   //
348   // --- Correlations
349   TProfile *hPedCorrProfhg[kNChannels], *hPedCorrProflg[kNChannels];
350   TF1 *ffunchg[kNChannels], *ffunclg[kNChannels];
351   char namhist4[50];
352   for(int i=0;i<kNChannels;i++) {
353      sprintf(namhist4,"ADCHRvsOOT%d_Prof",i);
354      hPedCorrProfhg[i] = hPedCorrhg[i]->ProfileX(namhist4,-1,-1,"S");
355      hPedCorrProfhg[i]->SetName(namhist4);
356      hPedCorrProfhg[i]->Fit("pol1","Q");
357      ffunchg[i] = hPedCorrProfhg[i]->GetFunction("pol1");
358      CorrCoeff0[i] = ffunchg[i]->GetParameter(0);
359      CorrCoeff1[i] = ffunchg[i]->GetParameter(1);
360      fprintf(fileShuttle,"\t%f\t%f\n",CorrCoeff0[i],CorrCoeff1[i]);
361      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",i, CorrCoeff0[i], i, CorrCoeff1[i]);
362   }    
363   for(int i=0;i<kNChannels;i++) {
364      sprintf(namhist4,"ADCLRvsOOT%d_Prof",i);
365      hPedCorrProflg[i] = hPedCorrlg[i]->ProfileX(namhist4,-1,-1,"S");
366      hPedCorrProflg[i]->SetName(namhist4);
367      hPedCorrProflg[i]->Fit("pol1","Q");
368      ffunclg[i] = hPedCorrProflg[i]->GetFunction("pol1");
369      CorrCoeff0[i+kNChannels] = ffunclg[i]->GetParameter(0);
370      CorrCoeff1[i+kNChannels] = ffunclg[i]->GetParameter(1);
371      fprintf(fileShuttle,"\t%f\t%f\n",CorrCoeff0[i+kNChannels],CorrCoeff1[i+kNChannels]);
372      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",
373      //         i+kNChannels, CorrCoeff0[i+kNChannels], i+kNChannels, CorrCoeff1[i+kNChannels]);
374   }    
375   //                                                   
376   fclose(fileShuttle);
377   //
378   for(Int_t j=0; j<kNChannels; j++){
379      delete hPedhg[j];
380      delete hPedOutOfTimehg[j];
381      delete hPedCorrhg[j];
382      delete hPedlg[j];
383      delete hPedOutOfTimelg[j];
384      delete hPedCorrlg[j];
385   }
386
387   /* write report */
388   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
389
390   /* close result file */
391   fclose(fp);
392   
393   /* report progress */
394   daqDA_progressReport(90);
395
396   /* store the result file on FES */
397   status = daqDA_FES_storeFile(fName,"ZDCPEDESTAL_data");
398   if(status){
399     printf("Failed to export file : %d\n",status);
400     return -1;
401   }
402
403   /* report progress */
404   daqDA_progressReport(100);
405
406
407   return status;
408 }