DAs upgraded to AliZDCRawStream class
[u/mrichter/AliRoot.git] / ZDC / ZDCPEDESTALda.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_PEDESTAL_RUN
17 DA Type: LDC
18 Number of events needed: no constraint (tipically ~10^3)
19 Input Files:  
20 Output Files: ZDCPedestal.dat, ZDCChMapping.dat
21 Trigger Types Used: Standalone Trigger
22
23 */
24 #define PEDDATA_FILE  "ZDCPedestal.dat"
25 #define MAPDATA_FILE  "ZDCChMapping.dat"
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <Riostream.h>
30
31 // DATE
32 #include <event.h>
33 #include <monitor.h>
34 #include <daqDA.h>
35
36 //ROOT
37 #include <TRandom.h>
38 #include <TH1F.h>
39 #include <TH2F.h>
40 #include <TProfile.h>
41 #include <TF1.h>
42 #include <TFile.h>
43 #include <TFitter.h>
44
45 //AliRoot
46 #include <AliRawReaderDate.h>
47 #include <AliRawEventHeaderBase.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   TFitter *minuitFit = new TFitter(4);
57   TVirtualFitter::SetFitter(minuitFit);
58
59   int status = 0;
60
61   /* log start of process */
62   printf("\n ZDC PEDESTAL program started\n");  
63
64   /* check that we got some arguments = list of files */
65   if (argc<2) {
66     printf("Wrong number of arguments\n");
67     return -1;
68   }
69
70   // --- Histograms for ADC pedestals 
71   //     [22 signal channels + 2 reference PTMs]  x 2 gain chains
72   //
73   TH1F::AddDirectory(0);
74   int const kNChannels = 24;
75   TH1F *hPedhg[kNChannels], *hPedOutOfTimehg[kNChannels];
76   TH2F *hPedCorrhg[kNChannels];
77   TH1F *hPedlg[kNChannels], *hPedOutOfTimelg[kNChannels];
78   TH2F *hPedCorrlg[kNChannels];
79   //
80   char namhist1hg[50], namhist2hg[50], namhist3hg[50];
81   char namhist1lg[50], namhist2lg[50], namhist3lg[50];
82   for(Int_t j=0; j<kNChannels; j++){
83      if(j<=4){ // ZNC
84        sprintf(namhist1hg,"PedZNChg_%d",j);
85        sprintf(namhist2hg,"PedZNChgOutOfTime_%d",j);
86        sprintf(namhist3hg,"PedCorrZNChg_%d",j);
87        //
88        sprintf(namhist1lg,"PedZNClg_%d",j);
89        sprintf(namhist2lg,"PedZNClgOutOfTime_%d",j);
90        sprintf(namhist3lg,"PedCorrZNClg_%d",j);
91      }
92      else if(j>=5 && j<=9){ // ZPC
93        sprintf(namhist1hg,"PedZPChg_%d",j-5);
94        sprintf(namhist2hg,"PedZPChgOutOfTime_%d",j-5);
95        sprintf(namhist3hg,"PedCorrZPChg_%d",j-5);
96        //
97        sprintf(namhist1lg,"PedZPClg_%d",j-5);
98        sprintf(namhist2lg,"PedZPClgOutOfTime_%d",j-5);
99        sprintf(namhist3lg,"PedCorrZPClg_%d",j-5);       
100      }
101      else if(j==10 || j==11){ // ZEM
102        sprintf(namhist1hg,"PedZEMhg_%d",j-10);
103        sprintf(namhist2hg,"PedZEMhgOutOfTime_%d",j-10);
104        sprintf(namhist3hg,"PedCorrZEMhg_%d",j-10);
105        //
106        sprintf(namhist1lg,"PedZEMlg_%d",j-10);
107        sprintf(namhist2lg,"PedZEMlgOutOfTime_%d",j-10);
108        sprintf(namhist3lg,"PedCorrZEMlg_%d",j-10);
109      }
110      else if(j>=12 && j<=16){ // ZNA
111        sprintf(namhist1hg,"PedZNAhg_%d",j-12);
112        sprintf(namhist2hg,"PedZNAhgOutOfTime_%d",j-12);
113        sprintf(namhist3hg,"PedCorrZNAhg_%d",j-12);
114        //
115        sprintf(namhist1lg,"PedZNAlg_%d",j-12);
116        sprintf(namhist2lg,"PedZNAlgOutOfTime_%d",j-12);
117        sprintf(namhist3lg,"PedCorrZNAlg_%d",j-12);
118      }
119      else if(j>=17 && j<=21){ // ZPA
120        sprintf(namhist1hg,"PedZPAhg_%d",j-17);
121        sprintf(namhist2hg,"PedZPAhgOutOfTime_%d",j-17);
122        sprintf(namhist3hg,"PedCorrZPAhg_%d",j-17);
123        //
124        sprintf(namhist1lg,"PedZPAlg_%d",j-17);
125        sprintf(namhist2lg,"PedZPAlgOutOfTime_%d",j-17);
126        sprintf(namhist3lg,"PedCorrZPAlg_%d",j-17);
127      }
128      else if(j==22 || j==24){ //Reference PMs
129        sprintf(namhist1hg,"PedRefhg_%d",j-22);
130        sprintf(namhist2hg,"PedRefhgOutOfTime_%d",j-22);
131        sprintf(namhist3hg,"PedCorrRefhg_%d",j-22);
132        //
133        sprintf(namhist1lg,"PedReflg_%d",j-22);
134        sprintf(namhist2lg,"PedReflgOutOfTime_%d",j-22);
135        sprintf(namhist3lg,"PedCorrReflg_%d",j-22);
136      }
137      // --- High gain chain histos
138      hPedhg[j] = new TH1F(namhist1hg, namhist1hg, 200,0., 200.);
139      hPedOutOfTimehg[j] = new TH1F(namhist2hg, namhist2hg, 200,0., 200.);
140      hPedCorrhg[j] = new TH2F(namhist3hg,namhist3hg,100,0.,200.,100,0.,200.);
141      // --- Low gain chain histos
142      hPedlg[j] = new TH1F(namhist1lg, namhist1lg, 100,0., 600.);
143      hPedOutOfTimelg[j] = new TH1F(namhist2lg, namhist2lg, 100,0., 600.);
144      hPedCorrlg[j] = new TH2F(namhist3lg,namhist3lg,100,0.,600.,100,0.,600.);
145   }
146
147
148   /* open result file */
149   FILE *fp=NULL;
150   fp=fopen("./result.txt","a");
151   if (fp==NULL) {
152     printf("Failed to open file\n");
153     return -1;
154   }
155   
156   FILE *mapFile4Shuttle;
157
158   /* report progress */
159   daqDA_progressReport(10);
160
161
162   /* init some counters */
163   int nevents_physics=0;
164   int nevents_total=0;
165
166   /* read the data files */
167   int n;
168   for(n=1;n<argc;n++){
169    
170     status=monitorSetDataSource( argv[n] );
171     if (status!=0) {
172       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
173       return -1;
174     }
175
176     /* report progress */
177     /* in this example, indexed on the number of files */
178     daqDA_progressReport(10+80*n/argc);
179
180     /* read the file */
181     for(;;) {
182       struct eventHeaderStruct *event;
183       eventTypeType eventT;
184
185       /* get next event */
186       status=monitorGetEventDynamic((void **)&event);
187       if(status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
188       if(status!=0) {
189         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
190         return -1;
191       }
192
193       /* retry if got no event */
194       if(event==NULL) {
195         break;
196       }
197       
198       // Initalize raw-data reading and decoding
199       AliRawReader *reader = new AliRawReaderDate((void*)event);
200       reader->Select("ZDC");
201       // --- Reading event header
202       //UInt_t evtype = reader->GetType();
203       //printf("\n\t ZDCPEDESTALda -> ev. type %d\n",evtype);
204       //printf("\t ZDCPEDESTALda -> run # %d\n",reader->GetRunNumber());
205       //
206       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
207         
208
209       /* use event - here, just write event id to result file */
210       eventT=event->eventType;
211       
212       Int_t ich=0, adcMod[48], adcCh[48], sigCode[48], det[48], sec[48];
213       
214       if(eventT==START_OF_DATA){
215                 
216         rawStreamZDC->SetSODReading(kTRUE);
217         
218         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
219         else{
220           while(rawStreamZDC->Next()){
221             if(rawStreamZDC->IsChMapping()){
222               adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);
223               adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);
224               sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
225               det[ich] = rawStreamZDC->GetDetectorFromMap(ich);
226               sec[ich] = rawStreamZDC->GetTowerFromMap(ich);
227               ich++;
228             }
229           }
230         }
231         // --------------------------------------------------------
232         // --- Writing ascii data file for the Shuttle preprocessor
233         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
234         for(Int_t i=0; i<ich; i++){
235            fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",i,
236              adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
237            //
238            //printf("ZDCPEDESTALDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
239            //      i,adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
240         }
241         fclose(mapFile4Shuttle);
242       }
243       
244       if(eventT==PHYSICS_EVENT){
245         fprintf(fp,"Run #%lu, event size: %lu, BC:%u, Orbit:%u, Period:%u\n",
246           (unsigned long)event->eventRunNb,
247           (unsigned long)event->eventSize,
248           EVENT_ID_GET_BUNCH_CROSSING(event->eventId),
249           EVENT_ID_GET_ORBIT(event->eventId),
250           EVENT_ID_GET_PERIOD(event->eventId));
251
252         // --- Reading data header
253         reader->ReadHeader();
254         const AliRawDataHeader* header = reader->GetDataHeader();
255         if(header){
256          UChar_t message = header->GetAttributes();
257          if(message & 0x20){ // PEDESTAL RUN
258             //printf("\t STANDALONE_PEDESTAL RUN raw data found\n");
259          }
260          else{
261             printf("\t NO STANDALONE_PEDESTAL RUN raw data found\n");
262             return -1;
263          }
264         }
265         else{
266            printf("\t ATTENTION! No Raw Data Header found!!!\n");
267            return -1;
268         }
269         
270         rawStreamZDC->SetSODReading(kTRUE);
271
272         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n"); 
273         //
274         // ----- Setting ch. mapping -----
275         for(Int_t jk=0; jk<48; jk++){
276           //printf("ZDCPEDESTALDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
277           //    jk,adcMod[jk],adcCh[jk],sigCode[jk],det[jk],sec[jk]);
278           rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
279           rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
280           rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
281           rawStreamZDC->SetMapDet(jk, det[jk]);
282           rawStreamZDC->SetMapTow(jk, sec[jk]);
283         }
284         //
285         Int_t iraw=0;
286         Int_t RawADChg[kNChannels], RawADCoothg[kNChannels];
287         Int_t RawADClg[kNChannels], RawADCootlg[kNChannels];
288         for(Int_t j=0; j<kNChannels; j++){
289            RawADChg[j]=0; RawADCoothg[j]=0;
290            RawADClg[j]=0; RawADCootlg[j]=0;
291         }
292         //
293         while(rawStreamZDC->Next()){
294          Int_t index=-1;
295          Int_t detector = rawStreamZDC->GetSector(0);
296          if(rawStreamZDC->IsADCDataWord() && (detector!=-1)){
297           if(rawStreamZDC->GetSector(1)!=5){ // Physics signals
298             if(detector==1) index = rawStreamZDC->GetSector(1); // *** ZNC
299             else if(detector==2) index = rawStreamZDC->GetSector(1)+5; // *** ZPC
300             else if(detector==3) index = rawStreamZDC->GetSector(1)+9; // *** ZEM
301             else if(detector==4) index = rawStreamZDC->GetSector(1)+12; // *** ZNA
302             else if(detector==5) index = rawStreamZDC->GetSector(1)+17; // *** ZPA
303           }
304           else{ // Reference PMs
305             index = (detector-1)/3+22;
306           }
307           //
308           if(index==-1) printf("\tERROR!!! iraw %d det %d quad %d res %d index %d ADC %d\n", 
309             iraw, detector, rawStreamZDC->GetSector(1), 
310             rawStreamZDC->GetADCGain(), index, rawStreamZDC->GetADCValue());
311           
312            //
313            if(iraw<2*kNChannels){ // --- In-time pedestals (1st 48 raw data)
314             if(rawStreamZDC->GetADCGain()==0){ 
315               hPedhg[index]->Fill(rawStreamZDC->GetADCValue()); 
316               RawADChg[index] = rawStreamZDC->GetADCValue();
317               //
318               //printf("\t filling histo hPedhg[%d]\n",index);
319             }
320             else{
321               hPedlg[index]->Fill(rawStreamZDC->GetADCValue()); 
322               RawADClg[index] = rawStreamZDC->GetADCValue();
323               //
324               //printf("\t filling histo hPedlg[%d]\n",index);
325             }
326            }
327            else{  // --- Out-of-time pedestals
328             if(rawStreamZDC->GetADCGain()==0){
329               hPedOutOfTimehg[index]->Fill(rawStreamZDC->GetADCValue());
330               RawADCoothg[index] = rawStreamZDC->GetADCValue();
331               //
332               //printf("\t filling histo hPedOutOfTimehg[%d]\n",index);
333             }
334             else{
335               hPedOutOfTimelg[index]->Fill(rawStreamZDC->GetADCValue());
336               RawADCootlg[index] = rawStreamZDC->GetADCValue();
337               //
338               //printf("\t filling histo hPedOutOfTimelg[%d]\n",index);
339             }
340            }
341             iraw++;
342           }//IsADCDataWord()
343           //
344           if(iraw == 4*kNChannels){ // Last ADC channel -> Filling correlation histos
345             for(Int_t k=0; k<kNChannels; k++){
346               hPedCorrhg[k]->Fill(RawADCoothg[k], RawADChg[k]);
347               hPedCorrlg[k]->Fill(RawADCootlg[k], RawADClg[k]);
348             }
349           }
350          }
351          //
352          nevents_physics++;
353          //
354          delete reader;
355          delete rawStreamZDC;
356
357       }//(if PHYSICS_EVENT) 
358       nevents_total++;
359
360       /* free resources */
361       free(event);
362     
363     }
364   }  
365   
366   /* Analysis of the histograms */
367   //
368   FILE *fileShuttle;
369   fileShuttle = fopen(PEDDATA_FILE,"w");
370   //
371   Float_t MeanPed[2*kNChannels], MeanPedWidth[2*kNChannels], 
372         MeanPedOOT[2*kNChannels], MeanPedWidthOOT[2*kNChannels];
373   // --- In-time pedestals
374   TF1 *ADCfunchg[kNChannels];
375   for(Int_t i=0; i<kNChannels; i++){
376      hPedhg[i]->Fit("gaus","Q");
377      ADCfunchg[i] = hPedhg[i]->GetFunction("gaus");
378      MeanPed[i] = (Double_t) ADCfunchg[i]->GetParameter(1);
379      MeanPedWidth[i] = (Double_t)  ADCfunchg[i]->GetParameter(2);
380      fprintf(fileShuttle,"\t%f\t%f\n",MeanPed[i],MeanPedWidth[i]);
381      printf("\t MeanPedhg[%d] = %f\n",i, MeanPed[i]);
382   }
383   TF1 *ADCfunclg[kNChannels];
384   for(Int_t i=0; i<kNChannels; i++){
385      hPedlg[i]->Fit("gaus","Q");
386      ADCfunclg[i] = hPedlg[i]->GetFunction("gaus");
387      MeanPed[i+kNChannels] = (Double_t)  ADCfunclg[i]->GetParameter(1);
388      MeanPedWidth[i+kNChannels] = (Double_t)  ADCfunclg[i]->GetParameter(2);
389      fprintf(fileShuttle,"\t%f\t%f\n",MeanPed[i+kNChannels],MeanPedWidth[i+kNChannels]);
390      printf("\t MeanPedlg[%d] = %f\n",i+kNChannels, MeanPed[i+kNChannels]);
391   }
392   // --- Out-of-time pedestals
393   TF1 *ADCootfunchg[kNChannels];
394   for(Int_t i=0; i<kNChannels; i++){
395      hPedOutOfTimehg[i]->Fit("gaus","Q");
396      ADCootfunchg[i] = hPedOutOfTimehg[i]->GetFunction("gaus");
397      MeanPedOOT[i] = (Double_t)  ADCootfunchg[i]->GetParameter(1);
398      MeanPedWidthOOT[i] = (Double_t)  ADCootfunchg[i]->GetParameter(2);
399      fprintf(fileShuttle,"\t%f\t%f\n",MeanPedOOT[i],MeanPedWidthOOT[i]);
400      printf("\t MeanPedOOThg[%d] = %f\n",i, MeanPedOOT[i]);
401   }
402   TF1 *ADCootfunclg[kNChannels];
403   for(Int_t i=0; i<kNChannels; i++){
404      hPedOutOfTimelg[i]->Fit("gaus","Q");
405      ADCootfunclg[i] = hPedOutOfTimelg[i]->GetFunction("gaus");
406      MeanPedOOT[i+kNChannels] = (Double_t)  ADCootfunclg[i]->GetParameter(1);
407      MeanPedWidthOOT[i+kNChannels] = (Double_t)  ADCootfunclg[i]->GetParameter(2);
408      fprintf(fileShuttle,"\t%f\t%f\n",MeanPedOOT[i+kNChannels],MeanPedWidthOOT[i+kNChannels]);
409      printf("\t MeanPedOOTlg[%d] = %f\n",i+kNChannels, MeanPedOOT[i+kNChannels]);
410   }
411   
412   // ***************************************************
413   //   Unless we have a narrow correlation to fit we
414   //    don't fit and store in-time vs. out-of-time
415   //    histograms -> mean pedstal subtracted!!!!!!
416   // ***************************************************
417   // --- Correlations
418 /*
419   Float_t CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
420   TProfile *hPedCorrProfhg[kNChannels], *hPedCorrProflg[kNChannels];
421   TF1 *ffunchg[kNChannels], *ffunclg[kNChannels];
422   char namhist4[50];
423   for(int i=0;i<kNChannels;i++) {
424      sprintf(namhist4,"ADCHRvsOOT%d_Prof",i);
425      hPedCorrProfhg[i] = hPedCorrhg[i]->ProfileX(namhist4,-1,-1,"S");
426      hPedCorrProfhg[i]->SetName(namhist4);
427      hPedCorrProfhg[i]->Fit("pol1","Q");
428      ffunchg[i] = hPedCorrProfhg[i]->GetFunction("pol1");
429      CorrCoeff0[i] = (Double_t)  ffunchg[i]->GetParameter(0);
430      CorrCoeff1[i] = (Double_t) ffunchg[i]->GetParameter(1);
431      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i,CorrCoeff0[i],CorrCoeff1[i]);
432      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",i, CorrCoeff0[i], i, CorrCoeff1[i]);
433   }    
434   for(int i=0;i<kNChannels;i++) {
435      sprintf(namhist4,"ADCLRvsOOT%d_Prof",i);
436      hPedCorrProflg[i] = hPedCorrlg[i]->ProfileX(namhist4,-1,-1,"S");
437      hPedCorrProflg[i]->SetName(namhist4);
438      hPedCorrProflg[i]->Fit("pol1","Q");
439      ffunclg[i] = hPedCorrProflg[i]->GetFunction("pol1");
440      CorrCoeff0[i+kNChannels] =  (Double_t) ffunclg[i]->GetParameter(0);
441      CorrCoeff1[i+kNChannels] =  (Double_t) ffunclg[i]->GetParameter(1);
442      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i+kNChannels,CorrCoeff0[i+kNChannels],CorrCoeff1[i+kNChannels]);
443      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",
444      //         i+kNChannels, CorrCoeff0[i+kNChannels], i+kNChannels, CorrCoeff1[i+kNChannels]);
445   }    
446 */
447   //                                                   
448   fclose(fileShuttle);
449   //
450   for(Int_t j=0; j<kNChannels; j++){
451      delete hPedhg[j];
452      delete hPedOutOfTimehg[j];
453      delete hPedCorrhg[j];
454      delete hPedlg[j];
455      delete hPedOutOfTimelg[j];
456      delete hPedCorrlg[j];
457   }
458
459   //delete minuitFit;
460   TVirtualFitter::SetFitter(0);
461
462   /* write report */
463   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
464
465   /* close result file */
466   fclose(fp);
467   
468   /* report progress */
469   daqDA_progressReport(90);
470
471   /* store the result files on FES */
472   // [1] File with mapping
473   status = daqDA_FES_storeFile(MAPDATA_FILE,MAPDATA_FILE);
474   if(status){
475     printf("Failed to export file %d to DAQ FES\n",status);
476     return -1;
477   }
478   // [2] File with pedestal data
479   status = daqDA_FES_storeFile(PEDDATA_FILE,PEDDATA_FILE);
480   if(status){
481     printf("Failed to export file %d to DAQ FES\n",status);
482     return -1;
483   }
484   
485   /* store the result files on DB */
486   status = daqDA_DB_storeFile(PEDDATA_FILE,PEDDATA_FILE);  
487   if(status){
488     printf("Failed to export file %d to DAQ DB\n",status);
489     return -1;
490   }
491
492   /* report progress */
493   daqDA_progressReport(100);
494
495   return status;
496 }