CPV DA for mad channel map calculation (S.Evdokimov)
[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
6Reference run:
7Run Type: PHYSICS
8DA Type: MON
9Number of events needed: ? events
10Input files: thr?_??.dat
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"
19
20//system
21#include <Riostream.h>
22#include <stdlib.h>
23#include <fstream>
24#include <string>
25
26//AliRoot
27#include "AliPHOSCpvRawDigiProducer.h"
28#include "AliPHOSCpvParam.h"
29#include "AliRawReaderDate.h"
30#include "AliBitPacking.h"
31#include "AliPHOSDigit.h"
32#include "AliPHOSGeometry.h"
33
34//ROOT
35#include "TROOT.h"
36#include "TPluginManager.h"
37#include "TSAXParser.h"
38#include "TTree.h"
39#include "TMath.h"
40#include "TString.h"
41#include "TFile.h"
42#include "TSystem.h"
43#include "TKey.h"
44#include "TH2S.h"
45#include "TH2F.h"
46#include "TObject.h"
47#include "TBenchmark.h"
48#include "TMath.h"
49#include "TRandom.h"
50
51void FillBadMap(TH2* DigMap, TH2* BadMap);
52
53int main( int argc, char **argv )
54{
55 gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
56 "*",
57 "TStreamerInfo",
58 "RIO",
59 "TStreamerInfo()");
60 Int_t status,print;
61 Int_t sigcut=3;
62 Bool_t turbo = kTRUE;
63
64 if (argc!=2) {
65 printf("Wrong number of arguments\n");
66 return -1;
67 }
68
69 // log start of process
70 printf("Cpv bad channel map DA program started\n");
71
72 /* report progress */
73 daqDA_progressReport(0);
74
75
76 /* retrieve pedestal tables from DAQ DB */
77 for(int iDDL = 0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL+=2){
78 if(iDDL!=4) continue; // only one module with DDL=4 by now
79 for (int iCC = 0; iCC<AliPHOSCpvParam::kNRows; iCC++){
80 status=daqDA_DB_getFile(Form("thr%d_%02d.dat", iDDL, iCC),Form("thr%d_%02d.dat", iDDL, iCC));
81 if(status!=0) {
82 printf("cannot retrieve file %s from DAQ DB. Exit.\n", Form("thr%d_%02d.dat", iDDL, iCC));
83 //return -1;
84 }
85 }
86 }
87
88 /* connecting to raw data */
89 status=monitorSetDataSource( argv[1] );
90 if (status!=0) {
91 printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
92 return -1;
93 }
94
95 /* declare monitoring program */
96 status=monitorDeclareMp( __FILE__ );
97 if (status!=0) {
98 printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
99 return -1;
100 }
101
102 /* define wait event timeout - 1s max */
103 monitorSetNowait();
104 monitorSetNoWaitNetworkTimeout(1000);
105
106 /* report progress */
107 daqDA_progressReport(5);
108
109
110 // init event counter
111 Int_t iPhysEvnt=0;
112 Int_t iTotEvnt =0;
113
114 // Reader
115 AliRawReader * reader;
116
117 //digiProducer
118 AliPHOSCpvRawDigiProducer* digiProducer = new AliPHOSCpvRawDigiProducer();
119 digiProducer->SetTurbo(turbo);
120 digiProducer->LoadPedFiles();
121 digiProducer->SetCpvMinAmp(0);
122
123 //digits
124 TClonesArray *digits = new TClonesArray("AliPHOSDigit",1);
125 digits->SetName("DIGITS");
126
127 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
128
129 //maps of digits and bad channels
130 TH2F* hMapOfDig[2*AliPHOSCpvParam::kNDDL];
131 TH2I *hBadChMap[2*AliPHOSCpvParam::kNDDL];
132 for(Int_t iDDL=0;iDDL<2*AliPHOSCpvParam::kNDDL;iDDL++){
133 hMapOfDig[iDDL] = new TH2F(Form("hMapOfDig%d",iDDL),Form("Map of digits with substructed pedestals, DDL = %d",iDDL),
134 AliPHOSCpvParam::kPadPcX,0,AliPHOSCpvParam::kPadPcX,
135 AliPHOSCpvParam::kPadPcY,0,AliPHOSCpvParam::kPadPcY);
136 hBadChMap[iDDL] = new TH2I(Form("hBadMap%d",iDDL),Form("Bad Channels Map, DDL= %d",iDDL),
137 AliPHOSCpvParam::kPadPcX,0,AliPHOSCpvParam::kPadPcX,
138 AliPHOSCpvParam::kPadPcY,0,AliPHOSCpvParam::kPadPcY);
139 }
140
141 /* report progress */
142 daqDA_progressReport(10);
143
144 /* main loop (infinite) */
145 for(;;) { // infinite loop
146 struct eventHeaderStruct *event;
147 eventTypeType eventT;
148
149 /* check shutdown condition */
150 if (daqDA_checkShutdown()) {break;}
151
152 // get next event
153 status=monitorGetEventDynamic((void **)&event);
154 if (status==MON_ERR_EOF) { // end of monitoring file has been reached
155 printf("End of monitoring file has been reached! \n");
156 break;
157 }
158
159 if (status!=0) {
160 printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
161 break;
162 }
163
164 // retry if got no event
165 if (event==NULL) continue;
166
167 // use event - here, just write event id to result file
168 eventT=event->eventType;
169 if (eventT==PHYSICS_EVENT) { //we use PHYSICS_EVENT for pedestal not CALIBRATION_EVENT???
170 iTotEvnt++;
171 reader = new AliRawReaderDate((void*)event);
172 digiProducer->LoadNewEvent(reader);
173 digiProducer->MakeDigits(digits);
174 if(digits->GetEntriesFast()>0) iPhysEvnt++;
175 for(Int_t i=0;i<digits->GetEntriesFast();i++)
176 {
177 Int_t relId[4];
178 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
179 geom->AbsToRelNumbering(digit->GetId(),relId);
180 //cout<<"relId : "<<relId[0]<<" "<<relId[1]<<" "<<relId[2]<<" "<<relId[3]<<endl;
181 if(
182 (relId[0]>=0&&relId[0]<=AliPHOSCpvParam::kNDDL)
183 &&(relId[1]==-1)//cpv
184 &&(relId[2]-1>=0&&relId[2]-1<AliPHOSCpvParam::kPadPcX)
185 &&(relId[3]-1>=0&&relId[3]-1<AliPHOSCpvParam::kPadPcY)
186 )
187 hMapOfDig[AliPHOSCpvParam::Mod2DDL(relId[0])]->Fill(relId[2]-1,relId[3]-1);
188 }
189 digits->Clear("C");
190 delete reader;
191 } // if PHYSICS_EVENT
192
193 free(event);
194
195 /* exit when last event received, no need to wait for TERM signal */
196 if (eventT==END_OF_RUN) {
197 printf("EOR event detected\n");
198 break;
199 }
200 }
201
202 Printf(" Received %d events, %d good events",iTotEvnt,iPhysEvnt);
203 /* report progress */
204 daqDA_progressReport(90);
205
206 TFile *fSave = TFile::Open("CpvBadMap.root","RECREATE");
207
208 for(int iDDL = 0; iDDL<2*AliPHOSCpvParam::kNDDL; iDDL+=2){
209 if(hMapOfDig[iDDL]->GetEntries()>0) {
210 FillBadMap(hMapOfDig[iDDL],hBadChMap[iDDL]);
211 fSave->WriteObject(hMapOfDig[iDDL],Form("hMapOfDig%d",iDDL));
212 fSave->WriteObject(hBadChMap[iDDL],Form("hBadChMap%d",iDDL));
213 }
214 }
215 fSave->Close();
216 //status = daqDA_DB_storeFile("CpvBadMap.root","CpvBadMap.root");
217 //if(status) printf("Failed to store CpvBadMap.root in DAQ DB!\n");
218 status = daqDA_FES_storeFile("CpvBadMap.root","CpvBadMap.root");
219 if(status) printf("Failed to store CpvBadMap.root in DAQ FXS!\n");
220
221
222 /* report progress */
223 daqDA_progressReport(100);
224
225
226 return status;
227}
228//==============================================================================
229void FillBadMap(TH2* hDigMap, TH2* hBadChMap){
230 Double_t meanOccupancy = hDigMap->GetEntries()/7680.;
231 Double_t nDigits=0,nChannels=0;
232 int x,y,iterationNumber=1;
233 int nBadChannelsTotal=0,nBadChannelsCurrentIteration=0;
234 //1st iteration
235 //cout<<"Iteration Number = "<<iterationNumber<<endl;
236 for(int ix = 1;ix<=128;ix++)
237 for(int iy = 1;iy<=60;iy++)
238 if(hDigMap->GetBinContent(ix,iy)>meanOccupancy*10.+1.){
239 nBadChannelsCurrentIteration++;
240 hBadChMap->Fill(ix-1,iy-1);
241 printf("Noisy channel found! DDL=4 x=%d y=%d\n",ix-1,iy-1);
242 }
243 //next iterations
244 while(nBadChannelsCurrentIteration!=0){
245 iterationNumber++;
246 //cout<<"Iteration Number = "<<iterationNumber<<endl;
247 nDigits=0;nChannels=0;
248 nBadChannelsTotal+=nBadChannelsCurrentIteration;
249 nBadChannelsCurrentIteration=0;
250 //1st -- calculate new mean occupancy excluding already badly marked channels
251 for(int ix = 1;ix<=128;ix++)
252 for(int iy = 1;iy<=60;iy++)
253 if(hBadChMap->GetBinContent(ix,iy)==0){
254 nDigits+=hDigMap->GetBinContent(ix,iy);
255 nChannels+=1;
256 }
257 //cout<<"nDigits = "<<nDigits<<" ; nChannels = "<<nChannels<<endl;
258 meanOccupancy=nDigits/nChannels;
259 //cout<<"meanOccupancy = "<<meanOccupancy<<endl;
260 //2nd -- spot new bad channels
261 for(int ix = 1;ix<=128;ix++)
262 for(int iy = 1;iy<=60;iy++)
263 if(hBadChMap->GetBinContent(ix,iy)==0&&hDigMap->GetBinContent(ix,iy)>meanOccupancy*10.+1.){
264 nBadChannelsCurrentIteration++;
265 hBadChMap->Fill(ix-1,iy-1);
266 printf("Noisy channel found! DDL=4 x=%d y=%d\n",ix-1,iy-1);
267 }
268
269 }
270 cout<<"Total number of bad channels: "<<nBadChannelsTotal<<endl;
271}
272