minor coverity defect: adding self-assignment protection
[u/mrichter/AliRoot.git] / T0 / T0Physda.cxx
1 /*
2 T0 DA for online calibration
3 Contact: AllaMaevskaya@cern.ch
4 Run Type: PHYSICS
5 DA Type: MON
6 Number of events needed: 10000 
7 Input Files: inPhys.dat, external parameters, T0/Calib/Slewing_Walk
8 Output Files: daPhys.root, to be exported to the DAQ FXS
9 Trigger types used: PHYSICS_EVENT
10 ------------------------------Alla
11 Now trigger type changed to CALICRATION_EVENT
12 to have data during test.
13 SOULD BE CHANGED BACK BEFORE BEAM
14 ------------------------------- Alla
15 */
16
17 #define FILE_OUT "daPhys.root"
18 #define FILE_IN "inPhys.dat"
19 #include <daqDA.h>
20 #include <event.h>
21 #include <monitor.h>
22  
23 #include <Riostream.h>
24 #include <stdio.h>
25 #include <stdlib.h>
26
27 //AliRoot
28 #include <AliRawReaderDate.h>
29 #include <AliRawReader.h>
30 #include <AliT0RawReader.h>
31 #include <AliT0CalibWalk.h>
32 #include <AliCDBManager.h>
33 #include <AliCDBEntry.h>
34 //ROOT
35 #include "TROOT.h"
36 #include "TPluginManager.h"
37 #include "TFile.h"
38 #include "TKey.h"
39 #include "TObject.h"
40 #include "TBenchmark.h"
41 #include "TString.h"
42 #include "TH1.h"
43 #include "TMath.h"
44
45
46 /* Main routine
47       Arguments: 
48       1- monitoring data source
49 */
50 int main(int argc, char **argv) {
51 //int main(){
52   int status;
53
54   /* magic line */
55   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
56                                         "*",
57                                         "TStreamerInfo",
58                                         "RIO",
59                                         "TStreamerInfo()");
60   
61   if(daqDA_DB_getFile(FILE_IN, FILE_IN)){
62     printf("Couldn't get input file >>inPhys.dat<< from DAQ_DB !!!\n");
63     return -1;
64   }
65   
66   
67   FILE *inp;
68   char c;
69   inp = fopen(FILE_IN, "r");
70   if(!inp){
71     printf("Input file >>inPhys.dat<< not found !!!\n");
72     return -1;
73   }
74   int  kcbx, kt0bx, knpmtA, knpmtC;
75   float kclx,kcmx, kt0lx, kt0hx;
76   
77   while((c=getc(inp))!=EOF) {
78     switch(c) {
79     case 'a': {fscanf(inp, "%d", &kcbx ); break;} //N of X bins hCFD1_CFD
80     case 'b': {fscanf(inp, "%f", &kclx ); break;} //Low x hCFD1_CFD
81     case 'c': {fscanf(inp, "%f", &kcmx ); break;} //High x hCFD1_CF
82     case 'd': {fscanf(inp, "%d", &knpmtC ); break;} //number of reference PMTC
83     case 'e': {fscanf(inp, "%d", &knpmtA ); break;} //number of reference PMTA
84     case 'f': {fscanf(inp, "%d", &kt0bx ); break;} //N of X bins hT0
85     case 'g': {fscanf(inp, "%f", &kt0lx ); break;} //Low x hT0
86     case 'k': {fscanf(inp, "%f", &kt0hx ); break;} //High x hT0
87     }
88   }
89   fclose(inp);
90
91   if (argc!=2) {
92     printf("Wrong number of arguments\n");
93     return -1;
94   }
95   
96
97   /* define data source : this is argument 1 */  
98   status=monitorSetDataSource( argv[1] );
99   if (status!=0) {
100     printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
101     return -1;
102   }
103   
104   
105   /* declare monitoring program */
106   status=monitorDeclareMp( __FILE__ );
107   if (status!=0) {
108     printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
109     return -1;
110   }
111   
112   
113   /* define wait event timeout - 1s max */
114   monitorSetNowait();
115   monitorSetNoWaitNetworkTimeout(1000);
116   
117   
118   /* log start of process */
119   printf("T0 monitoring program started\n");  
120   
121   // Get run number
122   if (getenv("DATE_RUN_NUMBER")==0) {
123     printf("DATE_RUN_NUMBER not properly set.\n");
124     return -1;
125   }
126   int runNr = atoi(getenv("DATE_RUN_NUMBER"));
127
128  // Get the necessary OCDB files from the DAQ detector DB
129   if (gSystem->AccessPathName("localOCDB/T0/Calib/Slewing_Walk/",kFileExists)) {
130     if (gSystem->mkdir("localOCDB/T0/Calib/Slewing_Walk/",kTRUE) != 0) {
131       printf("Failed to create directory: localOCDB/T0/Calib/Slewing_Walk/");
132       return -1;
133     }
134   }
135
136   status = daqDA_DB_getFile("T0/Calib/Slewing_Walk","localOCDB/T0/Calib/Slewing_Walk/Run0_999999999_v0_s0.root");
137   if (status) {
138     printf("Failed to get geometry file (GRP/Geometry/Data) from DAQdetDB, status=%d\n", status);
139     return -1;
140   }
141   TGraph *gr[24]; TGraph *gramp[24];
142   AliCDBManager *man = AliCDBManager::Instance();
143   man->SetDefaultStorage("local://localOCDB");
144   man->SetRun(runNr);
145   AliCDBEntry *entry = AliCDBManager::Instance()->Get("T0/Calib/Slewing_Walk");
146   if(entry) {
147     AliT0CalibWalk *fParam = (AliT0CalibWalk*)entry->GetObject();
148     for (Int_t i=0; i<24; i++) {
149       gr[i] = fParam->GetWalk(i); 
150       gramp[i] = fParam->GetQTC(i); 
151     }
152   }
153   Int_t chargeQT0[24], chargeQT1[24];
154   Float_t adc ,walk, amp;
155  
156   // Allocation of histograms - start
157
158   TH1F *hCFD1minCFD[24];  
159   TH1F *hCFD[24];  
160    
161   for(Int_t ic=0; ic<24; ic++) {
162     hCFD1minCFD[ic] = new TH1F(Form("CFD1minCFD%d",ic+1),"CFD-CFD",kcbx,kclx,kcmx);
163     hCFD[ic] = new TH1F(Form("CFD%d",ic+1),"CFD",kt0bx,kt0lx,kt0hx);
164   }
165   TH1F *hVertex = new TH1F("hVertex","T0 time",kt0bx,kt0lx,kt0hx);
166   
167    // Allocation of histograms - end
168
169
170  Int_t iev=0;
171   /* main loop (infinite) */
172   for(;;) {
173     struct eventHeaderStruct *event;
174     eventTypeType eventT;
175     
176     /* check shutdown condition */
177     if (daqDA_checkShutdown()) {break;}
178     
179     /* get next event (blocking call until timeout) */
180     status=monitorGetEventDynamic((void **)&event);
181     if (status==(int)MON_ERR_EOF) {
182       printf ("End of File detected\n");
183       break; /* end of monitoring file has been reached */
184     }
185     
186     if (status!=0) {
187       printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
188       break;
189     }
190     
191     /* retry if got no event */
192     if (event==NULL) {
193       continue;
194     }
195     
196     /* use event - here, just write event id to result file */
197     eventT=event->eventType;
198     
199     switch (event->eventType){
200       
201     case START_OF_RUN:
202       break;
203         
204     case END_OF_RUN:
205       break;
206       
207     case PHYSICS_EVENT:
208           //    case CALIBRATION_EVENT:
209       iev++;
210       
211       if(iev==1){
212         printf("First event - %i\n",iev);
213       }
214       
215       // Initalize raw-data reading and decoding
216       AliRawReader *reader = new AliRawReaderDate((void*)event);
217       
218       // Enable the following two lines in case of real-data
219       reader->RequireHeader(kTRUE);
220       AliT0RawReader *start = new AliT0RawReader(reader, kTRUE);
221       
222       // Read raw data
223       Int_t allData[110][5];
224       for(Int_t i0=0;i0<107;i0++)
225         for(Int_t j0=0;j0<5;j0++)
226           allData[i0][j0] = 0;
227       
228       if(start->Next()){
229         for (Int_t i=0; i<107; i++) {
230           for(Int_t iHit=0;iHit<5;iHit++){
231             allData[i][iHit]= start->GetData(i,iHit);
232           }
233         }
234       }
235       
236       // Fill the histograms
237       walk = adc = amp = -999;
238       for (Int_t in=0; in<12;  in++)
239         {
240           chargeQT0[in]=allData[2*in+25][0];
241           chargeQT1[in]=allData[2*in+26][0];
242         }       
243       for (Int_t in=12; in<24;  in++)
244         {
245           chargeQT0[in]=allData[2*in+57][0];
246           chargeQT1[in]=allData[2*in+58][0];
247         }
248      Float_t besttimeA=9999999;
249       Float_t besttimeC=9999999;
250       Float_t time[24]; 
251        Float_t meanShift[24];
252        for (Int_t ik = 0; ik<24; ik++)
253          {       
254            if( ( chargeQT0[ik] - chargeQT1[ik])>0)  {
255              adc = chargeQT0[ik] - chargeQT1[ik];
256              // cout<<ik <<"  "<<adc<<endl;
257            }
258            if(gramp[ik])
259              amp = gramp[ik]->Eval(Double_t(adc));
260            if(amp < 0.8) continue;
261            if(gr[ik]) 
262              walk = Int_t(gr[ik]->Eval(Double_t(adc) ) );
263            
264            if(ik<12 && allData[ik+1][0]>0 && allData[knpmtC][0]>0 ){
265              hCFD1minCFD[ik]->Fill(allData[ik+1][0]-allData[knpmtC][0]);
266              if( walk >-100) hCFD[ik]->Fill(allData[ik+1][0] - walk);
267              // cout<<ik<<" "<<allData[ik+1][0]<<" adc "<<adc<<" walk "<<walk<<endl;
268            }
269            
270            if(ik>11 && allData[ik+45][0]>0 && allData[56+knpmtA][0]>0 )
271              {
272              hCFD1minCFD[ik]->Fill(allData[ik+45][0]-allData[56+knpmtA][0]);
273               if( walk >-100) hCFD[ik]->Fill(allData[ik+45][0] - walk);
274               // cout<<ik<<" "<<allData[ik+1][0]<<" adc "<<adc<<" walk "<<walk<<endl;
275              }
276            if(iev == 10000) {   
277              meanShift[ik] =  hCFD1minCFD[ik]->GetMean(); 
278              if(ik==knpmtC || ik==(56+knpmtA)) meanShift[ik]=0;
279            }
280          }
281       //fill  mean time _ fast reconstruction
282       if (iev > 10000 )
283         {
284           for (Int_t in=0; in<12; in++)  
285             {
286               time[in] = allData[in+1][0] - meanShift[in]  ;
287               time[in+12] = allData[in+56+1][0] ;
288             }
289           for (Int_t ipmt=0; ipmt<12; ipmt++){
290             if(time[ipmt] > 1 ) {
291               if(time[ipmt]<besttimeC)
292                 besttimeC=time[ipmt]; //timeC
293             }
294           }
295            for ( Int_t ipmt=12; ipmt<24; ipmt++){
296              if(time[ipmt] > 1) {
297                if(time[ipmt]<besttimeA) 
298                  besttimeA=time[ipmt]; //timeA
299              }
300            }
301            if(besttimeA<9999999 &&besttimeC< 9999999) {
302              Float_t t0 =0.001* 24.4 * Float_t( besttimeA+besttimeC)/2.;
303              hVertex->Fill(t0);
304            }
305         }
306            
307      delete start;
308      start = 0x0;
309      delete reader;
310      reader= 0x0;
311       // End of fill histograms
312       
313     }
314     
315     /* free resources */
316     free(event);
317     
318     /* exit when last event received, no need to wait for TERM signal */
319     if (eventT==END_OF_RUN) {
320       printf("EOR event detected\n");
321       printf("Number of events processed - %i\n ",iev);         
322       break;
323     }
324   }
325   printf("After loop, before writing histos\n");
326   // write a file with the histograms
327
328   TFile hist(FILE_OUT,"RECREATE");
329
330   for(Int_t j=0;j<24;j++){
331     hCFD1minCFD[j]->SetDirectory(&hist);
332     hCFD1minCFD[j]->Write();
333     hCFD[j]->Write();
334
335   }
336   hVertex->Write();
337   hist.Close();
338   //delete hist;
339
340   status=0;
341
342   /* export file to FXS */
343   if (daqDA_FES_storeFile(FILE_OUT, "PHYSICS")) {
344     status=-2;
345   }
346
347   return status;
348 }
349
350