correct mask for V0 charge decoding in STU payload
[u/mrichter/AliRoot.git] / VZERO / VZEROPbPbda.cxx
1 /*********************************************************************************
2 - Contact:    Brigitte Cheynis     b.cheynis@ipnl.in2p3.fr
3 - Link:
4 - Raw data test file :          
5 - Reference run number : 137366       
6 - Run Type:   PHYSICS
7 - DA Type:    MON
8 - Number of events needed: >=2000
9 - Input Files:  argument list
10 - Output Files: FXS file     V0_EqualizationFactors.dat (Channel equalization factors)
11 - Trigger types used: PHYSICS_EVENT
12 **********************************************************************************/
13
14 /**********************************************************************************
15 *                                                                                 *
16 * VZERO Detector Algorithm for extracting channel equalization factors for Pb-Pb  *
17 *                                                                                 *
18 *                                                                                 *
19 ***********************************************************************************/
20
21 // DATE
22 #include "event.h"
23 #include "monitor.h"
24 #include "daqDA.h"
25
26 //AliRoot
27 #include <AliVZERORawStream.h>
28 #include <AliRawReaderDate.h>
29 #include <AliRawReader.h>
30 #include <AliDAQ.h>
31
32 // standard
33 #include <stdio.h>
34 #include <stdlib.h>
35
36 //ROOT
37 #include "TROOT.h"
38 #include "TPluginManager.h"
39 #include <TFile.h>
40 #include <TH1F.h>
41 #include <TMath.h>
42
43 Int_t GetOfflineChannel(Int_t channel);
44
45 /* Main routine --- Arguments: monitoring data source */
46       
47 int main(int argc, char **argv) {
48
49 /* magic line from Cvetan */
50   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
51                     "*",
52                     "TStreamerInfo",
53                     "RIO",
54                     "TStreamerInfo()");
55   int status;
56   if (argc!=2) {
57      printf("Wrong number of arguments\n");
58      return -1;
59   }
60   
61 //___________________________________________________
62 // Get parameters from V00DAEqualFactors.config file
63   Int_t    neventsMin = 2000;
64
65   Int_t    kStartClock = 9;  // First clock in the search for max adc
66   Int_t    kEndClock = 11;   // Last clock in the search for max adc
67   Int_t    kNPreClocks = 6;  // Number of clock before max used in the charge sum
68   Int_t    kNPostClocks = 1; // Number of clock after max used in the charge sum
69
70   UChar_t  kNFlagsCut   = 63;// Minimum number of TDC flags requested, replaces the trigger selection
71
72   Int_t    kNBins = 10000;
73   Float_t  kRange = 0.1;
74
75   status = daqDA_DB_getFile("V00DAEqualFactors.config","./V00DAEqualFactors.config");
76   if (status) {
77     printf("Failed to get Config file (V00DAEqualFactors.config) from DAQ DB, status=%d\n", status);
78     printf("Take default values of parameters for pedestal calculation \n");
79   } else {
80     /* open the config file and retrieve cuts */
81     FILE *fpConfig = fopen("V00DAEqualFactors.config","r");
82     int res = fscanf(fpConfig,"%d %d %d %d %d %hhu %d %f",
83                      &neventsMin,&kStartClock,&kEndClock,&kNPreClocks,&kNPostClocks,&kNFlagsCut,&kNBins,&kRange);
84     if(res!=8) {
85       printf("Failed to get values from Config file (V00DAEqualFactors.config): wrong file format - 7 integers and 1 float are expected - \n");
86     }
87     fclose(fpConfig);
88   }
89   
90   printf("Minimum number of events to process = %d; First LHC Clock = %d; Last LHC Clock = %d; N Pre Clock = %d; N Post Clock = %d; Minimum number of TDC flags = %hhu; Number of histogram bins = %d; Histogram range = %.3f\n",
91          neventsMin,kStartClock, kEndClock, kNPreClocks, kNPostClocks, kNFlagsCut, kNBins, kRange);
92
93   TH1D *fMedian[64];
94   for(Int_t j = 0; j < 64; ++j) fMedian[j] = new TH1D(Form("fMedian_%d",j),"Slopes weighted median, channel par channel",kNBins,0,kRange);
95
96   Bool_t fFirst = kTRUE;
97   Float_t fPrevTotCharge = 0;
98   Float_t fPrevadc[64];
99   for(Int_t j = 0; j < 64; ++j) fPrevadc[j] = 0;
100
101 //___________________________________________________ 
102    
103   /* define data source : this is argument 1 */  
104   status=monitorSetDataSource( argv[1] );
105   if (status!=0) {
106      printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
107      return -1;
108   }
109
110   /* declare monitoring program */
111   status=monitorDeclareMp( __FILE__ );
112   if (status!=0) {
113     printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
114     return -1;
115   }
116
117   /* define wait event timeout - 1s max */
118   monitorSetNowait();
119   monitorSetNoWaitNetworkTimeout(1000);
120   
121   /* init counters on events */
122   int neventsPhysicsAll=0;
123   int neventsPhysics=0;
124   int neventsTotal=0;
125
126
127   /* loop on events (infinite) */
128   for(;;) {
129       struct eventHeaderStruct *event;
130       eventTypeType eventT;
131
132       /* check shutdown condition */
133       if (daqDA_checkShutdown()) {break;}   
134
135       /* get next event (blocking call until timeout) */
136       status=monitorGetEventDynamic((void **)&event);
137       if (status==MON_ERR_EOF) {
138           printf ("End of File detected\n");
139       break; /* end of monitoring file has been reached */
140       }
141     
142       if (status!=0) {
143            printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
144       break;
145       }
146
147       /* retry if got no event */
148       if (event==NULL) continue;
149                
150       /* decode event */
151       eventT=event->eventType;
152         
153       switch (event->eventType){
154       
155       case START_OF_RUN:
156            break;
157       
158       case END_OF_RUN:
159            printf("End Of Run detected\n");
160            break;
161       
162       case PHYSICS_EVENT:
163                          
164            AliRawReader *rawReader = new AliRawReaderDate((void*)event);
165   
166            AliVZERORawStream* rawStream  = new AliVZERORawStream(rawReader); 
167            if (rawStream->Next()) {     
168              neventsPhysicsAll++;
169              UShort_t triggers = rawStream->GetTriggerInputs();
170              UChar_t nFlags = 0;
171              for(Int_t i = 0; i < 64; ++i) {
172                if (rawStream->GetBBFlag(i,10)) nFlags++;
173              }
174              if (nFlags >= kNFlagsCut) { // Check if the number of minimum number of TDC flags is reached (trigger selection)
175
176                neventsPhysics++;
177
178                Float_t adc[64];
179                Float_t totCharge = 0.;
180                for(Int_t i = 0; i < 64; ++i) {
181                  adc[i] = 0.;
182                  Float_t maxadc=0.;
183                  Int_t imax=-1;
184                  for(Int_t j = kStartClock; j <= kEndClock; ++j) {
185                    Float_t charge = (Float_t)(rawStream->GetPedestal(i,j));
186                    if(charge > maxadc){
187                      maxadc = charge;
188                      imax   = j;
189                    }
190                  }
191
192                  if (imax != -1) {
193                    Int_t start = imax - kNPreClocks;
194                    if (start < 0) start = 0;
195                    Int_t end = imax + kNPostClocks;
196                    if (end > 20) end = 20;
197                    for(Int_t iClock = start; iClock <= end; iClock++) {
198                      adc[i] += (Float_t)(rawStream->GetPedestal(i,iClock));
199                    }
200                  }
201                  totCharge += adc[i];
202                }
203
204                if (fFirst) {
205                  fFirst = kFALSE;
206                  fPrevTotCharge = totCharge;
207                  for(int i = 0; i < 64; ++i) fPrevadc[i] = adc[i];
208                }
209                else {
210                  fFirst = kTRUE;
211                  Float_t deltaTotCharge = totCharge - fPrevTotCharge;
212                  Float_t weight = deltaTotCharge*deltaTotCharge;
213                  if (weight > 1) {
214                    for(int i = 0; i < 64; ++i) {
215                      fMedian[i]->Fill((adc[i]-fPrevadc[i])/deltaTotCharge,weight);
216                    }
217                  }
218                }
219
220              }
221            }   // End : if rawstream
222            delete rawStream;
223            rawStream = 0x0;      
224            delete rawReader;
225            rawReader = 0x0;                                                                      
226       } // end of switch on event type 
227         
228       neventsTotal++;
229       /* free resources */
230       free(event);
231         
232       /* exit when last event received, no need to wait for TERM signal */
233       if (eventT==END_OF_RUN) {
234         printf("End Of Run event detected\n");
235         break;
236       }
237
238   }  // loop over events
239   
240   printf("%d physics events processed (out of %d physics and %d total events)\n",neventsPhysics,neventsPhysicsAll,neventsTotal);
241     
242 //___________________________________________________________________________
243 //  Computes regression parameters
244 // charge_i = p0 + charge_tot * p1
245
246   if(neventsPhysics>neventsMin){
247     /* open result file to be exported to FES */
248     FILE *fp=NULL;
249     fp=fopen("./V0_EqualizationFactors.dat","w");
250     if (fp==NULL) {
251       printf("Failed to open local result file\n");
252       return -1;}
253
254     Double_t beta[64];
255     Double_t q = 0.5;
256     for(int i = 0; i < 64; ++i) fMedian[i]->GetQuantiles(1,&beta[i],&q);
257
258     for(Int_t i=0; i<64; i++) {
259       fprintf(fp," %d %.5f\n",GetOfflineChannel(i), beta[i]*64.);                                      
260       printf(" %d %.5f\n",GetOfflineChannel(i), beta[i]*64.);                                  
261     }
262
263     /* close local result file and FXS result file*/
264     fclose(fp);
265    
266     /* export result file to FES */
267     status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
268     if (status)    {
269       printf("Failed to export file : %d\n",status);
270       return -1; }
271
272     /* store result file into Online DB */
273     status=daqDA_DB_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
274     if (status)    {
275       printf("Failed to store file into Online DB: %d\n",status);
276       return -1; }
277   }
278   else {
279     // Take the last run's file from the DB
280     status=daqDA_DB_getFile("V00DAEqualFactors","./V0_EqualizationFactors.dat");
281     if (status)    {
282       printf("Failed to get file from Online DB: %d\n",status);
283       return -1; }
284
285     /* export result file to FES */
286     status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
287     if (status)    {
288       printf("Failed to export file : %d\n",status);
289       return -1; }
290   }
291
292   return status;
293 }
294
295  Int_t GetOfflineChannel(Int_t channel) {
296
297 // Channel mapping Online - Offline:
298  
299  Int_t fOfflineChannel[64] = {39, 38, 37, 36, 35, 34, 33, 32, 
300                               47, 46, 45, 44, 43, 42, 41, 40, 
301                               55, 54, 53, 52, 51, 50, 49, 48, 
302                               63, 62, 61, 60, 59, 58, 57, 56,
303                                7,  6,  5,  4,  3,  2,  1,  0, 
304                               15, 14, 13, 12, 11, 10,  9,  8,
305                               23, 22, 21, 20, 19, 18, 17, 16, 
306                               31, 30, 29, 28, 27, 26, 25, 24};
307  return fOfflineChannel[channel];                     
308 }