9c997bf17592190871586d822c5d11bf417597bb
[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
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <Riostream.h>
28
29 // DATE
30 #include <event.h>
31 #include <monitor.h>
32 #include <daqDA.h>
33
34 //ROOT
35 #include <TRandom.h>
36 #include <TH1F.h>
37 #include <TH2F.h>
38 #include <TProfile.h>
39 #include <TF1.h>
40 #include <TFile.h>
41 #include <TFitter.h>
42
43 //AliRoot
44 #include <AliRawReaderDate.h>
45 #include <AliRawEventHeaderBase.h>
46 #include <AliZDCRawStream.h>
47
48
49 /* Main routine
50       Arguments: list of DATE raw data files
51 */
52 int main(int argc, char **argv) {
53   
54   TFitter *minuitFit = new TFitter(4);
55   TVirtualFitter::SetFitter(minuitFit);
56
57   int status = 0;
58
59   /* log start of process */
60   printf("ZDC PEDESTAL program started\n");  
61
62   /* check that we got some arguments = list of files */
63   if (argc<2) {
64     printf("Wrong number of arguments\n");
65     return -1;
66   }
67
68   // --- Histograms for ADC pedestals 
69   //     [22 signal channels + 2 reference PTMs]  x 2 gain chains
70   //
71   TH1F::AddDirectory(0);
72   int const kNChannels = 24;
73   TH1F *hPedhg[kNChannels], *hPedOutOfTimehg[kNChannels];
74   TH2F *hPedCorrhg[kNChannels];
75   TH1F *hPedlg[kNChannels], *hPedOutOfTimelg[kNChannels];
76   TH2F *hPedCorrlg[kNChannels];
77   //
78   char namhist1hg[50], namhist2hg[50], namhist3hg[50];
79   char namhist1lg[50], namhist2lg[50], namhist3lg[50];
80   for(Int_t j=0; j<kNChannels; j++){
81      if(j<5){ // ZN1
82        sprintf(namhist1hg,"PedZN1hg_%d",j);
83        sprintf(namhist2hg,"PedZN1hgOutOfTime_%d",j);
84        sprintf(namhist3hg,"PedCorrZN1hg_%d",j);
85        //
86        sprintf(namhist1lg,"PedZN1lg_%d",j);
87        sprintf(namhist2lg,"PedZN1lgOutOfTime_%d",j);
88        sprintf(namhist3lg,"PedCorrZN1lg_%d",j);
89      }
90      else if(j>=5 && j<10){ // ZP1
91        sprintf(namhist1hg,"PedZP1hg_%d",j-5);
92        sprintf(namhist2hg,"PedZP1hgOutOfTime_%d",j-5);
93        sprintf(namhist3hg,"PedCorrZP1hg_%d",j-5);
94        //
95        sprintf(namhist1lg,"PedZP1lg_%d",j-5);
96        sprintf(namhist2lg,"PedZP1lgOutOfTime_%d",j-5);
97        sprintf(namhist3lg,"PedCorrZP1lg_%d",j-5);       
98      }
99      else if(j>=10 && j<12){ // ZEM
100        sprintf(namhist1hg,"PedZEMhg_%d",j-10);
101        sprintf(namhist2hg,"PedZEMhgOutOfTime_%d",j-10);
102        sprintf(namhist3hg,"PedCorrZEMhg_%d",j-10);
103        //
104        sprintf(namhist1lg,"PedZEMlg_%d",j-10);
105        sprintf(namhist2lg,"PedZEMlgOutOfTime_%d",j-10);
106        sprintf(namhist3lg,"PedCorrZEMlg_%d",j-10);
107      }
108      else if(j>=12 && j<17){ // ZN2
109        sprintf(namhist1hg,"PedZN2hg_%d",j-12);
110        sprintf(namhist2hg,"PedZN2hgOutOfTime_%d",j-12);
111        sprintf(namhist3hg,"PedCorrZN2hg_%d",j-12);
112        //
113        sprintf(namhist1lg,"PedZN2lg_%d",j-12);
114        sprintf(namhist2lg,"PedZN2lgOutOfTime_%d",j-12);
115        sprintf(namhist3lg,"PedCorrZN2lg_%d",j-12);
116      }
117      else if(j>=17 && j<22){ // ZP2
118        sprintf(namhist1hg,"PedZP2hg_%d",j-17);
119        sprintf(namhist2hg,"PedZP2hgOutOfTime_%d",j-17);
120        sprintf(namhist3hg,"PedCorrZP2hg_%d",j-17);
121        //
122        sprintf(namhist1lg,"PedZP2lg_%d",j-17);
123        sprintf(namhist2lg,"PedZP2lgOutOfTime_%d",j-17);
124        sprintf(namhist3lg,"PedCorrZP2lg_%d",j-17);
125      }
126      else if(j>=22 && j<24){ //Reference PMs
127        sprintf(namhist1hg,"PedRefhg_%d",j-22);
128        sprintf(namhist2hg,"PedRefhgOutOfTime_%d",j-22);
129        sprintf(namhist3hg,"PedCorrRefhg_%d",j-22);
130        //
131        sprintf(namhist1lg,"PedReflg_%d",j-22);
132        sprintf(namhist2lg,"PedReflgOutOfTime_%d",j-22);
133        sprintf(namhist3lg,"PedCorrReflg_%d",j-22);
134      }
135      // --- High gain chain histos
136      hPedhg[j] = new TH1F(namhist1hg, namhist1hg, 200,0., 200.);
137      hPedOutOfTimehg[j] = new TH1F(namhist2hg, namhist2hg, 200,0., 200.);
138      hPedCorrhg[j] = new TH2F(namhist3hg,namhist3hg,100,0.,200.,100,0.,200.);
139      // --- Low gain chain histos
140      hPedlg[j] = new TH1F(namhist1lg, namhist1lg, 100,0., 600.);
141      hPedOutOfTimelg[j] = new TH1F(namhist2lg, namhist2lg, 100,0., 600.);
142      hPedCorrlg[j] = new TH2F(namhist3lg,namhist3lg,100,0.,600.,100,0.,600.);
143   }
144
145
146   /* open result file */
147   FILE *fp=NULL;
148   fp=fopen("./result.txt","a");
149   if (fp==NULL) {
150     printf("Failed to open file\n");
151     return -1;
152   }
153   
154   FILE *mapFile4Shuttle;
155   const char *mapfName = "ZDCChMapping.dat";
156   
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       if(eventT==START_OF_DATA){
214                 
215         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
216         else{
217           while(rawStreamZDC->Next()){
218             if(rawStreamZDC->IsChMapping()){
219               adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);
220               adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);
221               sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
222               det[ich] = rawStreamZDC->GetDetectorFromMap(ich);
223               sec[ich] = rawStreamZDC->GetTowerFromMap(ich);
224               ich++;
225             }
226           }
227         }
228         // --------------------------------------------------------
229         // --- Writing ascii data file for the Shuttle preprocessor
230         mapFile4Shuttle = fopen(mapfName,"w");
231         for(Int_t i=0; i<ich; i++){
232            fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",i,
233              adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
234            //
235            //printf("ZDCPEDESTALDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
236            //      i,adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
237         }
238         fclose(mapFile4Shuttle);
239       }
240       
241       if(eventT==PHYSICS_EVENT){
242         fprintf(fp,"Run #%lu, event size: %lu, BC:%u, Orbit:%u, Period:%u\n",
243           (unsigned long)event->eventRunNb,
244           (unsigned long)event->eventSize,
245           EVENT_ID_GET_BUNCH_CROSSING(event->eventId),
246           EVENT_ID_GET_ORBIT(event->eventId),
247           EVENT_ID_GET_PERIOD(event->eventId));
248
249         // --- Reading data header
250         reader->ReadHeader();
251         const AliRawDataHeader* header = reader->GetDataHeader();
252         if(header){
253          UChar_t message = header->GetAttributes();
254          if(message & 0x20){ // PEDESTAL RUN
255             //printf("\t STANDALONE_PEDESTAL RUN raw data found\n");
256          }
257          else{
258             printf("\t NO STANDALONE_PEDESTAL RUN raw data found\n");
259             return -1;
260          }
261         }
262         else{
263            printf("\t ATTENTION! No Raw Data Header found!!!\n");
264            return -1;
265         }
266
267         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n"); 
268         //
269         // ----- Setting ch. mapping -----
270         for(Int_t jk=0; jk<48; jk++){
271           rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
272           rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
273           rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
274           rawStreamZDC->SetMapDet(jk, det[jk]);
275           rawStreamZDC->SetMapTow(jk, sec[jk]);
276         }
277         //
278         Int_t iraw=0;
279         Int_t RawADChg[kNChannels], RawADCoothg[kNChannels];
280         Int_t RawADClg[kNChannels], RawADCootlg[kNChannels];
281         for(Int_t j=0; j<kNChannels; j++){
282            RawADChg[j]=0; RawADCoothg[j]=0;
283            RawADClg[j]=0; RawADCootlg[j]=0;
284         }
285         //
286         while(rawStreamZDC->Next()){
287          Int_t index=-1;
288          if(rawStreamZDC->IsADCDataWord()){
289           if(rawStreamZDC->GetSector(1)!=5){ // Physics signals
290             if(rawStreamZDC->GetSector(0)==1) index = rawStreamZDC->GetSector(1); // *** ZNC
291             else if(rawStreamZDC->GetSector(0)==2) index = rawStreamZDC->GetSector(1)+5; // *** ZPC
292             else if(rawStreamZDC->GetSector(0)==3) index = rawStreamZDC->GetSector(1)+9; // *** ZEM
293             else if(rawStreamZDC->GetSector(0)==4) index = rawStreamZDC->GetSector(1)+12; // *** ZNA
294             else if(rawStreamZDC->GetSector(0)==5) index = rawStreamZDC->GetSector(1)+17; // *** ZPA
295           }
296           else{ // Reference PMs
297             index = (rawStreamZDC->GetSector(0)-1)/3+22;
298           }
299           //
300           /*printf("\t iraw %d index %d det %d quad %d res %d ADC %d\n", iraw, index,
301                 rawStreamZDC->GetSector(0), rawStreamZDC->GetSector(1), 
302                 rawStreamZDC->GetADCGain(), rawStreamZDC->GetADCValue());
303           */
304            //
305            if(iraw<2*kNChannels){ // --- In-time pedestals (1st 48 raw data)
306             if(rawStreamZDC->GetADCGain()==0){ 
307               hPedhg[index]->Fill(rawStreamZDC->GetADCValue()); 
308               RawADChg[index] = rawStreamZDC->GetADCValue();
309               //
310               //printf("\t filling histo hPedhg[%d]\n",index);
311             }
312             else{
313               hPedlg[index]->Fill(rawStreamZDC->GetADCValue()); 
314               RawADClg[index] = rawStreamZDC->GetADCValue();
315               //
316               //printf("\t filling histo hPedlg[%d]\n",index);
317             }
318            }
319            else{  // --- Out-of-time pedestals
320             if(rawStreamZDC->GetADCGain()==0){
321               hPedOutOfTimehg[index]->Fill(rawStreamZDC->GetADCValue());
322               RawADCoothg[index] = rawStreamZDC->GetADCValue();
323               //
324               //printf("\t filling histo hPedOutOfTimehg[%d]\n",index);
325             }
326             else{
327               hPedOutOfTimelg[index]->Fill(rawStreamZDC->GetADCValue());
328               RawADCootlg[index] = rawStreamZDC->GetADCValue();
329               //
330               //printf("\t filling histo hPedOutOfTimelg[%d]\n",index);
331             }
332            }
333             iraw++;
334           }//IsADCDataWord()
335           //
336           if(iraw == 4*kNChannels){ // Last ADC channel -> Filling correlation histos
337             for(Int_t k=0; k<kNChannels; k++){
338               hPedCorrhg[k]->Fill(RawADCoothg[k], RawADChg[k]);
339               hPedCorrlg[k]->Fill(RawADCootlg[k], RawADClg[k]);
340             }
341           }
342          }
343          //
344          nevents_physics++;
345          //
346          delete reader;
347          delete rawStreamZDC;
348
349       }//(if PHYSICS_EVENT) 
350       nevents_total++;
351
352       /* free resources */
353       free(event);
354     
355     }
356   }  
357   
358   /* Analysis of the histograms */
359   //
360   FILE *fileShuttle;
361   const char *fName = "ZDCPedestal.dat";
362   fileShuttle = fopen(fName,"w");
363   //
364   Float_t MeanPed[2*kNChannels], MeanPedWidth[2*kNChannels], 
365         MeanPedOOT[2*kNChannels], MeanPedWidthOOT[2*kNChannels],
366         CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
367   // --- In-time pedestals
368   TF1 *ADCfunchg[kNChannels];
369   for(Int_t i=0; i<kNChannels; i++){
370      hPedhg[i]->Fit("gaus","Q");
371      ADCfunchg[i] = hPedhg[i]->GetFunction("gaus");
372      MeanPed[i] = (Double_t) ADCfunchg[i]->GetParameter(1);
373      MeanPedWidth[i] = (Double_t)  ADCfunchg[i]->GetParameter(2);
374      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i,MeanPed[i],MeanPedWidth[i]);
375      //printf("\t MeanPed[%d] = %f\n",i, MeanPed[i]);
376   }
377   TF1 *ADCfunclg[kNChannels];
378   for(Int_t i=0; i<kNChannels; i++){
379      hPedlg[i]->Fit("gaus","Q");
380      ADCfunclg[i] = hPedlg[i]->GetFunction("gaus");
381      MeanPed[i+kNChannels] = (Double_t)  ADCfunclg[i]->GetParameter(1);
382      MeanPedWidth[i+kNChannels] = (Double_t)  ADCfunclg[i]->GetParameter(2);
383      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i+kNChannels,MeanPed[i+kNChannels],MeanPedWidth[i+kNChannels]);
384      //printf("\t MeanPed[%d] = %f\n",i+kNChannels, MeanPed[i+kNChannels]);
385   }
386   // --- Out-of-time pedestals
387   TF1 *ADCootfunchg[kNChannels];
388   for(Int_t i=0; i<kNChannels; i++){
389      hPedOutOfTimehg[i]->Fit("gaus","Q");
390      ADCootfunchg[i] = hPedOutOfTimehg[i]->GetFunction("gaus");
391      MeanPedOOT[i] = (Double_t)  ADCootfunchg[i]->GetParameter(1);
392      MeanPedWidthOOT[i] = (Double_t)  ADCootfunchg[i]->GetParameter(2);
393      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i,MeanPedOOT[i],MeanPedWidthOOT[i]);
394      //printf("\t MeanPedOOT[%d] = %f\n",i, MeanPedOOT[i]);
395   }
396   TF1 *ADCootfunclg[kNChannels];
397   for(Int_t i=0; i<kNChannels; i++){
398      hPedOutOfTimelg[i]->Fit("gaus","Q");
399      ADCootfunclg[i] = hPedOutOfTimelg[i]->GetFunction("gaus");
400      MeanPedOOT[i+kNChannels] = (Double_t)  ADCootfunclg[i]->GetParameter(1);
401      MeanPedWidthOOT[i+kNChannels] = (Double_t)  ADCootfunclg[i]->GetParameter(2);
402      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i+kNChannels,MeanPedOOT[i+kNChannels],MeanPedWidthOOT[i+kNChannels]);
403      //printf("\t MeanPedOOT[%d] = %f\n",i+kNChannels, MeanPedOOT[i+kNChannels]);
404   }
405   //
406   // --- Correlations
407   TProfile *hPedCorrProfhg[kNChannels], *hPedCorrProflg[kNChannels];
408   TF1 *ffunchg[kNChannels], *ffunclg[kNChannels];
409   char namhist4[50];
410   for(int i=0;i<kNChannels;i++) {
411      sprintf(namhist4,"ADCHRvsOOT%d_Prof",i);
412      hPedCorrProfhg[i] = hPedCorrhg[i]->ProfileX(namhist4,-1,-1,"S");
413      hPedCorrProfhg[i]->SetName(namhist4);
414      hPedCorrProfhg[i]->Fit("pol1","Q");
415      ffunchg[i] = hPedCorrProfhg[i]->GetFunction("pol1");
416      CorrCoeff0[i] = (Double_t)  ffunchg[i]->GetParameter(0);
417      CorrCoeff1[i] = (Double_t) ffunchg[i]->GetParameter(1);
418      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i,CorrCoeff0[i],CorrCoeff1[i]);
419      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",i, CorrCoeff0[i], i, CorrCoeff1[i]);
420   }    
421   for(int i=0;i<kNChannels;i++) {
422      sprintf(namhist4,"ADCLRvsOOT%d_Prof",i);
423      hPedCorrProflg[i] = hPedCorrlg[i]->ProfileX(namhist4,-1,-1,"S");
424      hPedCorrProflg[i]->SetName(namhist4);
425      hPedCorrProflg[i]->Fit("pol1","Q");
426      ffunclg[i] = hPedCorrProflg[i]->GetFunction("pol1");
427      CorrCoeff0[i+kNChannels] =  (Double_t) ffunclg[i]->GetParameter(0);
428      CorrCoeff1[i+kNChannels] =  (Double_t) ffunclg[i]->GetParameter(1);
429      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i+kNChannels,CorrCoeff0[i+kNChannels],CorrCoeff1[i+kNChannels]);
430      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",
431      //         i+kNChannels, CorrCoeff0[i+kNChannels], i+kNChannels, CorrCoeff1[i+kNChannels]);
432   }    
433   //                                                   
434   fclose(fileShuttle);
435   //
436   for(Int_t j=0; j<kNChannels; j++){
437      delete hPedhg[j];
438      delete hPedOutOfTimehg[j];
439      delete hPedCorrhg[j];
440      delete hPedlg[j];
441      delete hPedOutOfTimelg[j];
442      delete hPedCorrlg[j];
443   }
444
445   //delete minuitFit;
446   TVirtualFitter::SetFitter(0);
447
448   /* write report */
449   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
450
451   /* close result file */
452   fclose(fp);
453   
454   /* report progress */
455   daqDA_progressReport(90);
456
457   /* store the result files on FES */
458   status = daqDA_FES_storeFile(mapfName,"ZDCCHMAPPING_data");
459   if(status){
460     printf("Failed to export file : %d\n",status);
461     return -1;
462   }
463   //
464   status = daqDA_FES_storeFile(fName,"ZDCPEDESTAL_data");
465   if(status){
466     printf("Failed to export file : %d\n",status);
467     return -1;
468   }
469
470   /* report progress */
471   daqDA_progressReport(100);
472
473   return status;
474 }