9bac53da76b710223519a4018c4dcdc0c8686723
[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: 
7 Run Type:  PHYSICS
8 DA Type: MON
9 Number of events needed: ? events
10 Input files: thr?_??.dat
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
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
51 void FillBadMap(TH2* DigMap, TH2* BadMap);
52
53 int 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 //==============================================================================
229 void 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