Initial version of channel equalization DA for Pb-Pb data.
[u/mrichter/AliRoot.git] / VZERO / VZEROPbPbda.cxx
CommitLineData
2ffa87b2 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
43Int_t GetOfflineChannel(Int_t channel);
44
45/* Main routine --- Arguments: monitoring data source */
46
47int 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}