Changes for report #71958: ZDC: changes needs code update in STEER
[u/mrichter/AliRoot.git] / ZDC / ZDCCALIBMBda.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_MB runs
13
14 Contact: Chiara.Oppedisano@to.infn.it
15 Link: 
16 Run Type: CALIBRATION_MB_RUN
17 DA Type: LDC
18 Number of events needed: at least 10^4
19 Input Files: ZDCPedestal.dat
20 Output Files: ZDCMBCalib.root, ZDCChMapping.dat
21 Trigger Types Used: Standalone Trigger
22
23 */
24 #define PEDDATA_FILE  "ZDCPedestal.dat"
25 #define MAPDATA_FILE  "ZDCChMapping.dat"
26 #define MBCALIBDATA_FILE   "ZDCMBCalib.root"
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 <TROOT.h>
39 #include <TPluginManager.h>
40 #include <TH1F.h>
41 #include <TH2F.h>
42 #include <TProfile.h>
43 #include <TF1.h>
44 #include <TFile.h>
45 #include <TFitter.h>
46 #include "TMinuitMinimizer.h"
47
48 //AliRoot
49 #include <AliRawReaderDate.h>
50 #include <AliRawEventHeaderBase.h>
51 #include <AliZDCRawStream.h>
52
53
54 /* Main routine
55       Arguments: 
56       1- monitoring data source
57 */
58 int main(int argc, char **argv) {
59   
60
61   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
62                                         "*",
63                                         "TStreamerInfo",
64                                         "RIO",
65                                         "TStreamerInfo()"); 
66
67   TMinuitMinimizer m; 
68   gROOT->GetPluginManager()->AddHandler("ROOT::Math::Minimizer", "Minuit","TMinuitMinimizer",
69       "Minuit", "TMinuitMinimizer(const char *)");
70   TVirtualFitter::SetDefaultFitter("Minuit");
71
72   int status = 0;
73   int const kNModules = 10;
74   int const kNChannels = 24;
75   int const kNScChannels = 32;
76   Int_t kFirstADCGeo=0, kLastADCGeo=3;
77       
78   Int_t iMod=-1;
79   Int_t modGeo[kNModules], modType[kNModules],modNCh[kNModules];
80   for(Int_t kl=0; kl<kNModules; kl++){
81      modGeo[kl]=modType[kl]=modNCh[kl]=0;
82   }
83   
84   Int_t ich=0;
85   Int_t adcMod[2*kNChannels], adcCh[2*kNChannels], sigCode[2*kNChannels];
86   Int_t det[2*kNChannels], sec[2*kNChannels];
87   for(Int_t y=0; y<2*kNChannels; y++){
88     adcMod[y]=adcCh[y]=sigCode[y]=det[y]=sec[y]=0;
89   }
90   
91   Int_t iScCh=0;
92   Int_t scMod[kNScChannels], scCh[kNScChannels], scSigCode[kNScChannels];
93   Int_t scDet[kNScChannels], scSec[kNScChannels];
94   for(Int_t y=0; y<kNScChannels; y++){
95     scMod[y]=scCh[y]=scSigCode[y]=scDet[y]=scSec[y]=0;
96   }
97       
98   Int_t itdcCh=0;
99   Int_t tdcMod[kNScChannels], tdcCh[kNScChannels], tdcSigCode[kNScChannels];
100   Int_t tdcDet[kNScChannels], tdcSec[kNScChannels];
101   for(Int_t y=0; y<kNScChannels; y++){
102     tdcMod[y]=tdcCh[y]=tdcSigCode[y]=tdcDet[y]=tdcSec[y]=-1;
103   }
104
105   /* log start of process */
106   printf("ZDC EMD program started\n");  
107
108   /* check that we got some arguments = list of files */
109   if (argc<2) {
110     printf("Wrong number of arguments\n");
111     return -1;
112   }
113   
114   // --- Preparing histos for EM dissociation spectra
115   //
116   TH2F * hZDCvsZEM  = new TH2F("hZDCvsZEM","EZDCTot vs. EZEM",100,0.,5.,100,0.,800.);
117   TH2F * hZDCCvsZEM = new TH2F("hZDCCvsZEM","EZDCC vs. EZEM",100,0.,5.,100,0.,400.);
118   TH2F * hZDCAvsZEM = new TH2F("hZDCAvsZEM","EZDCA vs. EZEM",100,0.,5.,100,0.,400.);
119   
120   /* open result file */
121   FILE *fp=NULL;
122   fp=fopen("./result.txt","a");
123   if (fp==NULL) {
124     printf("Failed to open file\n");
125     return -1;
126   }
127   
128   FILE *mapFile4Shuttle;
129
130   // *** To analyze CALIBRATION_MB events a pedestal data file is needed!!!
131   // *** -> check if a pedestal run has been analyzed
132   int read = 0;
133   read = daqDA_DB_getFile(PEDDATA_FILE,PEDDATA_FILE);
134   if(read){
135     printf("\t ERROR!!! ZDCPedestal.dat file NOT FOUND in DAQ db!!!\n");
136     return -1;
137   }
138   else printf("\t ZDCPedestal.dat file retrieved from DAQ db\n");
139   
140   FILE *filePed = fopen(PEDDATA_FILE,"r");
141   if (filePed==NULL) {
142     printf("\t ERROR!!! Can't open ZDCPedestal.dat file!!!\n");
143     return -1;
144   }
145
146   // 144 = 48 in-time + 48 out-of-time + 48 correlations
147   Float_t readValues[2][6*kNChannels];
148   Float_t MeanPedhg[kNChannels], MeanPedlg[kNChannels];
149   Float_t CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
150   // ***************************************************
151   //   Unless we have a narrow correlation to fit we
152   //    don't fit and store in-time vs. out-of-time
153   //    histograms -> mean pedstal subtracted!!!!!!
154   // ***************************************************
155   //
156   for(int jj=0; jj<6*kNChannels; jj++){
157     for(int ii=0; ii<2; ii++){
158        fscanf(filePed,"%f",&readValues[ii][jj]);
159     }
160     if(jj<kNChannels){
161       MeanPedhg[jj] = readValues[0][jj];
162       //printf("\t MeanPedhg[%d] = %1.1f\n",jj, MeanPedhg[jj]);
163     }
164     else if(jj>=kNChannels && jj<2*kNChannels){
165       MeanPedlg[jj-kNChannels] = readValues[0][jj];
166       //printf("\t MeanPedlg[%d] = %1.1f\n",jj-kNChannels, MeanPedlg[jj-kNChannels]);
167     }
168     else if(jj>4*kNChannels){
169       CorrCoeff0[jj-4*kNChannels] = readValues[0][jj]; 
170       CorrCoeff1[jj-4*kNChannels] = readValues[1][jj];;
171     }
172   }
173
174   /* report progress */
175   daqDA_progressReport(10);
176
177
178   /* init some counters */
179   int nevents_physics=0;
180   int nevents_total=0;
181
182   struct eventHeaderStruct *event;
183   eventTypeType eventT;
184
185   /* read the data files */
186   int n;
187   for(n=1;n<argc;n++){
188    
189     status=monitorSetDataSource( argv[n] );
190     if (status!=0) {
191       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
192       return -1;
193     }
194
195     /* report progress */
196     /* in this example, indexed on the number of files */
197     daqDA_progressReport(10+80*n/argc);
198
199     /* read the file */
200     for(;;){
201
202       /* get next event */
203       status=monitorGetEventDynamic((void **)&event);
204       if(status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
205       if(status!=0) {
206         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
207         return -1;
208       }
209
210       /* retry if got no event */
211       if(event==NULL) {
212         break;
213       }
214       
215       // Initalize raw-data reading and decoding
216       AliRawReader *reader = new AliRawReaderDate((void*)event);
217       reader->Select("ZDC");
218       // --- Reading event header
219       //UInt_t evtype = reader->GetType();
220       //printf("\n\t ZDCEMDda -> ev. type %d\n",evtype);
221       //printf("\t ZDCEMDda -> run # %d\n",reader->GetRunNumber());
222       //
223       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
224         
225
226       /* use event - here, just write event id to result file */
227       eventT=event->eventType;
228       
229       if(eventT==START_OF_DATA){
230                         
231         rawStreamZDC->SetSODReading(kTRUE);
232         
233         // --------------------------------------------------------
234         // --- Writing ascii data file for the Shuttle preprocessor
235         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
236         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
237         else{
238           while((rawStreamZDC->Next())){
239             if(rawStreamZDC->IsHeaderMapping()){ // mapping header
240                iMod++;
241                modGeo[iMod]  = rawStreamZDC->GetADCModule();
242                modType[iMod] = rawStreamZDC->GetModType();
243                modNCh[iMod]  = rawStreamZDC->GetADCNChannels();
244             }
245             if(rawStreamZDC->IsChMapping()){ 
246               if(modType[iMod]==1){ // ADC mapping ----------------------
247                 adcMod[ich]  = rawStreamZDC->GetADCModFromMap(ich);
248                 adcCh[ich]   = rawStreamZDC->GetADCChFromMap(ich);
249                 sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
250                 det[ich]     = rawStreamZDC->GetDetectorFromMap(ich);
251                 sec[ich]     = rawStreamZDC->GetTowerFromMap(ich);
252                 ich++;
253               }
254               else if(modType[iMod]==2){ //VME scaler mapping --------------------
255                 scMod[iScCh]     = rawStreamZDC->GetScalerModFromMap(iScCh);
256                 scCh[iScCh]      = rawStreamZDC->GetScalerChFromMap(iScCh);
257                 scSigCode[iScCh] = rawStreamZDC->GetScalerSignFromMap(iScCh);
258                 scDet[iScCh]     = rawStreamZDC->GetScDetectorFromMap(iScCh);
259                 scSec[iScCh]     = rawStreamZDC->GetScTowerFromMap(iScCh);
260                 iScCh++;
261               }
262               else if(modType[iMod]==6 && modGeo[iMod]==4){ // ZDC TDC mapping --------------------
263                 tdcMod[itdcCh]     = rawStreamZDC->GetTDCModFromMap(itdcCh);
264                 tdcCh[itdcCh]      = rawStreamZDC->GetTDCChFromMap(itdcCh);
265                 tdcSigCode[itdcCh] = rawStreamZDC->GetTDCSignFromMap(itdcCh);
266                 itdcCh++;
267               }
268             }
269           }
270           // Writing data on output FXS file
271           for(Int_t is=0; is<2*kNChannels; is++){
272              fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",
273                is,adcMod[is],adcCh[is],sigCode[is],det[is],sec[is]);
274              //printf("  CalibMB DA -> %d ADC: mod %d ch %d, code %d det %d, sec %d\n",
275              //  is,adcMod[is],adcCh[is],sigCode[is],det[is],sec[is]);
276           }
277           for(Int_t is=0; is<kNScChannels; is++){
278              fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",
279                is,scMod[is],scCh[is],scSigCode[is],scDet[is],scSec[is]);
280              //printf("  CalibMB DA -> %d Scaler: mod %d ch %d, code %d det %d, sec %d\n",
281              //  is,scMod[is],scCh[is],scSigCode[is],scDet[is],scSec[is]);
282           }
283           for(Int_t is=0; is<kNScChannels; is++){
284              fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\n",
285                is,tdcMod[is],tdcCh[is],tdcSigCode[is]);
286              //if(tdcMod[is]!=-1) printf("  Mapping DA -> %d TDC: mod %d ch %d, code %d\n",
287              //  is,tdcMod[is],tdcCh[is],tdcSigCode[is]);
288           }
289           for(Int_t is=0; is<kNModules; is++){
290              fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\n",
291              modGeo[is],modType[is],modNCh[is]);
292              //printf("  CalibMB DA -> Module mapping: geo %d type %d #ch %d\n",
293              //  modGeo[is],modType[is],modNCh[is]);
294           }
295           
296         }
297         fclose(mapFile4Shuttle);
298       }// SOD event
299     
300     if(eventT==PHYSICS_EVENT){
301       // --- Reading data header
302       reader->ReadHeader();
303       const AliRawDataHeader* header = reader->GetDataHeader();
304       if(header){
305          UChar_t message = header->GetAttributes();
306          if((message & 0xf0) == 0x60){ // DEDICATED CALIBRATION MB RUN
307             //printf("\t STANDALONE_CALIBRAION_MB_RUN raw data found\n");
308             continue;
309          }
310          else{
311             printf("\t NO STANDALONE_MB_RUN raw data found\n");
312             return -1;
313          }
314       }
315       else{
316          printf("\t ATTENTION! No Raw Data Header found!!!\n");
317          return -1;
318       }
319       
320       rawStreamZDC->SetSODReading(kTRUE);
321
322       if (!rawStreamZDC->Next()) printf(" \t No raw data found!! ");
323       //
324       // ----- Setting ch. mapping -----
325       for(Int_t jk=0; jk<2*kNChannels; jk++){
326         rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
327         rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
328         rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
329         rawStreamZDC->SetMapDet(jk, det[jk]);
330         rawStreamZDC->SetMapTow(jk, sec[jk]);
331       }
332       //
333       Float_t rawZDC=0., rawZEM=0.;
334       Float_t corrZDC=0., corrZEM=0.;
335       Float_t corrZDCTot=0., corrZDCC=0., corrZDCA=0., corrZEMTot=0.;
336       //
337       while(rawStreamZDC->Next()){
338         Int_t detector = rawStreamZDC->GetSector(0);
339         Int_t quad = rawStreamZDC->GetSector(1);
340         Int_t index=-1;
341         
342         if( (rawStreamZDC->IsADCDataWord()) && !(rawStreamZDC->IsUnderflow())
343              && !(rawStreamZDC->IsOverflow()) && (detector!=-1)
344              && ((rawStreamZDC->GetADCGain()) == 0) && // Selecting HIGH RES ch.s
345              (rawStreamZDC->GetADCModule()>=kFirstADCGeo) && (rawStreamZDC->GetADCModule()<=kLastADCGeo)){
346
347             //printf("  IsADCWord %d, IsUnderflow %d, IsOverflow %d\n",
348             //  rawStreamZDC->IsADCDataWord(),rawStreamZDC->IsUnderflow(),rawStreamZDC->IsOverflow());
349           
350             if(quad!=5){ // Physics signals
351               if(detector==1)      index = quad;   // *** ZNC
352               else if(detector==2) index = quad+5; // *** ZPC
353               else if(detector==3) index = quad+9; // *** ZEM
354               else if(detector==4) index = quad+12;// *** ZNA
355               else if(detector==5) index = quad+17;// *** ZPA
356             }
357             else continue;
358             //
359             if(index==-1) printf("ERROR in ZDCCALIBMBda.cxx -> det %d quad %d index %d ADC %d\n", 
360               detector, quad, index, rawStreamZDC->GetADCValue());
361             
362             // Mean pedestal subtraction 
363             Float_t Pedestal = MeanPedhg[index];
364             // Pedestal subtraction from correlation with out-of-time signals
365             //Float_t Pedestal = CorrCoeff0[index]+CorrCoeff1[index]*MeanPedOOT[index];
366             //
367             if(index!=-1 && quad!=5){ 
368               //
369               if(detector==3){
370                rawZEM  = (Float_t) rawStreamZDC->GetADCValue();
371                corrZEM = (rawStreamZDC->GetADCValue()) - Pedestal;
372                corrZEMTot += corrZEM;
373               }
374               else{
375                rawZDC  = (Float_t) rawStreamZDC->GetADCValue();
376                corrZDC = rawZDC - Pedestal;
377                if(detector==1 || detector==2) corrZDCC += corrZDC;
378                else corrZDCA += corrZDC;
379                corrZDCTot += corrZDC;
380               }
381             }
382         }//IsADCDataWord()
383         
384        } // Next()
385         
386        hZDCvsZEM ->Fill(corrZDCTot, corrZEMTot);
387        hZDCCvsZEM->Fill(corrZDCC, corrZEMTot);
388        hZDCAvsZEM->Fill(corrZDCA, corrZEMTot);
389
390        nevents_physics++;
391     }//(if PHYSICS_EVENT)
392     
393     nevents_total++;
394
395    }
396
397    /* free resources */
398    free(event);
399   }
400     
401   //
402   TFile *histofile = new TFile(MBCALIBDATA_FILE,"RECREATE");
403   histofile->cd();
404   hZDCvsZEM ->Write();
405   hZDCCvsZEM->Write();
406   hZDCAvsZEM->Write();
407   histofile->Close();
408   
409   if(hZDCvsZEM)  delete hZDCvsZEM ;
410   if(hZDCCvsZEM) delete hZDCCvsZEM;
411   if(hZDCAvsZEM) delete hZDCAvsZEM;
412   
413   /* write report */
414   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
415
416   /* close result file */
417   fclose(fp);
418   
419   /* report progress */
420   daqDA_progressReport(90);
421
422   /* store the result file on FES */
423   status = daqDA_FES_storeFile(MAPDATA_FILE, "MAPPING");
424   if(status){
425     printf("Failed to export file : %d\n",status);
426     return -1;
427   }
428   //
429   status = daqDA_FES_storeFile(MBCALIBDATA_FILE, "MBCALIB");
430   if(status){
431     printf("Failed to export file : %d\n",status);
432     return -1;
433   }
434  
435   /* report progress */
436   daqDA_progressReport(100);
437
438   return status;
439 }