Alghorithm to calculate intercalib. coeff. implemented
[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 CALIBRATION_EMD 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 ENCALIBDATA_FILE   "ZDCEnergyCalib.dat"
27 #define TOWCALIBDATA_FILE  "ZDCTowerCalib.dat"
28
29 #include <stdio.h>
30 #include <Riostream.h>
31 #include <Riostream.h>
32
33 // DATE
34 #include <daqDA.h>
35 #include <event.h>
36 #include <monitor.h>
37
38 //ROOT
39 #include <TROOT.h>
40 #include <TPluginManager.h>
41 #include <TH1F.h>
42 #include <TH2F.h>
43 #include <TProfile.h>
44 #include <TF1.h>
45 #include <TFile.h>
46 #include <TFitter.h>
47 #include "TMinuitMinimizer.h"
48
49 //AliRoot
50 #include <AliRawReaderDate.h>
51 #include <AliRawEventHeaderBase.h>
52 #include <AliZDCRawStream.h>
53
54
55 /* Main routine
56       Arguments: 
57       1- monitoring data source
58 */
59 int main(int argc, char **argv) {
60   
61
62   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
63                                         "*",
64                                         "TStreamerInfo",
65                                         "RIO",
66                                         "TStreamerInfo()"); 
67
68   TMinuitMinimizer m; 
69   gROOT->GetPluginManager()->AddHandler("ROOT::Math::Minimizer", "Minuit","TMinuitMinimizer",
70       "Minuit", "TMinuitMinimizer(const char *)");
71   TVirtualFitter::SetDefaultFitter("Minuit");
72
73   int status = 0;
74   // No. of ZDC cabled ch.
75   int const kNChannels = 24;
76   int const kNScChannels = 32;
77
78   /* log start of process */
79   printf("ZDC EMD program started\n");  
80
81   /* check that we got some arguments = list of files */
82   if (argc<2) {
83     printf("Wrong number of arguments\n");
84     return -1;
85   }
86   
87   // --- Preparing histos for EM dissociation spectra
88   //
89   TH1F* histoEMDRaw[4];
90   TH1F* histoEMDCorr[4];
91   //
92   char namhistr[50], namhistc[50];
93   for(Int_t i=0; i<4; i++) {
94      if(i==0){
95        sprintf(namhistr,"ZN%d-EMDRaw",i+1);
96        sprintf(namhistc,"ZN%d-EMDCorr",i+1);
97      }
98      else if(i==1){
99        sprintf(namhistr,"ZP%d-EMDRaw",i);
100        sprintf(namhistc,"ZP%d-EMDCorr",i);
101      }
102      else if(i==2){
103        sprintf(namhistr,"ZN%d-EMDRaw",i);
104        sprintf(namhistc,"ZN%d-EMDCorr",i);
105      }
106      else if(i==3){
107        sprintf(namhistr,"ZP%d-EMDRaw",i-1);
108        sprintf(namhistc,"ZP%d-EMDCorr",i-1);
109      }
110      histoEMDRaw[i] = new TH1F(namhistr,namhistr,100,0.,4000.);
111      histoEMDCorr[i] = new TH1F(namhistc,namhistc,100,0.,4000.);
112   }
113
114    // --- Preparing histos for tower inter-calibration
115   //
116   TH1F* histZNCtow[4]; TH1F* histZPCtow[4];
117   TH1F* histZNAtow[4]; TH1F* histZPAtow[4];
118   //
119   char namhistznc[50], namhistzpc[50];
120   char namhistzna[50], namhistzpa[50];
121   for(Int_t i=0; i<4; i++) {
122      sprintf(namhistznc,"ZNC-tow%d",i+1);
123      sprintf(namhistzpc,"ZPC-tow%d",i+1);
124      sprintf(namhistzna,"ZNA-tow%d",i+1);
125      sprintf(namhistzpa,"ZPA-tow%d",i+1);
126      //
127      histZNCtow[i] = new TH1F(namhistznc,namhistznc,100,0.,4000.);
128      histZPCtow[i] = new TH1F(namhistzpc,namhistzpc,100,0.,4000.);
129      histZNAtow[i] = new TH1F(namhistzna,namhistzna,100,0.,4000.);
130      histZPAtow[i] = new TH1F(namhistzpa,namhistzpa,100,0.,4000.);
131   }
132   
133   /* open result file */
134   FILE *fp=NULL;
135   fp=fopen("./result.txt","a");
136   if (fp==NULL) {
137     printf("Failed to open file\n");
138     return -1;
139   }
140   
141   FILE *mapFile4Shuttle;
142
143   // *** To analyze EMD events you MUST have a pedestal data file!!!
144   // *** -> check if a pedestal run has been analyzed
145   int read = 0;
146   read = daqDA_DB_getFile(PEDDATA_FILE,PEDDATA_FILE);
147   if(read){
148     printf("\t ERROR!!! ZDCPedestal.dat file NOT FOUND in DAQ db!!!\n");
149     return -1;
150   }
151   else printf("\t ZDCPedestal.dat file retrieved from DAQ db\n");
152   
153   FILE *filePed = fopen(PEDDATA_FILE,"r");
154   if (filePed==NULL) {
155     printf("\t ERROR!!! Can't open ZDCPedestal.dat file!!!\n");
156     return -1;
157   }
158
159   // 144 = 48 in-time + 48 out-of-time + 48 correlations
160   Float_t readValues[2][3*2*kNChannels];
161   Float_t MeanPed[2*kNChannels];
162   Float_t CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
163   // ***************************************************
164   //   Unless we have a narrow correlation to fit we
165   //    don't fit and store in-time vs. out-of-time
166   //    histograms -> mean pedstal subtracted!!!!!!
167   // ***************************************************
168   //
169   for(int jj=0; jj<6*kNChannels; jj++){
170     for(int ii=0; ii<2; ii++){
171        fscanf(filePed,"%f",&readValues[ii][jj]);
172     }
173     if(jj<kNChannels && jj<2*kNChannels){
174       MeanPed[jj] = readValues[0][jj];
175       //printf("\t MeanPedhg[%d] = %1.1f\n",jj, MeanPedhg[jj]);
176     }
177     else if(jj>2*kNChannels && jj>4*kNChannels){
178       CorrCoeff0[jj-4*kNChannels] = readValues[0][jj]; 
179       CorrCoeff1[jj-4*kNChannels] = readValues[1][jj];;
180     }
181   }
182
183   /* report progress */
184   daqDA_progressReport(10);
185
186
187   /* init some counters */
188   int nevents_physics=0;
189   int nevents_total=0;
190
191   struct eventHeaderStruct *event;
192   eventTypeType eventT;
193
194   /* read the data files */
195   int n;
196   for(n=1;n<argc;n++){
197    
198     status=monitorSetDataSource( argv[n] );
199     if (status!=0) {
200       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
201       return -1;
202     }
203
204     /* report progress */
205     /* in this example, indexed on the number of files */
206     daqDA_progressReport(10+80*n/argc);
207
208     /* read the file */
209     for(;;){
210
211       /* get next event */
212       status=monitorGetEventDynamic((void **)&event);
213       if(status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
214       if(status!=0) {
215         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
216         return -1;
217       }
218
219       /* retry if got no event */
220       if(event==NULL) {
221         break;
222       }
223       
224       // Initalize raw-data reading and decoding
225       AliRawReader *reader = new AliRawReaderDate((void*)event);
226       reader->Select("ZDC");
227       // --- Reading event header
228       //UInt_t evtype = reader->GetType();
229       //printf("\n\t ZDCEMDda -> ev. type %d\n",evtype);
230       //printf("\t ZDCEMDda -> run # %d\n",reader->GetRunNumber());
231       //
232       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
233         
234
235       /* use event - here, just write event id to result file */
236       eventT=event->eventType;
237       
238       Int_t ich=0;
239       Int_t adcMod[2*kNChannels], adcCh[2*kNChannels], sigCode[2*kNChannels];
240       Int_t det[2*kNChannels], sec[2*kNChannels];
241       for(Int_t y=0; y<2*kNChannels; y++){
242         adcMod[y]=adcCh[y]=sigCode[y]=det[y]=sec[y]=0;
243       }
244       
245       Int_t iScCh=0;
246       Int_t scMod[kNScChannels], scCh[kNScChannels], scSigCode[kNScChannels];
247       Int_t scDet[kNScChannels], scSec[kNScChannels];
248       for(Int_t y=0; y<kNScChannels; y++){
249         scMod[y]=scCh[y]=scSigCode[y]=scDet[y]=scSec[y]=0;
250       }
251       //
252       Int_t modNum=-1, modType=-1;
253       
254       if(eventT==START_OF_DATA){
255                 
256         rawStreamZDC->SetSODReading(kTRUE);
257         
258         // --------------------------------------------------------
259         // --- Writing ascii data file for the Shuttle preprocessor
260         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
261         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
262         else{
263           while((rawStreamZDC->Next())){
264             if(rawStreamZDC->IsHeaderMapping()){ // mapping header
265                modNum = rawStreamZDC->GetADCModule();
266                modType = rawStreamZDC->GetModType();
267             }
268             if(rawStreamZDC->IsChMapping()){ 
269               if(modType==1){ // ADC mapping ----------------------
270                 adcMod[ich]  = rawStreamZDC->GetADCModFromMap(ich);
271                 adcCh[ich]   = rawStreamZDC->GetADCChFromMap(ich);
272                 sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
273                 det[ich]     = rawStreamZDC->GetDetectorFromMap(ich);
274                 sec[ich]     = rawStreamZDC->GetTowerFromMap(ich);
275                 //
276                 fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",
277                   ich,adcMod[ich],adcCh[ich],sigCode[ich],det[ich],sec[ich]);
278                 //
279                 //printf("  Mapping in DA -> %d ADC: mod %d ch %d, code %d det %d, sec %d\n",
280                 //  ich,adcMod[ich],adcCh[ich],sigCode[ich],det[ich],sec[ich]);
281                 //
282                 ich++;
283               }
284               else if(modType==2){ //VME scaler mapping --------------------
285                 scMod[iScCh]     = rawStreamZDC->GetScalerModFromMap(iScCh);
286                 scCh[iScCh]      = rawStreamZDC->GetScalerChFromMap(iScCh);
287                 scSigCode[iScCh] = rawStreamZDC->GetScalerSignFromMap(iScCh);
288                 scDet[iScCh]     = rawStreamZDC->GetScDetectorFromMap(iScCh);
289                 scSec[iScCh]    = rawStreamZDC->GetScTowerFromMap(iScCh);
290                 //
291                 fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",
292                   iScCh,scMod[iScCh],scCh[iScCh],scSigCode[iScCh],scDet[iScCh],scSec[iScCh]);
293                 //
294                 //printf("  Mapping in DA -> %d Scaler: mod %d ch %d, code %d det %d, sec %d\n",
295                 //  iScCh,scMod[iScCh],scCh[iScCh],scSigCode[iScCh],scDet[iScCh],scSec[iScCh]);
296                 //
297                 iScCh++;
298               }
299             }
300           }
301         }
302         fclose(mapFile4Shuttle);
303       }// SOD event
304     
305     if(eventT==PHYSICS_EVENT){
306       // --- Reading data header
307       reader->ReadHeader();
308       const AliRawDataHeader* header = reader->GetDataHeader();
309       if(header){
310          UChar_t message = header->GetAttributes();
311          if(message & 0x70){ // DEDICATED EMD RUN
312             //printf("\t STANDALONE_EMD_RUN raw data found\n");
313             continue;
314          }
315          else{
316             printf("\t NO STANDALONE_EMD_RUN raw data found\n");
317             return -1;
318          }
319       }
320       else{
321          printf("\t ATTENTION! No Raw Data Header found!!!\n");
322          return -1;
323       }
324       
325       rawStreamZDC->SetSODReading(kTRUE);
326
327       if (!rawStreamZDC->Next()) printf(" \t No raw data found!! ");
328       //
329       // ----- Setting ch. mapping -----
330       for(Int_t jk=0; jk<2*kNChannels; jk++){
331         rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
332         rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
333         rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
334         rawStreamZDC->SetMapDet(jk, det[jk]);
335         rawStreamZDC->SetMapTow(jk, sec[jk]);
336       }
337       //
338       Float_t ZDCRawADC[4], ZDCCorrADC[4], ZDCCorrADCSum[4];
339       for(Int_t g=0; g<4; g++){
340            ZDCCorrADCSum[g] = 0.;
341            ZDCRawADC[g] = 0.;
342       }
343       //
344       while(rawStreamZDC->Next()){
345         Int_t det = rawStreamZDC->GetSector(0);
346         Int_t quad = rawStreamZDC->GetSector(1);
347         
348         if(rawStreamZDC->IsADCDataWord() && !(rawStreamZDC->IsUnderflow())
349              && !(rawStreamZDC->IsOverflow()) && det!=-1
350              && (rawStreamZDC->GetADCGain() == 1)){ // Selecting LOW RES ch.s
351
352           //printf("  IsADCWord %d, IsUnderflow %d, IsOverflow %d\n",
353           //  rawStreamZDC->IsADCDataWord(),rawStreamZDC->IsUnderflow(),rawStreamZDC->IsOverflow());
354           
355           // Taking LOW RES channels -> ch.+kNChannels !!!!
356           Int_t DetIndex=999, PedIndex=999;
357           // Not ZEM not PMRef && low gain chain
358           if((det!=3) && (quad!=5) && (rawStreamZDC->GetADCGain()==1)){ 
359             if(det == 1){
360               DetIndex = det-1;
361               PedIndex = quad+kNChannels; 
362             }
363             else if(det==2){
364               DetIndex = det-1;
365               PedIndex = quad+5+kNChannels;
366             }
367             else if(det == 4){
368               DetIndex = det-2;
369               PedIndex = quad+12+kNChannels;
370             }
371             else if(det == 5){
372               DetIndex = det-2;
373               PedIndex = quad+17+kNChannels;
374             }
375             // Mean pedestal subtraction 
376             Float_t Pedestal = MeanPed[PedIndex];
377             // Pedestal subtraction from correlation with out-of-time signals
378             //Float_t Pedestal = CorrCoeff0[PedIndex]+CorrCoeff1[PedIndex]*MeanPedOOT[PedIndex];
379             //
380             if(DetIndex!=999 || PedIndex!=999){ 
381               //
382               ZDCRawADC[DetIndex] += (Float_t) rawStreamZDC->GetADCValue();
383               //
384               //
385               ZDCCorrADC[DetIndex] = (rawStreamZDC->GetADCValue()) - Pedestal;
386               ZDCCorrADCSum[DetIndex] += ZDCCorrADC[DetIndex];
387               //
388               /*printf("\t det %d quad %d res %d pedInd %d detInd %d:"
389                  "ADCCorr = %d, ZDCCorrADCSum = %d\n", 
390                  det,quad,rawStreamZDC->GetADCGain(),PedIndex,DetIndex, 
391                  (Int_t) ZDCCorrADC[DetIndex],(Int_t) ZDCCorrADCSum[DetIndex]);
392               */              
393             }
394             // Not common nor ref. PM && low gain chain
395             if((quad!=0) && (rawStreamZDC->GetADCGain()==1)){
396               Float_t corrADCval = (rawStreamZDC->GetADCValue()) - Pedestal;
397               if(det==1)      histZNCtow[quad]->Fill(corrADCval);
398               else if(det==2) histZPCtow[quad]->Fill(corrADCval);
399               else if(det==4) histZNAtow[quad]->Fill(corrADCval);
400               else if(det==5) histZPAtow[quad]->Fill(corrADCval);
401             }
402
403             if(DetIndex==999 || PedIndex==999) 
404                 printf(" WARNING! Detector a/o pedestal index are WRONG!!!\n");
405  
406           }
407         }//IsADCDataWord()
408         
409        }
410        //
411        nevents_physics++;
412        //
413        for(Int_t j=0; j<4; j++){
414           histoEMDRaw[j]->Fill(ZDCRawADC[j]);
415           histoEMDCorr[j]->Fill(ZDCCorrADCSum[j]);
416        }
417     }//(if PHYSICS_EVENT)
418     
419     nevents_total++;
420
421    }
422
423    /* free resources */
424    free(event);
425   }
426     
427   /* Analysis of the histograms */
428   //
429   FILE *fileShuttle1 = fopen(ENCALIBDATA_FILE,"w");
430   //
431   Int_t BinMax[4];
432   Float_t YMax[4];
433   Int_t NBinsx[4];
434   Float_t MeanFitVal[4];
435   TF1 *fitfun[4];
436   for(Int_t k=0; k<4; k++){
437      if(histoEMDCorr[k]->GetEntries() == 0){
438        printf("\n WARNING! Empty histos -> ending DA WITHOUT writing output\n\n");
439        return -1;
440      } 
441      BinMax[k] = histoEMDCorr[k]->GetMaximumBin();
442      YMax[k] = (histoEMDCorr[k]->GetXaxis())->GetXmax();
443      NBinsx[k] = (histoEMDCorr[k]->GetXaxis())->GetNbins();
444 //     printf("\n\t Det%d -> BinMax = %d, ChXMax = %f\n", k+1, BinMax[k], BinMax[k]*YMax[k]/NBinsx[k]);
445      histoEMDCorr[k]->Fit("gaus","Q","",BinMax[k]*YMax[k]/NBinsx[k]*0.7,BinMax[k]*YMax[k]/NBinsx[k]*1.25);
446      fitfun[k] = histoEMDCorr[k]->GetFunction("gaus");
447      MeanFitVal[k] = (Float_t) (fitfun[k]->GetParameter(1));
448      //printf("\n\t Mean Value from gaussian fit = %f\n", MeanFitVal[k]);
449   }
450   //
451   Float_t CalibCoeff[6];     
452   //Float_t icoeff[5];
453   //
454   for(Int_t j=0; j<6; j++){
455      if(j<4){
456        CalibCoeff[j] = MeanFitVal[j];
457        fprintf(fileShuttle1,"\t%f\n",CalibCoeff[j]);
458      }
459      // ZEM energy calib. coeff. = 1
460      else if(j==4 || j==5){
461        CalibCoeff[j] = 1.; 
462        fprintf(fileShuttle1,"\t%f\n",CalibCoeff[j]);
463      }
464   }
465   fclose(fileShuttle1);
466  
467   FILE *fileShuttle2 = fopen(TOWCALIBDATA_FILE,"w");
468   Float_t meanvalznc[4], meanvalzpc[4], meanvalzna[4], meanvalzpa[4];
469   for(Int_t j=0; j<4; j++){
470      if(histZNCtow[j]->GetEntries() == 0){
471        printf("\n WARNING! Empty histos -> ending DA WITHOUT writing output\n\n");
472        return -1;
473      } 
474      meanvalznc[j] = histZNCtow[j]->GetMean();
475      meanvalzpc[j] = histZPCtow[j]->GetMean();
476      meanvalzna[j] = histZNAtow[j]->GetMean();
477      meanvalzpa[j] = histZPAtow[j]->GetMean();
478      
479      // Note -> For the moment the inter-calibration coeff. are set to 1 
480      /*for(Int_t k=0; k<4; k++){  
481        icoeff[k] = 1.;
482        fprintf(fileShuttle2,"\t%f",icoeff[k]);
483        if(k==4) fprintf(fileShuttle2,"\n");
484      }*/
485   }
486   if(meanvalznc[1]!=0 && meanvalznc[2]!=0 && meanvalznc[3]!=0 && 
487      meanvalzpc[1]!=0 && meanvalzpc[2]!=0 && meanvalzpc[3]!=0 &&
488      meanvalzna[1]!=0 && meanvalzna[2]!=0 && meanvalzna[3]!=0 &&
489      meanvalzpa[1]!=0 && meanvalzpa[2]!=0 && meanvalzpa[3]!=0){
490     fprintf(fileShuttle2,"\t%f\t%f\t%f\t%f\n",
491         1.0,meanvalznc[0]/meanvalznc[1],meanvalznc[0]/meanvalznc[2],meanvalznc[0]/meanvalznc[3]);
492     fprintf(fileShuttle2,"\t%f\t%f\t%f\t%f\n",
493         1.0,meanvalzpc[0]/meanvalzpc[1],meanvalzpc[0]/meanvalzpc[2],meanvalzpc[0]/meanvalzpc[3]);
494     fprintf(fileShuttle2,"\t%f\t%f\t%f\t%f\n",
495         1.0,meanvalzna[0]/meanvalzna[1],meanvalzpc[0]/meanvalzna[2],meanvalzpc[0]/meanvalzna[3]);
496     fprintf(fileShuttle2,"\t%f\t%f\t%f\t%f\n",
497         1.0,meanvalzpa[0]/meanvalzpa[1],meanvalzpc[0]/meanvalzpa[2],meanvalzpc[0]/meanvalzpa[3]);
498   }
499   else{
500     printf("\n Tower intercalib. coeff. CAN'T be calculated (some mean values are ZERO)!!!\n\n");
501     return -1;
502   }
503   fclose(fileShuttle2);
504   
505   for(Int_t ij=0; ij<4; ij++){
506     delete histoEMDRaw[ij];
507     delete histoEMDCorr[ij];
508   }
509   
510   //delete minuitFit;
511   TVirtualFitter::SetFitter(0);
512
513   /* write report */
514   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
515
516   /* close result file */
517   fclose(fp);
518   
519   /* report progress */
520   daqDA_progressReport(90);
521
522   /* store the result file on FES */
523   status = daqDA_FES_storeFile(MAPDATA_FILE, "MAPPING");
524   if(status){
525     printf("Failed to export file : %d\n",status);
526     return -1;
527   }
528   //
529   status = daqDA_FES_storeFile(ENCALIBDATA_FILE, "EMDENERGYCALIB");
530   if(status){
531     printf("Failed to export file : %d\n",status);
532     return -1;
533   }
534   //
535   status = daqDA_FES_storeFile(TOWCALIBDATA_FILE, "EMDTOWERCALIB");
536   if(status){
537     printf("Failed to export file : %d\n",status);
538     return -1;
539   }
540
541   /* report progress */
542   daqDA_progressReport(100);
543
544   return status;
545 }