PHOS/DA/CMakeLists.txt
[u/mrichter/AliRoot.git] / PHOS / DA / PHOSCPVBCMda.cxx
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 }