754f0ff24410f7a6de3c6c811c2bb845ebfbde2e
[u/mrichter/AliRoot.git] / ZDC / ZDCEMDda.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_EMD_RUN
17 DA Type: LDC
18 Number of events needed: at least ~5*10^3
19 Input Files: ZDCPedestal.dat
20 Output Files: ZDCEMDCalib.dat, ZDCChMapping.dat
21 Trigger Types Used: Standalone Trigger
22
23 */
24 #define PEDDATA_FILE  "ZDCPedestal.dat"
25 #define MAPDATA_FILE  "ZDCChMapping.dat"
26 #define EMDDATA_FILE  "ZDCEMDCalib.dat"
27
28 #include <stdio.h>
29 #include <Riostream.h>
30 #include <Riostream.h>
31
32 // DATE
33 #include <daqDA.h>
34 #include <event.h>
35 #include <monitor.h>
36
37 //ROOT
38 #include <TRandom.h>
39 #include <TH1F.h>
40 #include <TH2F.h>
41 #include <TProfile.h>
42 #include <TF1.h>
43 #include <TFile.h>
44 #include <TFitter.h>
45
46 //AliRoot
47 #include <AliRawReaderDate.h>
48 #include <AliRawEventHeaderBase.h>
49 #include <AliZDCRawStream.h>
50
51
52 /* Main routine
53       Arguments: 
54       1- monitoring data source
55 */
56 int main(int argc, char **argv) {
57   
58   TFitter *minuitFit = new TFitter(4);
59   TVirtualFitter::SetFitter(minuitFit);
60
61   int status = 0;
62
63   /* log start of process */
64   printf("ZDC EMD program started\n");  
65
66   /* check that we got some arguments = list of files */
67   if (argc<2) {
68     printf("Wrong number of arguments\n");
69     return -1;
70   }
71   
72   //
73   // --- Preparing histos for EM dissociation spectra
74   //
75   TH1F* histoEMDRaw[4];
76   TH1F* histoEMDCorr[4];
77
78   char namhistr[50], namhistc[50];
79   for(Int_t i=0; i<4; i++) {
80      if(i==0){
81        sprintf(namhistr,"ZN%d-EMDRaw",i+1);
82        sprintf(namhistc,"ZN%d-EMDCorr",i+1);
83      }
84      else if(i==1){
85        sprintf(namhistr,"ZP%d-EMDRaw",i);
86        sprintf(namhistc,"ZP%d-EMDCorr",i);
87      }
88      else if(i==2){
89        sprintf(namhistr,"ZN%d-EMDRaw",i);
90        sprintf(namhistc,"ZN%d-EMDCorr",i);
91      }
92      else if(i==3){
93        sprintf(namhistr,"ZP%d-EMDRaw",i-1);
94        sprintf(namhistc,"ZP%d-EMDCorr",i-1);
95      }
96      histoEMDRaw[i] = new TH1F(namhistr,namhistr,100,0.,4000.);
97      histoEMDCorr[i] = new TH1F(namhistc,namhistc,100,0.,4000.);
98   }
99
100   /* open result file */
101   FILE *fp=NULL;
102   fp=fopen("./result.txt","a");
103   if (fp==NULL) {
104     printf("Failed to open file\n");
105     return -1;
106   }
107   
108   FILE *mapFile4Shuttle;
109
110   // *** To analyze LASER events you MUST have a pedestal data file!!!
111   // *** -> check if a pedestal run has been analyzied
112   int read = 0;
113   read = daqDA_DB_getFile(PEDDATA_FILE,PEDDATA_FILE);
114   if(read){
115     printf("\t ERROR!!! ZDCPedestal.dat file NOT FOUND in DAQ db!!!\n");
116     return -1;
117   }
118   else printf("\t ZDCPedestal.dat file retrieved from DAQ db\n");
119   
120   FILE *filePed = fopen(PEDDATA_FILE,"r");
121   if (filePed==NULL) {
122     printf("\t ERROR!!! Can't open ZDCPedestal.dat file!!!\n");
123     return -1;
124   }
125
126   // 144 = 48 in-time + 48 out-of-time + 48 correlations
127   Float_t readValues[2][144], MeanPed[44], MeanPedWidth[44], 
128         MeanPedOOT[44], MeanPedWidthOOT[44];
129   // ***************************************************
130   //   Unless we have a narrow correlation to fit we
131   //    don't fit and store in-time vs. out-of-time
132   //    histograms -> mean pedstal subtracted!!!!!!
133   // ***************************************************
134   //Float_t CorrCoeff0[44], CorrCoeff1[44];
135   //
136   for(int jj=0; jj<144; jj++){
137     for(int ii=0; ii<2; ii++){
138        fscanf(filePed,"%f",&readValues[ii][jj]);
139     }
140     if(jj<48){
141       MeanPed[jj] = readValues[0][jj];
142       MeanPedWidth[jj] = readValues[1][jj];
143       //printf("\t MeanPed[%d] = %1.1f\n",jj, MeanPed[jj]);
144     }
145     else if(jj>48 && jj<96){
146       MeanPedOOT[jj-48] = readValues[0][jj];
147       MeanPedWidthOOT[jj-48] = readValues[1][jj];
148     }
149     /*else if(jj>144){
150       CorrCoeff0[jj-96] = readValues[0][jj]; 
151       CorrCoeff1[jj-96] = readValues[1][jj];;
152     }
153     */
154   }
155
156   /* report progress */
157   daqDA_progressReport(10);
158
159
160   /* init some counters */
161   int nevents_physics=0;
162   int nevents_total=0;
163
164   /* read the data files */
165   int n;
166   for(n=1;n<argc;n++){
167    
168     status=monitorSetDataSource( argv[n] );
169     if (status!=0) {
170       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
171       return -1;
172     }
173
174     /* report progress */
175     /* in this example, indexed on the number of files */
176     daqDA_progressReport(10+80*n/argc);
177
178     /* read the file */
179     for(;;){
180       struct eventHeaderStruct *event;
181       eventTypeType eventT;
182
183       /* get next event */
184       status=monitorGetEventDynamic((void **)&event);
185       if(status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
186       if(status!=0) {
187         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
188         return -1;
189       }
190
191       /* retry if got no event */
192       if(event==NULL) {
193         break;
194       }
195       
196       // Initalize raw-data reading and decoding
197       AliRawReader *reader = new AliRawReaderDate((void*)event);
198       reader->Select("ZDC");
199       // --- Reading event header
200       //UInt_t evtype = reader->GetType();
201       //printf("\n\t ZDCEMDda -> ev. type %d\n",evtype);
202       //printf("\t ZDCEMDda -> run # %d\n",reader->GetRunNumber());
203       //
204       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
205         
206
207       /* use event - here, just write event id to result file */
208       eventT=event->eventType;
209       
210       Int_t ich=0, adcMod[48], adcCh[48], sigCode[48], det[48], sec[48];
211       if(eventT==START_OF_DATA){
212                 
213         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
214         else{
215           while(rawStreamZDC->Next()){
216             if(rawStreamZDC->IsChMapping()){
217               adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);
218               adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);
219               sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
220               det[ich] = rawStreamZDC->GetDetectorFromMap(ich);
221               sec[ich] = rawStreamZDC->GetTowerFromMap(ich);
222               ich++;
223             }
224           }
225         }
226         // --------------------------------------------------------
227         // --- Writing ascii data file for the Shuttle preprocessor
228         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
229         for(Int_t i=0; i<ich; i++){
230            fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",i,
231              adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
232            //
233            //printf("ZDCPEDESTALDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
234            //      i,adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
235         }
236         fclose(mapFile4Shuttle);
237       }
238     
239     if(eventT==PHYSICS_EVENT){
240       // --- Reading data header
241       reader->ReadHeader();
242       const AliRawDataHeader* header = reader->GetDataHeader();
243       if(header){
244          UChar_t message = header->GetAttributes();
245          if(message & 0x70){ // DEDICATED EMD RUN
246             //printf("\t STANDALONE_EMD_RUN raw data found\n");
247             continue;
248          }
249          else{
250             printf("\t NO STANDALONE_EMD_RUN raw data found\n");
251             return -1;
252          }
253       }
254       else{
255          printf("\t ATTENTION! No Raw Data Header found!!!\n");
256          return -1;
257       }
258
259       if (!rawStreamZDC->Next()) printf(" \t No raw data found!! ");
260       //
261       // ----- Setting ch. mapping -----
262       for(Int_t jk=0; jk<48; jk++){
263         rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
264         rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
265         rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
266         rawStreamZDC->SetMapDet(jk, det[jk]);
267         rawStreamZDC->SetMapTow(jk, sec[jk]);
268       }
269       //
270       Float_t ZDCRawADC[4], ZDCCorrADC[4], ZDCCorrADCSum[4];
271       for(Int_t g=0; g<4; g++){
272            ZDCCorrADCSum[g] = 0.;
273            ZDCRawADC[g] = 0.;
274       }
275       //
276       while(rawStreamZDC->Next()){
277         if(rawStreamZDC->IsADCDataWord()){
278           Int_t DetIndex=999, PedIndex=999;
279           if(rawStreamZDC->GetSector(0) == 1 || rawStreamZDC->GetSector(0) == 2){
280             DetIndex = rawStreamZDC->GetSector(0)-1;
281             PedIndex = (rawStreamZDC->GetSector(0)+1)+4*rawStreamZDC->GetSector(1);
282           }
283           else if(rawStreamZDC->GetSector(0) == 4 || rawStreamZDC->GetSector(0) == 5){
284             DetIndex = rawStreamZDC->GetSector(0)-2;
285             PedIndex = (rawStreamZDC->GetSector(0)-2)+4*rawStreamZDC->GetSector(1)+24;
286           }
287           //
288           if(rawStreamZDC->GetADCGain() == 1){ //EMD -> LR ADC
289             //
290             ZDCRawADC[DetIndex] += (Float_t) rawStreamZDC->GetADCValue();
291             // Mean pedestal subtraction 
292             Float_t Pedestal = MeanPed[PedIndex];
293             // Pedestal subtraction from correlation with out-of-time signals
294             //Float_t Pedestal = CorrCoeff0[PedIndex]+CorrCoeff1[PedIndex]*MeanPedOOT[PedIndex];
295             //
296             ZDCCorrADC[DetIndex] = (rawStreamZDC->GetADCValue()) - Pedestal;
297             ZDCCorrADCSum[DetIndex] += ZDCCorrADC[DetIndex];
298             //
299             /*printf("\t det %d quad %d res %d pedInd %d ADCCorr %d ZDCCorrADCSum[%d] = %d\n", 
300                rawStreamZDC->GetSector(0),rawStreamZDC->GetSector(1),
301                rawStreamZDC->GetADCGain(),PedIndex,  
302                (Int_t) (rawStreamZDC->GetADCValue() - Pedestal), DetIndex, 
303                (Int_t) ZDCCorrADCSum[DetIndex]);
304             */
305           }   
306         }//IsADCDataWord()
307          //
308        }
309        //
310        nevents_physics++;
311        //
312        for(Int_t j=0; j<4; j++){
313           histoEMDRaw[j]->Fill(ZDCRawADC[j]);
314           histoEMDCorr[j]->Fill(ZDCCorrADCSum[j]);
315        }
316     }
317     
318     nevents_total++;
319
320     /* free resources */
321     free(event);
322     
323     /* exit when last event received, no need to wait for TERM signal */
324     if (eventT==END_OF_RUN) {
325       printf("EOR event detected\n");
326       break;
327     }
328
329    }
330   }
331     
332   /* Analysis of the histograms */
333   //
334   FILE *fileShuttle;
335   fileShuttle = fopen(EMDDATA_FILE,"w");
336   //
337   Int_t BinMax[4];
338   Float_t YMax[4];
339   Int_t NBinsx[4];
340   Float_t MeanFitVal[4];
341   TF1 *fitfun[4];
342   for(Int_t k=0; k<4; k++){
343      BinMax[k] = histoEMDCorr[k]->GetMaximumBin();
344      YMax[k] = (histoEMDCorr[k]->GetXaxis())->GetXmax();
345      NBinsx[k] = (histoEMDCorr[k]->GetXaxis())->GetNbins();
346 //     printf("\n\t Det%d -> BinMax = %d, ChXMax = %f\n", k+1, BinMax[k], BinMax[k]*YMax[k]/NBinsx[k]);
347      histoEMDCorr[k]->Fit("gaus","Q","",BinMax[k]*YMax[k]/NBinsx[k]*0.7,BinMax[k]*YMax[k]/NBinsx[k]*1.25);
348      fitfun[k] = histoEMDCorr[k]->GetFunction("gaus");
349      MeanFitVal[k] = (Float_t) (fitfun[k]->GetParameter(1));
350      //printf("\n\t Mean Value from gaussian fit = %f\n", MeanFitVal[k]);
351   }
352   //
353   Float_t CalibCoeff[6];     
354   Float_t icoeff[5];
355   //
356   for(Int_t j=0; j<10; j++){
357      if(j<4){
358        CalibCoeff[j] = MeanFitVal[j];
359        fprintf(fileShuttle,"\t%f\n",CalibCoeff[j]);
360      }
361      // ZEM calib. coeff. = 1
362      else if(j==4 || j==5){
363        CalibCoeff[j] = 1.; 
364        fprintf(fileShuttle,"\t%f\n",CalibCoeff[j]);
365      }
366      // Note -> For the moment the inter-calibration
367      //      coefficients are set to 1 
368      else if(j>5){
369        for(Int_t k=0; k<5; k++){  
370          icoeff[k] = 1.;
371          fprintf(fileShuttle,"\t%f",icoeff[k]);
372          if(k==4) fprintf(fileShuttle,"\n");
373        }
374      }
375   }
376   //                                                   
377   fclose(fileShuttle);
378   
379   for(Int_t ij=0; ij<4; ij++){
380     delete histoEMDRaw[ij];
381     delete histoEMDCorr[ij];
382   }
383   
384   //delete minuitFit;
385   TVirtualFitter::SetFitter(0);
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(MAPDATA_FILE, MAPDATA_FILE);
398   if(status){
399     printf("Failed to export file : %d\n",status);
400     return -1;
401   }
402   //
403   status = daqDA_FES_storeFile(EMDDATA_FILE, EMDDATA_FILE);
404   if(status){
405     printf("Failed to export file : %d\n",status);
406     return -1;
407   }
408
409   /* report progress */
410   daqDA_progressReport(100);
411
412   return status;
413 }