CMake: Renaming ADPEDESTAL to AD0PEDESTAL to match DAQ setting
[u/mrichter/AliRoot.git] / PHOS / DA / PHOSCPVBCMda.cxx
CommitLineData
115e0859 1/*
2CPV BCM DA for processing physics runs and producing bad channel map for further use at reconstruction time.
3
4Contact: Sergey Evdokimov <sevdokim@cern.ch>
5Link: https://twiki.cern.ch/twiki/bin/view/ALICE/CPVda
dfd5ff54 6Reference run: 214340 (/afs/cern.ch/user/s/sevdokim/public/CPV_run214340_standalone.raw)
115e0859 7Run Type: PHYSICS
8DA Type: MON
dfd5ff54 9Number of events needed: ~100000 events
10Input files: thr?_??.dat CpvPeds.root PHOSCPVBCMda.cfg
115e0859 11Output files: CpvBadMap.root
12Trigger types used: PHYSICS_EVENT
13*/
14
15//daqDA
16#include "event.h"
17#include "monitor.h"
18#include "daqDA.h"
dfd5ff54 19//AMORE monitoring framework
20#include <AmoreDA.h>
115e0859 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
dfd5ff54 53Double_t FillNoisyMap(TH2* DigMap, TH2* BadMap); //returns mean occupancy when all bad channels are excluded
54void FillDeadMap(TH2* PedMap, TH2* DigMap, TH2* BadMap);
115e0859 55
56int main( int argc, char **argv )
57{
58 gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
59 "*",
60 "TStreamerInfo",
61 "RIO",
62 "TStreamerInfo()");
dfd5ff54 63 Int_t status,statusPeds=0,print,statusBadMap=0;
115e0859 64 Int_t sigcut=3;
dfd5ff54 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
115e0859 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
dfd5ff54 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
115e0859 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));
dfd5ff54 100 //return 11;//error code 11 (cannot retrive thr.dat from DAQ DB)
115e0859 101 }
102 }
103 }
104
dfd5ff54 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
115e0859 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();
dfd5ff54 148 digiProducer->SetCpvMinAmp(minAmpl);
115e0859 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];
dfd5ff54 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;
115e0859 163 }
164
dfd5ff54 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 }
115e0859 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
dfd5ff54 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)
115e0859 262 TFile *fSave = TFile::Open("CpvBadMap.root","RECREATE");
263
dfd5ff54 264 Double_t meanOccupancy = 0;
265 amore::da::AmoreDA* myAmore = new amore::da::AmoreDA(amore::da::AmoreDA::kSender);
115e0859 266 for(int iDDL = 0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL+=2){
267 if(hMapOfDig[iDDL]->GetEntries()>0) {
dfd5ff54 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]);
115e0859 272 fSave->WriteObject(hMapOfDig[iDDL],Form("hMapOfDig%d",iDDL));
dfd5ff54 273 //send digit maps to amore
274 myAmore->Send(Form("hMapOfDig%d",iDDL),hMapOfDig[iDDL]);
115e0859 275 fSave->WriteObject(hBadChMap[iDDL],Form("hBadChMap%d",iDDL));
276 }
277 }
278 fSave->Close();
dfd5ff54 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");
115e0859 297
dfd5ff54 298 }
115e0859 299
300 /* report progress */
301 daqDA_progressReport(100);
302
303
304 return status;
305}
306//==============================================================================
dfd5ff54 307Double_t FillNoisyMap(TH2* hDigMap, TH2* hBadChMap){
308 if(!hDigMap)return 0;
309 if(!hBadChMap)return 0;
115e0859 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 }
dfd5ff54 350 cout<<"Total number of noisy channels (DDL = 4): "<<nBadChannelsTotal<<endl;
351 return meanOccupancy;
352}
353//==============================================================================
354void 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;
115e0859 373}