]>
Commit | Line | Data |
---|---|---|
1 | /* | |
2 | CPV BCM DA for processing physics runs and producing bad channel map for further use at reconstruction time. | |
3 | ||
4 | Contact: Sergey Evdokimov <sevdokim@cern.ch> | |
5 | Link: https://twiki.cern.ch/twiki/bin/view/ALICE/CPVda | |
6 | Reference run: 214340 (/afs/cern.ch/user/s/sevdokim/public/CPV_run214340_standalone.raw) | |
7 | Run Type: PHYSICS | |
8 | DA Type: MON | |
9 | Number of events needed: ~100000 events | |
10 | Input files: thr?_??.dat CpvPeds.root PHOSCPVBCMda.cfg | |
11 | Output files: CpvBadMap.root | |
12 | Trigger types used: PHYSICS_EVENT | |
13 | */ | |
14 | ||
15 | //daqDA | |
16 | #include "event.h" | |
17 | #include "monitor.h" | |
18 | #include "daqDA.h" | |
19 | //AMORE monitoring framework | |
20 | #include <AmoreDA.h> | |
21 | ||
22 | //system | |
23 | #include <Riostream.h> | |
24 | #include <stdlib.h> | |
25 | #include <fstream> | |
26 | #include <string> | |
27 | ||
28 | //AliRoot | |
29 | #include "AliPHOSCpvRawDigiProducer.h" | |
30 | #include "AliPHOSCpvParam.h" | |
31 | #include "AliRawReaderDate.h" | |
32 | #include "AliBitPacking.h" | |
33 | #include "AliPHOSDigit.h" | |
34 | #include "AliPHOSGeometry.h" | |
35 | ||
36 | //ROOT | |
37 | #include "TROOT.h" | |
38 | #include "TPluginManager.h" | |
39 | #include "TSAXParser.h" | |
40 | #include "TTree.h" | |
41 | #include "TMath.h" | |
42 | #include "TString.h" | |
43 | #include "TFile.h" | |
44 | #include "TSystem.h" | |
45 | #include "TKey.h" | |
46 | #include "TH2S.h" | |
47 | #include "TH2F.h" | |
48 | #include "TObject.h" | |
49 | #include "TBenchmark.h" | |
50 | #include "TMath.h" | |
51 | #include "TRandom.h" | |
52 | ||
53 | Double_t FillNoisyMap(TH2* DigMap, TH2* BadMap); //returns mean occupancy when all bad channels are excluded | |
54 | void FillDeadMap(TH2* PedMap, TH2* DigMap, TH2* BadMap); | |
55 | ||
56 | int main( int argc, char **argv ) | |
57 | { | |
58 | gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo", | |
59 | "*", | |
60 | "TStreamerInfo", | |
61 | "RIO", | |
62 | "TStreamerInfo()"); | |
63 | Int_t status,statusPeds=0,print,statusBadMap=0; | |
64 | Int_t sigcut=3; | |
65 | Int_t minAmpl = 10;//minimal amplitude for consideration, to be read from DAQ DB | |
66 | Int_t minOccupancy = 10;//min occupancy for publishing in OCDB, to be read from DAQ DB | |
67 | Bool_t turbo = kTRUE; | |
68 | ||
69 | if (argc!=2) { | |
70 | printf("Wrong number of arguments\n"); | |
71 | return -1; | |
72 | } | |
73 | ||
74 | // log start of process | |
75 | printf("Cpv bad channel map DA program started\n"); | |
76 | ||
77 | /* report progress */ | |
78 | daqDA_progressReport(0); | |
79 | ||
80 | /* retrieve configuration file from DAQ DB */ | |
81 | status=daqDA_DB_getFile("PHOSCPVGAINda.cfg", "PHOSCPVGAINda.cfg"); | |
82 | if(!status) { | |
83 | char buf[500]; | |
84 | FILE * fConf = fopen("PHOSCPVGAINda.cfg","r"); | |
85 | while(fgets(buf, 500, fConf)){ | |
86 | if(buf[0]=='#') continue;//comment indicator | |
87 | if(strstr(buf,"minOccupancy")) sscanf(buf,"%*s %d",&minOccupancy); | |
88 | if(strstr(buf,"minAmpl")) sscanf(buf,"%*s %d",&minAmpl); | |
89 | } | |
90 | } | |
91 | ||
92 | ||
93 | /* retrieve pedestal tables from DAQ DB */ | |
94 | for(int iDDL = 0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL+=2){ | |
95 | if(iDDL!=4) continue; // only one module with DDL=4 by now | |
96 | for (int iCC = 0; iCC<AliPHOSCpvParam::kNRows; iCC++){ | |
97 | status=daqDA_DB_getFile(Form("thr%d_%02d.dat", iDDL, iCC),Form("thr%d_%02d.dat", iDDL, iCC)); | |
98 | if(status!=0) { | |
99 | printf("cannot retrieve file %s from DAQ DB. Exit.\n", Form("thr%d_%02d.dat", iDDL, iCC)); | |
100 | //return 11;//error code 11 (cannot retrive thr.dat from DAQ DB) | |
101 | } | |
102 | } | |
103 | } | |
104 | ||
105 | /* retrieve pedestals in root format to find dead channels */ | |
106 | statusPeds=daqDA_DB_getFile("CpvPeds.root", "CpvPeds.root"); | |
107 | if(statusPeds) { | |
108 | printf("cannot retrieve CpvPeds.root from DAQ DB! No dead channels will be found."); | |
109 | //return 12; //error code 12 (cannot retrive CpvPeds.root from DAQ DB) | |
110 | } | |
111 | ||
112 | /* retrieve bad map from DAQ DB to see if we have some statistics saved form previous runs */ | |
113 | statusBadMap=daqDA_DB_getFile("CpvBadMap.root", "CpvBadMap.root"); | |
114 | ||
115 | /* connecting to raw data */ | |
116 | status=monitorSetDataSource( argv[1] ); | |
117 | if (status!=0) { | |
118 | printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status)); | |
119 | return -1; | |
120 | } | |
121 | ||
122 | /* declare monitoring program */ | |
123 | status=monitorDeclareMp( __FILE__ ); | |
124 | if (status!=0) { | |
125 | printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status)); | |
126 | return -1; | |
127 | } | |
128 | ||
129 | /* define wait event timeout - 1s max */ | |
130 | monitorSetNowait(); | |
131 | monitorSetNoWaitNetworkTimeout(1000); | |
132 | ||
133 | /* report progress */ | |
134 | daqDA_progressReport(5); | |
135 | ||
136 | ||
137 | // init event counter | |
138 | Int_t iPhysEvnt=0; | |
139 | Int_t iTotEvnt =0; | |
140 | ||
141 | // Reader | |
142 | AliRawReader * reader; | |
143 | ||
144 | //digiProducer | |
145 | AliPHOSCpvRawDigiProducer* digiProducer = new AliPHOSCpvRawDigiProducer(); | |
146 | digiProducer->SetTurbo(turbo); | |
147 | digiProducer->LoadPedFiles(); | |
148 | digiProducer->SetCpvMinAmp(minAmpl); | |
149 | ||
150 | //digits | |
151 | TClonesArray *digits = new TClonesArray("AliPHOSDigit",1); | |
152 | digits->SetName("DIGITS"); | |
153 | ||
154 | const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ; | |
155 | ||
156 | //maps of digits and bad channels | |
157 | TH2F* hMapOfDig[2*AliPHOSCpvParam::kNDDL]; | |
158 | TH2I *hBadChMap[2*AliPHOSCpvParam::kNDDL]; | |
159 | ||
160 | for (int i = 0;i<2*AliPHOSCpvParam::kNDDL;i++){ | |
161 | hMapOfDig[i]= 0x0; | |
162 | hBadChMap[i]= 0x0; | |
163 | } | |
164 | ||
165 | //any previously gained statistics? | |
166 | TFile *fPreviousStatistics = 0x0; | |
167 | if(!statusBadMap) fPreviousStatistics = TFile::Open("CpvBadMap.root"); | |
168 | ||
169 | for(Int_t iDDL=0;iDDL<2*AliPHOSCpvParam::kNDDL;iDDL++){ | |
170 | if(!statusBadMap){//we have some statistics from previous runs | |
171 | if(fPreviousStatistics->Get(Form("hMapOfDig%d",iDDL))) | |
172 | hMapOfDig[iDDL] = new TH2F(*(TH2F*)(fPreviousStatistics->Get(Form("hMapOfDig%d",iDDL)))); | |
173 | } | |
174 | if(!hMapOfDig[iDDL]) | |
175 | hMapOfDig[iDDL] = new TH2F(Form("hMapOfDig%d",iDDL),Form("Map of digits with substructed pedestals, DDL = %d",iDDL), | |
176 | AliPHOSCpvParam::kPadPcX,0,AliPHOSCpvParam::kPadPcX, | |
177 | AliPHOSCpvParam::kPadPcY,0,AliPHOSCpvParam::kPadPcY); | |
178 | hBadChMap[iDDL] = new TH2I(Form("hBadMap%d",iDDL),Form("Bad Channels Map, DDL= %d",iDDL), | |
179 | AliPHOSCpvParam::kPadPcX,0,AliPHOSCpvParam::kPadPcX, | |
180 | AliPHOSCpvParam::kPadPcY,0,AliPHOSCpvParam::kPadPcY); | |
181 | ||
182 | } | |
183 | /* report progress */ | |
184 | daqDA_progressReport(10); | |
185 | ||
186 | /* main loop (infinite) */ | |
187 | for(;;) { // infinite loop | |
188 | struct eventHeaderStruct *event; | |
189 | eventTypeType eventT; | |
190 | ||
191 | /* check shutdown condition */ | |
192 | if (daqDA_checkShutdown()) {break;} | |
193 | ||
194 | // get next event | |
195 | status=monitorGetEventDynamic((void **)&event); | |
196 | if (status==MON_ERR_EOF) { // end of monitoring file has been reached | |
197 | printf("End of monitoring file has been reached! \n"); | |
198 | break; | |
199 | } | |
200 | ||
201 | if (status!=0) { | |
202 | printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status)); | |
203 | break; | |
204 | } | |
205 | ||
206 | // retry if got no event | |
207 | if (event==NULL) continue; | |
208 | ||
209 | // use event - here, just write event id to result file | |
210 | eventT=event->eventType; | |
211 | if (eventT==PHYSICS_EVENT) { //we use PHYSICS_EVENT for pedestal not CALIBRATION_EVENT??? | |
212 | iTotEvnt++; | |
213 | reader = new AliRawReaderDate((void*)event); | |
214 | digiProducer->LoadNewEvent(reader); | |
215 | digiProducer->MakeDigits(digits); | |
216 | if(digits->GetEntriesFast()>0) iPhysEvnt++; | |
217 | for(Int_t i=0;i<digits->GetEntriesFast();i++) | |
218 | { | |
219 | Int_t relId[4]; | |
220 | AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ; | |
221 | geom->AbsToRelNumbering(digit->GetId(),relId); | |
222 | //cout<<"relId : "<<relId[0]<<" "<<relId[1]<<" "<<relId[2]<<" "<<relId[3]<<endl; | |
223 | if( | |
224 | (relId[0]>=0&&relId[0]<=AliPHOSCpvParam::kNDDL) | |
225 | &&(relId[1]==-1)//cpv | |
226 | &&(relId[2]-1>=0&&relId[2]-1<AliPHOSCpvParam::kPadPcX) | |
227 | &&(relId[3]-1>=0&&relId[3]-1<AliPHOSCpvParam::kPadPcY) | |
228 | ) | |
229 | hMapOfDig[AliPHOSCpvParam::Mod2DDL(relId[0])]->Fill(relId[2]-1,relId[3]-1); | |
230 | } | |
231 | digits->Clear("C"); | |
232 | delete reader; | |
233 | } // if PHYSICS_EVENT | |
234 | ||
235 | free(event); | |
236 | ||
237 | /* exit when last event received, no need to wait for TERM signal */ | |
238 | if (eventT==END_OF_RUN) { | |
239 | printf("EOR event detected\n"); | |
240 | break; | |
241 | } | |
242 | } | |
243 | ||
244 | Printf(" Received %d events, %d good events",iTotEvnt,iPhysEvnt); | |
245 | /* report progress */ | |
246 | daqDA_progressReport(90); | |
247 | ||
248 | TH2* hPedMap[2*AliPHOSCpvParam::kNDDL]; | |
249 | //prepare ped maps for dead channels search | |
250 | TFile* fPeds = TFile::Open("CpvPeds.root"); | |
251 | for(int iDDL = 0; iDDL< 2*AliPHOSCpvParam::kNDDL; iDDL+=2){ | |
252 | if(fPeds->Get(Form("fPedMeanMap%d",iDDL))) | |
253 | hPedMap[iDDL] = new TH2F(*(TH2F*)(fPeds->Get(Form("fPedMeanMap%d",iDDL)))); | |
254 | else | |
255 | hPedMap[iDDL] = 0x0; | |
256 | } | |
257 | fPeds->Close(); | |
258 | ||
259 | ||
260 | ||
261 | //find noisy channels (i.e. channelOccupancy > 10*meanOccupancy) | |
262 | TFile *fSave = TFile::Open("CpvBadMap.root","RECREATE"); | |
263 | ||
264 | Double_t meanOccupancy = 0; | |
265 | amore::da::AmoreDA* myAmore = new amore::da::AmoreDA(amore::da::AmoreDA::kSender); | |
266 | for(int iDDL = 0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL+=2){ | |
267 | if(hMapOfDig[iDDL]->GetEntries()>0) { | |
268 | Double_t Occupancy = FillNoisyMap(hMapOfDig[iDDL],hBadChMap[iDDL]); | |
269 | if(meanOccupancy>0&&Occupancy<meanOccupancy) meanOccupancy = Occupancy; | |
270 | if(meanOccupancy==0) meanOccupancy = Occupancy; | |
271 | if(Occupancy>minOccupancy) FillDeadMap(hPedMap[iDDL],hMapOfDig[iDDL],hBadChMap[iDDL]); | |
272 | fSave->WriteObject(hMapOfDig[iDDL],Form("hMapOfDig%d",iDDL)); | |
273 | //send digit maps to amore | |
274 | myAmore->Send(Form("hMapOfDig%d",iDDL),hMapOfDig[iDDL]); | |
275 | fSave->WriteObject(hBadChMap[iDDL],Form("hBadChMap%d",iDDL)); | |
276 | } | |
277 | } | |
278 | fSave->Close(); | |
279 | cout<< "meanOccupancy = "<<meanOccupancy<<"; minOccupancy = "<<minOccupancy<<endl; | |
280 | ||
281 | if(meanOccupancy>minOccupancy){//send file to FES if only statistics is enough | |
282 | status = daqDA_FES_storeFile("CpvBadMap.root","CpvBadMap.root"); | |
283 | if(status) printf("Failed to store CpvBadMap.root in DAQ FXS!\n"); | |
284 | //store dummy file in DAQ DB | |
285 | TFile *fDummy = TFile::Open("dummy.root","RECREATE"); | |
286 | fDummy->Close(); | |
287 | status = daqDA_DB_storeFile("dummy.root","CpvBadMap.root"); | |
288 | if(status) printf("Failed to store dummy.root as CpvBadMap.root in DAQ DB!\n"); | |
289 | //send bad map to amore as well | |
290 | for(int iDDL = 0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL+=2) | |
291 | if(hMapOfDig[iDDL]->GetEntries()>0) | |
292 | myAmore->Send(Form("hBadChMap%d",iDDL),hBadChMap[iDDL]); | |
293 | } | |
294 | else{//store file with current statistics in DAQ DB for further use. | |
295 | status = daqDA_DB_storeFile("CpvBadMap.root","CpvBadMap.root"); | |
296 | if(status) printf("Failed to store CpvBadMap.root in DAQ DB!\n"); | |
297 | ||
298 | } | |
299 | ||
300 | /* report progress */ | |
301 | daqDA_progressReport(100); | |
302 | ||
303 | ||
304 | return status; | |
305 | } | |
306 | //============================================================================== | |
307 | Double_t FillNoisyMap(TH2* hDigMap, TH2* hBadChMap){ | |
308 | if(!hDigMap)return 0; | |
309 | if(!hBadChMap)return 0; | |
310 | Double_t meanOccupancy = hDigMap->GetEntries()/7680.; | |
311 | Double_t nDigits=0,nChannels=0; | |
312 | int x,y,iterationNumber=1; | |
313 | int nBadChannelsTotal=0,nBadChannelsCurrentIteration=0; | |
314 | //1st iteration | |
315 | //cout<<"Iteration Number = "<<iterationNumber<<endl; | |
316 | for(int ix = 1;ix<=128;ix++) | |
317 | for(int iy = 1;iy<=60;iy++) | |
318 | if(hDigMap->GetBinContent(ix,iy)>meanOccupancy*10.+1.){ | |
319 | nBadChannelsCurrentIteration++; | |
320 | hBadChMap->Fill(ix-1,iy-1); | |
321 | printf("Noisy channel found! DDL=4 x=%d y=%d\n",ix-1,iy-1); | |
322 | } | |
323 | //next iterations | |
324 | while(nBadChannelsCurrentIteration!=0){ | |
325 | iterationNumber++; | |
326 | //cout<<"Iteration Number = "<<iterationNumber<<endl; | |
327 | nDigits=0;nChannels=0; | |
328 | nBadChannelsTotal+=nBadChannelsCurrentIteration; | |
329 | nBadChannelsCurrentIteration=0; | |
330 | //1st -- calculate new mean occupancy excluding already badly marked channels | |
331 | for(int ix = 1;ix<=128;ix++) | |
332 | for(int iy = 1;iy<=60;iy++) | |
333 | if(hBadChMap->GetBinContent(ix,iy)==0){ | |
334 | nDigits+=hDigMap->GetBinContent(ix,iy); | |
335 | nChannels+=1; | |
336 | } | |
337 | //cout<<"nDigits = "<<nDigits<<" ; nChannels = "<<nChannels<<endl; | |
338 | meanOccupancy=nDigits/nChannels; | |
339 | //cout<<"meanOccupancy = "<<meanOccupancy<<endl; | |
340 | //2nd -- spot new bad channels | |
341 | for(int ix = 1;ix<=128;ix++) | |
342 | for(int iy = 1;iy<=60;iy++) | |
343 | if(hBadChMap->GetBinContent(ix,iy)==0&&hDigMap->GetBinContent(ix,iy)>meanOccupancy*10.+1.){ | |
344 | nBadChannelsCurrentIteration++; | |
345 | hBadChMap->Fill(ix-1,iy-1); | |
346 | printf("Noisy channel found! DDL=4 x=%d y=%d\n",ix-1,iy-1); | |
347 | } | |
348 | ||
349 | } | |
350 | cout<<"Total number of noisy channels (DDL = 4): "<<nBadChannelsTotal<<endl; | |
351 | return meanOccupancy; | |
352 | } | |
353 | //============================================================================== | |
354 | void FillDeadMap(TH2* hPedMap, TH2* hDigMap, TH2* hBadChMap){ | |
355 | if(!hPedMap) return; | |
356 | if(!hBadChMap) return; | |
357 | if(!hDigMap)return; | |
358 | Int_t nDeadTotal = 0; | |
359 | for(int ix = 1;ix<=128;ix++) | |
360 | for(int iy = 1;iy<=60;iy++){ | |
361 | if(hPedMap->GetBinContent(ix,iy) < 1. && hBadChMap->GetBinContent(ix,iy)==0) { | |
362 | nDeadTotal++; | |
363 | hBadChMap->Fill(ix-1,iy-1); | |
364 | printf("Dead channel found! DDL=4 x=%d y=%d\n",ix-1,iy-1); | |
365 | } | |
366 | if(hDigMap->GetBinContent(ix,iy) < 1. && hBadChMap->GetBinContent(ix,iy)==0) { | |
367 | nDeadTotal++; | |
368 | hBadChMap->Fill(ix-1,iy-1); | |
369 | printf("Dead channel found! DDL=4 x=%d y=%d\n",ix-1,iy-1); | |
370 | } | |
371 | } | |
372 | cout<<"Total number of noisy channels (DDL = 4): "<<nDeadTotal<<endl; | |
373 | } |