Initial version of channel equalization DA for Pb-Pb data.
[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
64   Int_t    kStartClock = 9;  // First clock in the search for max adc
65   Int_t    kEndClock = 11;   // Last clock in the search for max adc
66   Int_t    kNPreClocks = 6;  // Number of clock before max used in the charge sum
67   Int_t    kNPostClocks = 1; // Number of clock after max used in the charge sum
68
69   UShort_t    kTriggerAcc = 64;    // Trigger mask for accepted events (64 = CTA1 & CTC1)
70   UShort_t    kTriggerRej = 256;   // Trigger mask for rejected events (256 = CTA2 & CTC2)
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 %u %u %d %f",
83                      &kClockStart,&kClockStop,&kNPreClock,&kNPostClock,&kTriggerAcc,&kTriggerRej,&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("First LHC Clock = %d; Last LHC Clock = %d; N Pre Clock = %d ; N Post Clock = %d; Trigger mask for accepted events = %u; Trigger mask for rejected events = %u; Number of histogram bins = %d; Histogram range = %.3f\n",
91          kClockStart, kClockStop, kNPreClock, kNPostClock, kTriggerAcc, kTriggerRej, 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 //___________________________________________________ 
97
98   /* open result file to be exported to FES */
99   FILE *fp=NULL;
100   fp=fopen("./V0_EqualizationFactors.dat","w");
101   if (fp==NULL) {
102       printf("Failed to open local result file\n");
103       return -1;}
104    
105   /* define data source : this is argument 1 */  
106   status=monitorSetDataSource( argv[1] );
107   if (status!=0) {
108      printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
109      return -1;
110   }
111
112   /* declare monitoring program */
113   status=monitorDeclareMp( __FILE__ );
114   if (status!=0) {
115     printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
116     return -1;
117   }
118
119   /* define wait event timeout - 1s max */
120   monitorSetNowait();
121   monitorSetNoWaitNetworkTimeout(1000);
122   
123   /* init counters on events */
124   int neventsPhysics=0;
125   int neventsTotal=0;
126
127
128   /* loop on events (infinite) */
129   for(;;) {
130       struct eventHeaderStruct *event;
131       eventTypeType eventT;
132
133       /* check shutdown condition */
134       if (daqDA_checkShutdown()) {break;}   
135
136       /* get next event (blocking call until timeout) */
137       status=monitorGetEventDynamic((void **)&event);
138       if (status==MON_ERR_EOF) {
139           printf ("End of File detected\n");
140       break; /* end of monitoring file has been reached */
141       }
142     
143       if (status!=0) {
144            printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
145       break;
146       }
147
148       /* retry if got no event */
149       if (event==NULL) continue;
150                
151       /* decode event */
152       eventT=event->eventType;
153         
154       switch (event->eventType){
155       
156       case START_OF_RUN:
157            break;
158       
159       case END_OF_RUN:
160            printf("End Of Run detected\n");
161            break;
162       
163       case PHYSICS_EVENT:
164                          
165            AliRawReader *rawReader = new AliRawReaderDate((void*)event);
166   
167            AliVZERORawStream* rawStream  = new AliVZERORawStream(rawReader); 
168            if (rawStream->Next()) {     
169              UShort_t triggers = rawStream->GetTriggerInputs();
170              if (((triggers & kTriggerAcc) == kTriggerAcc) &&  // Check if the requested trigger(s) is fired
171                  ((triggers & kTriggerRej) == 0)) { // Check if the requested trigger(s) is NOT fired
172                neventsPhysics++;
173
174                Float_t adc[64];
175                Float_t totCharge = 0.;
176                for(Int_t i = 0; i < 64; ++i) {
177                  adc[i] = 0.;
178                  Float_t maxadc=0.;
179                  Int_t imax=-1;
180                  for(Int_t j = kStartClock; j <= kEndClock; ++j) {
181                    Float_t charge = (Float_t)(rawStream->GetPedestal(i,j));
182                    if(charge > maxadc){
183                      maxadc = charge;
184                      imax   = j;
185                    }
186                  }
187
188                  if (imax != -1) {
189                    Int_t start = imax - kNPreClocks;
190                    if (start < 0) start = 0;
191                    Int_t end = imax + kNPostClocks;
192                    if (end > 20) end = 20;
193                    for(Int_t iClock = start; iClock <= end; iClock++) {
194                      adc[i] += (Float_t)(rawStream->GetPedestal(i,iClock));
195                    }
196                  }
197                  totCharge += adc[i];
198                }
199
200                if (fFirst) {
201                  fFirst = kFALSE;
202                  fPrevTotCharge = totCharge;
203                  for(int i = 0; i < 64; ++i) fPrevadc[i] = adc[i];
204                }
205                else {
206                  fFirst = kTRUE;
207                  Float_t deltaTotCharge = totCharge - fPrevTotCharge;
208                  Float_t weight = deltaTotCharge*deltaTotCharge;
209                  if (weight > 1) {
210                    for(int i = 0; i < 64; ++i) {
211                      fMedian[i]->Fill((adc[i]-fPrevadc[i])/deltaTotCharge,weight);
212                    }
213                  }
214                }
215
216              }
217            }   // End : if rawstream
218            delete rawStream;
219            rawStream = 0x0;      
220            delete rawReader;
221            rawReader = 0x0;                                                                      
222       } // end of switch on event type 
223         
224       neventsTotal++;
225       /* free resources */
226       free(event);
227         
228       /* exit when last event received, no need to wait for TERM signal */
229       if (eventT==END_OF_RUN) {
230         printf("End Of Run event detected\n");
231         break;
232       }
233
234   }  // loop over events
235   
236   printf("%d physics events processed\n",neventsPhysics);
237     
238 //___________________________________________________________________________
239 //  Computes regression parameters
240 // charge_i = p0 + charge_tot * p1
241
242   Double_t beta[64];
243   Double_t q = 0.5;
244   for(int i = 0; i < 64; ++i) fMedian[i]->GetQuantiles(1,&beta[i],&q);
245
246   for(Int_t i=0; i<64; i++) {
247     fprintf(fp," %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);                                
248     printf(" %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);                                    
249   }
250  
251 //________________________________________________________________________
252    
253   /* close local result file and FXS result file*/
254   fclose(fp);
255
256   if(neventsPhysics>2000){
257     /* export result file to FES */
258     status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
259     if (status)    {
260       printf("Failed to export file : %d\n",status);
261       return -1; }
262
263     /* store result file into Online DB */
264     status=daqDA_DB_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
265     if (status)    {
266       printf("Failed to store file into Online DB: %d\n",status);
267       return -1; }
268   }
269
270   return status;
271 }
272
273  Int_t GetOfflineChannel(Int_t channel) {
274
275 // Channel mapping Online - Offline:
276  
277  Int_t fOfflineChannel[64] = {39, 38, 37, 36, 35, 34, 33, 32, 
278                               47, 46, 45, 44, 43, 42, 41, 40, 
279                               55, 54, 53, 52, 51, 50, 49, 48, 
280                               63, 62, 61, 60, 59, 58, 57, 56,
281                                7,  6,  5,  4,  3,  2,  1,  0, 
282                               15, 14, 13, 12, 11, 10,  9,  8,
283                               23, 22, 21, 20, 19, 18, 17, 16, 
284                               31, 30, 29, 28, 27, 26, 25, 24};
285  return fOfflineChannel[channel];                     
286 }