DAs upgraded to AliZDCRawStream class
[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         rawStreamZDC->SetSODReading(kTRUE);
195                 
196         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
197         else{
198           while(rawStreamZDC->Next()){
199             if(rawStreamZDC->IsChMapping()){
200               adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);
201               adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);
202               sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
203               det[ich] = rawStreamZDC->GetDetectorFromMap(ich);
204               sec[ich] = rawStreamZDC->GetTowerFromMap(ich);
205               ich++;
206             }
207           }
208         }
209         // --------------------------------------------------------
210         // --- Writing ascii data file for the Shuttle preprocessor
211         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
212         for(Int_t i=0; i<ich; i++){
213            fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",i,
214              adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
215            //
216            //printf("ZDCLASERDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
217            //      i,adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
218         }
219         fclose(mapFile4Shuttle);
220       }
221
222       /* use event - here, just write event id to result file */
223       eventT=event->eventType;
224     
225       if(eventT==PHYSICS_EVENT){
226         //
227         // --- Reading data header
228         reader->ReadHeader();
229         const AliRawDataHeader* header = reader->GetDataHeader();
230         if(header) {
231          UChar_t message = header->GetAttributes();
232          if(message & 0x20){ // DEDICATED LASER RUN
233             //printf("\t STANDALONE_LASER_RUN raw data found\n");
234             continue;
235          }
236          else{
237             printf("\t NO STANDALONE_LASER_RUN raw data found\n");
238             return -1;
239          }
240         }
241         else{
242            printf("\t ATTENTION! No Raw Data Header found!!!\n");
243            return -1;
244         }
245
246         rawStreamZDC->SetSODReading(kTRUE);
247
248         if (!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
249         //
250         // ----- Setting ch. mapping -----
251         for(Int_t jk=0; jk<48; jk++){
252           rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
253           rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
254           rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
255           rawStreamZDC->SetMapDet(jk, det[jk]);
256           rawStreamZDC->SetMapTow(jk, sec[jk]);
257         }
258         //
259         while(rawStreamZDC->Next()){
260           Int_t index=-1;
261           // Getting data only for reference PMTs (sector[1]=5)
262           if((rawStreamZDC->IsADCDataWord()) && (rawStreamZDC->GetSector(1)==5)){
263             index = rawStreamZDC->GetADCChannel();
264             Float_t Pedestal = MeanPed[index];
265             Float_t CorrADC = rawStreamZDC->GetADCValue() - Pedestal;
266             
267             // ==== HIGH GAIN CHAIN
268             if(rawStreamZDC->GetADCGain() == 0){
269               // %%%%% PMRef chain side C
270               if(rawStreamZDC->GetSector(0)==1) hPMRefChg->Fill(CorrADC);
271               // %%%%% PMRef side A
272               else if(rawStreamZDC->GetSector(0)==4) hPMRefAhg->Fill(CorrADC);
273             }
274             // ==== LOW GAIN CHAIN
275             else{
276               // %%%%% PMRef chain side C
277               if(rawStreamZDC->GetSector(0)==1) hPMRefClg->Fill(CorrADC);
278               // %%%%% PMRef side A
279               else if(rawStreamZDC->GetSector(0)==4) hPMRefAlg->Fill(CorrADC);
280             }
281           }//IsADCDataWord()
282           //
283          }
284          //
285          nevents_physics++;
286          //
287          delete reader;
288          delete rawStreamZDC;
289
290       }//(if PHYSICS_EVENT) 
291       nevents_total++;
292
293       /* free resources */
294       free(event);
295     
296     }
297   }  
298   
299   /* Analysis of the histograms */
300   //
301   Int_t maxBinRef[4], nBinRef[4];
302   Float_t xMaxRef[4], maxXvalRef[4], xlowRef[4]; 
303   Float_t meanRef[2], sigmaRef[2];
304   TF1 *funRef[4];
305   
306   // ~~~~~~~~ PM Ref side C high gain chain ~~~~~~~~
307   maxBinRef[0] = hPMRefChg->GetMaximumBin();
308   nBinRef[0] = (hPMRefChg->GetXaxis())->GetNbins();
309   xMaxRef[0] = (hPMRefChg->GetXaxis())->GetXmax();
310   maxXvalRef[0] = maxBinRef[0]*xMaxRef[0]/nBinRef[0];
311   // 
312   if(maxXvalRef[0]-100.<0.) {xlowRef[0]=0.;}
313   else xlowRef[0] = maxXvalRef[0];
314   hPMRefChg->Fit("gaus","Q","",xlowRef[0],maxXvalRef[0]+100.);
315   funRef[0] = hPMRefChg->GetFunction("gaus");
316   meanRef[0] = (Float_t) (funRef[0]->GetParameter(1));
317   sigmaRef[0] = (Float_t) (funRef[0]->GetParameter(2));
318   
319   // ~~~~~~~~ PM Ref side A high gain chain ~~~~~~~~
320   maxBinRef[1] = hPMRefAhg->GetMaximumBin();
321   nBinRef[1] = (hPMRefAhg->GetXaxis())->GetNbins();
322   xMaxRef[1] = (hPMRefAhg->GetXaxis())->GetXmax();
323   maxXvalRef[1] = maxBinRef[1]*xMaxRef[1]/nBinRef[1];
324   //
325   if(maxXvalRef[1]-100.<0.) {xlowRef[1]=0.;}
326   else xlowRef[1] = maxXvalRef[1];
327   hPMRefAhg->Fit("gaus","Q","",xlowRef[1],maxXvalRef[1]+100.);
328   funRef[1] = hPMRefAhg->GetFunction("gaus");
329   meanRef[1] = (Float_t) (funRef[1]->GetParameter(1));
330   sigmaRef[1] = (Float_t) (funRef[1]->GetParameter(2));
331   
332   // ~~~~~~~~ PM Ref side C low gain chain ~~~~~~~~
333   maxBinRef[2] = hPMRefClg->GetMaximumBin();
334   nBinRef[2] = (hPMRefClg->GetXaxis())->GetNbins();
335   xMaxRef[2] = (hPMRefClg->GetXaxis())->GetXmax();
336   maxXvalRef[2] = maxBinRef[2]*xMaxRef[2]/nBinRef[2];
337   //
338   if(maxXvalRef[2]-100.<0.) {xlowRef[2]=0.;}
339   else xlowRef[2] = maxXvalRef[2];
340   hPMRefClg->Fit("gaus","Q","",xlowRef[2],maxXvalRef[2]+100.);
341   funRef[2] = hPMRefClg->GetFunction("gaus");
342   meanRef[2] = (Float_t) (funRef[2]->GetParameter(1));
343   sigmaRef[2] = (Float_t) (funRef[2]->GetParameter(2));
344   
345   // ~~~~~~~~ PM Ref side A low gain chain ~~~~~~~~
346   maxBinRef[3] = hPMRefAlg->GetMaximumBin();
347   nBinRef[3] = (hPMRefAlg->GetXaxis())->GetNbins();
348   xMaxRef[3] = (hPMRefAlg->GetXaxis())->GetXmax();
349   maxXvalRef[3] = maxBinRef[3]*xMaxRef[3]/nBinRef[3];
350   //
351   if(maxXvalRef[3]-100.<0.) {xlowRef[3]=0.;}
352   else xlowRef[3] = maxXvalRef[3];
353   hPMRefAlg->Fit("gaus","Q","",xlowRef[3],maxXvalRef[3]+100.);
354   funRef[3] = hPMRefAlg->GetFunction("gaus");
355   meanRef[3] = (Float_t) (funRef[3]->GetParameter(1));
356   sigmaRef[3] = (Float_t) (funRef[3]->GetParameter(2));
357   //
358   FILE *fileShuttle;
359   fileShuttle = fopen(LASDATA_FILE,"w");
360   for(Int_t i=0; i<4; i++)  fprintf(fileShuttle,"\t%f\t%f\n",meanRef[i], sigmaRef[i]); 
361   //                                                   
362   fclose(fileShuttle);
363   //
364   delete hPMRefChg;
365   delete hPMRefAhg;
366   delete hPMRefClg;
367   delete hPMRefAlg;
368
369   //delete minuitFit;
370   TVirtualFitter::SetFitter(0);
371
372   /* write report */
373   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
374
375   /* close result file */
376   fclose(fp);
377   
378   /* report progress */
379   daqDA_progressReport(90);
380
381   /* store the result file on FES */
382   status = daqDA_FES_storeFile(MAPDATA_FILE,MAPDATA_FILE);
383   if(status){
384     printf("Failed to export file : %d\n",status);
385     return -1;
386   }
387   //
388   status = daqDA_FES_storeFile(LASDATA_FILE,LASDATA_FILE);
389   if(status){
390     printf("Failed to export file : %d\n",status);
391     return -1;
392   }
393
394   /* report progress */
395   daqDA_progressReport(100);
396
397
398   return status;
399 }