fix in calling of gaussian spread function
[u/mrichter/AliRoot.git] / ITS / ITSSPDVertexDiamondda.cxx
1 /*
2 Contact: cvetan.cheshkov@cern.ch
3 Link: http://alisoft.cern.ch/viewvc/trunk/ITS/ITSSPDVertexDiamondda.cxx?root=AliRoot&view=log , /afs/cern.ch/user/c/cheshkov/public/ITS/VD_da_test.date , /afs/cern.ch/user/c/cheshkov/public/08000058338016.30.root.date.gz
4 Reference Run: 58338
5 Run Type: PHYSICS
6 DA Type: MON
7 Number of events needed: 100
8 Input Files: GRP/Geometry/Data , ITS/Align/Data , spd_noisy_ocdb , spd_dead_ocdb , TRIGGER/SPD/PITConditions (all the files are taken from SPD daqDetDB)
9 Output Files: SPDVertexDiamondDA.root
10 Trigger types used: PHYSICS, SPD-F0 
11  */
12
13 #define OUTPUT_FILE "SPDVertexDiamondDA.root"
14
15 extern "C" {
16 #include "daqDA.h"
17 }
18
19 #include "event.h"
20 #include "monitor.h"
21
22 #define ALI_AMORE 
23
24 #ifdef ALI_AMORE
25 #include <AmoreDA.h>
26 //int amore::da::Updated(char const*) {}
27 #endif
28
29 #include <TPluginManager.h>
30 #include <TROOT.h>
31 #include <TH1.h>
32 #include <TH2.h>
33 #include <TPaveText.h>
34 #include <TSystem.h>
35 #include <TGeoGlobalMagField.h>
36 #include <TTimeStamp.h>
37
38 #include "AliLog.h"
39 #include "AliMagF.h"
40 #include "AliRawReaderDate.h"
41 #include "AliCDBManager.h"
42 #include "AliITSMeanVertexer.h"
43 #include <cstdlib>
44
45 int main(int argc, char **argv) {
46
47  gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
48    "*",
49    "TStreamerInfo",
50    "RIO",
51    "TStreamerInfo()"); 
52
53  int status;
54  if (argc<2) {
55   printf("Wrong number of arguments\n");
56   return -1;
57  }
58  // get mean vertex quality cuts
59  status = daqDA_DB_getFile("ITSSPD_VertexQualityTuning_DA.config","./ITSSPD_VertexQualityTuning_DA.config");
60  if (status) {
61   printf("Failed to get config file (ITSSPD_VertexQualityTuning_DA.config) from DAQ DB, status=%d\n", status);
62   return -1;
63  }
64  /* open the config file and retrieve running parameters */
65  Float_t errX  = -1;
66  Float_t r     = -1;
67  UInt_t minClInner = 999 ;
68  UInt_t maxClInner = 999 ;
69  Int_t nEvFirstLoop = 0;
70  Int_t nEvAUTOSAVE = 0; 
71  Int_t zFiducialRegion=0;
72  Int_t timeWindowExport=0; 
73  Int_t timeErrWindowExport=0; 
74
75
76  char name[10][10];
77
78  FILE *fpConfig = fopen("ITSSPD_VertexQualityTuning_DA.config","r");
79  fscanf(fpConfig,"%s %f\n %s %f\n %s %d\n %s %d \n %s %d \n %s %d \n %s %d \n %s %d \n %s %d",&name[0], &errX, &name[1], &r, &name[2], &minClInner,&name[3], &maxClInner, &name[4],&nEvFirstLoop,&name[5],&nEvAUTOSAVE,&name[6],&zFiducialRegion, &name[7], &timeWindowExport, &name[8], &timeErrWindowExport);
80
81  fclose(fpConfig);
82
83  printf("\n\n Mean Vertex quality cuts : \n- errX = %f\n- r = %f\n- minSPD0 = %d maxSPD0 = %d\n- nEventsFirstLoop = %d nEventsAUTOSAVE = %d \n- zFiducialRegion = %d\n- timeWindowExport = %d \n- timeErrWindowExport = %d \n\n\n",errX,r,minClInner,maxClInner,nEvFirstLoop,nEvAUTOSAVE,zFiducialRegion, timeWindowExport, timeErrWindowExport);
84
85  /* define data source : this is argument 1 */  
86  status=monitorSetDataSource( argv[1] );
87  if (status!=0) {
88   printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
89   return -1;
90  }
91
92  /* declare monitoring program */
93  status=monitorDeclareMp( __FILE__ );
94  if (status!=0) {
95   printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
96   return -1;
97  }
98
99  /* define wait event timeout - 1s max */
100  monitorSetNowait();
101  monitorSetNoWaitNetworkTimeout(1000);
102
103  /* log start of process */
104  printf("Vertex-Diamond SPD DA started\n");  
105
106  /* init some counters */
107  int nevents_with_vertex = 0;
108  int nevents_physics=0;
109  int nevents_total=0;
110
111  struct eventHeaderStruct *event;
112  eventTypeType eventT;
113
114  // Get run number
115  if (getenv("DATE_RUN_NUMBER")==0) {
116   printf("DATE_RUN_NUMBER not properly set.\n");
117   return -1;
118  }
119  int runNr = atoi(getenv("DATE_RUN_NUMBER"));
120
121  // Get the necessary OCDB files from the DAQ detector DB
122  if (gSystem->AccessPathName("localOCDB/GRP/Geometry/Data",kFileExists)) {
123   if (gSystem->mkdir("localOCDB/GRP/Geometry/Data",kTRUE) != 0) {
124    printf("Failed to create directory: localOCDB/GRP/Geometry/Data");
125    return -1;
126   }
127  }
128  status = daqDA_DB_getFile("GRP/Geometry/Data","localOCDB/GRP/Geometry/Data/Run0_999999999_v0_s0.root");
129  if (status) {
130   printf("Failed to get geometry file (GRP/Geometry/Data) from DAQdetDB, status=%d\n", status);
131   return -1;
132  }
133
134  if (gSystem->AccessPathName("localOCDB/ITS/Align/Data",kFileExists)) {
135   if (gSystem->mkdir("localOCDB/ITS/Align/Data",kTRUE) != 0) {
136    printf("Failed to create directory: localOCDB/ITS/Align/Data");
137    return -1;
138   }
139  }
140  status = daqDA_DB_getFile("ITS/Align/Data","localOCDB/ITS/Align/Data/Run0_999999999_v0_s0.root");
141  if (status) {
142   printf("Failed to get its-alignment file (ITS/Align/Data) from DAQdetDB, status=%d\n", status);
143   return -1;
144  }
145
146  if (gSystem->AccessPathName("localOCDB/ITS/Calib/SPDNoisy",kFileExists)) {
147   if (gSystem->mkdir("localOCDB/ITS/Calib/SPDNoisy",kTRUE) != 0) {
148    printf("Failed to create directory: localOCDB/ITS/Calib/SPDNoisy");
149    return -1;
150   }
151  }
152  status = daqDA_DB_getFile("spd_noisy_ocdb","localOCDB/ITS/Calib/SPDNoisy/Run0_999999999_v0_s0.root");
153  if (status) {
154   printf("Failed to get spd file (spd_noisy_ocdb) from DAQdetDB, status=%d\n", status);
155   return -1;
156  }
157
158  if (gSystem->AccessPathName("localOCDB/ITS/Calib/SPDDead",kFileExists)) {
159   if (gSystem->mkdir("localOCDB/ITS/Calib/SPDDead",kTRUE) != 0) {
160    printf("Failed to create directory: localOCDB/ITS/Calib/SPDDead");
161    return -1;
162   }
163  }
164  status = daqDA_DB_getFile("spd_dead_ocdb","localOCDB/ITS/Calib/SPDDead/Run0_999999999_v0_s0.root");
165  if (status) {
166   printf("Failed to get spd file (spd_dead_ocdb) from DAQdetDB, status=%d\n", status);
167   return -1;
168  }
169
170  if (gSystem->AccessPathName("localOCDB/ITS/Calib/SPDSparseDead",kFileExists)) {
171   if (gSystem->mkdir("localOCDB/ITS/Calib/SPDSparseDead",kTRUE) != 0) {
172    printf("Failed to create directory: localOCDB/ITS/Calib/SPDSparseDead");
173    return -1;
174   }
175  }
176
177  status = daqDA_DB_getFile("spd_sparsedead_ocdb","localOCDB/ITS/Calib/SPDSparseDead/Run0_999999999_v0_s0.root");
178  if (status) {
179   printf("Failed to get spd file (spd_sparsedead_ocdb) from DAQdetDB, status=%d\n", status);
180   return -1;
181  }
182
183  if (gSystem->AccessPathName("localOCDB/TRIGGER/SPD/PITConditions",kFileExists)) {
184   if (gSystem->mkdir("localOCDB/TRIGGER/SPD/PITConditions",kTRUE) != 0) {
185    printf("Failed to create directory: localOCDB/TRIGGER/SPD/PITConditions");
186    return -1;
187   }
188  }
189  status = daqDA_DB_getFile("TRIGGER/SPD/PITConditions","localOCDB/TRIGGER/SPD/PITConditions/Run0_999999999_v0_s0.root");
190  if (status) {
191   printf("Failed to get spd trigger file (TRIGGER/SPD/PITConditions) from DAQdetDB, status=%d\n", status);
192   return -1;
193  }
194
195  status = daqDA_DB_getFile("mfchebKGI_sym.root","localOCDB/mfchebKGI_sym.root");
196  if (status) {
197   printf("Failed to get spd file (mfchebKGI_sym.root) from DAQdetDB, status=%d\n", status);
198   return -1;
199  }
200
201  // Global initializations
202
203  // The B filed is required in AliITSClusterFinderV2SPD
204  // for the Lorentz angle correction. B set to 0.      
205  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 0., 0., AliMagF::k5kGUniform,AliMagF::kBeamTypeAA,-1,2,15,"localOCDB/mfchebKGI_sym.root"));  
206  AliLog::SetGlobalLogLevel(AliLog::kError);
207  AliCDBManager *man = AliCDBManager::Instance();
208  man->SetDefaultStorage("local://localOCDB");
209  man->SetRun(runNr);
210
211  // Init mean vertexer
212  AliITSMeanVertexer *mv = new AliITSMeanVertexer();
213
214  // Optionally change the Z fiducial region here (default +/- 40 cm
215   mv->SetZFiducialRegion(zFiducialRegion);
216
217  if (!mv->Init()) {
218   printf("Initialization of mean vertexer object failed ! Check the log for details");
219   return -1;
220  }
221
222  mv->SetCutOnErrX(errX);
223  mv->SetCutOnR(r);
224  mv->SetCutOnCls(minClInner,maxClInner);
225
226  gSystem->Exec(Form("rm %s",OUTPUT_FILE));
227
228  TTimeStamp *timeStamp;
229  timeStamp->Set();
230  Int_t t1 = timeStamp->GetSec(); 
231  Int_t t2=0; 
232
233 // Initialization of AMORE sender
234 #ifdef ALI_AMORE
235  amore::da::AmoreDA vtxAmore(amore::da::AmoreDA::kSender);
236 #endif
237  Int_t amore_status = 0;
238  /* main loop (infinite) */
239  for(;;) {
240
241   /* check shutdown condition */
242   if (daqDA_checkShutdown()) {break;}
243
244   /* get next event (blocking call until timeout) */
245   status=monitorGetEventDynamic((void **)&event);
246   if (status==MON_ERR_EOF) {
247    printf ("End of File detected\n");
248    break; /* end of monitoring file has been reached */
249   }
250
251   if (status!=0) {
252    printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
253    break;
254   }
255
256   /* retry if got no event */
257   if (event==NULL) {
258    continue;
259   }
260
261   nevents_total++;
262   eventT=event->eventType;
263   switch (event->eventType){
264
265    /* START OF RUN */
266    case START_OF_RUN:
267     break;
268     /* END START OF RUN */
269
270     /* END OF RUN */
271    case END_OF_RUN:
272     break;
273
274    case PHYSICS_EVENT:
275     nevents_physics++;
276     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
277
278     // Run mean-vertexer reco
279     if (mv->Reconstruct(rawReader)) nevents_with_vertex++;
280     if(nevents_physics < nEvFirstLoop) continue;
281     // Auto save
282     if ((nevents_physics%nEvAUTOSAVE) == 0) {
283      TH2F *histo = ((TH2F*)mv->GetVertexXY());
284      histo->SetStats(kFALSE);
285      histo->SetTitle("");
286      histo->GetListOfFunctions()->SetOwner(kTRUE);
287      if(histo->GetListOfFunctions()->GetEntries()<1) histo->GetListOfFunctions()->Add(new TPaveText(-5,4.5,5.,6.2,"br"));
288      for(Int_t i=0; i<histo->GetListOfFunctions()->GetEntries(); i++){
289       TString funcName = histo->GetListOfFunctions()->At(i)->ClassName();
290       if(funcName.Contains("TPaveText")){
291        TPaveText *p = (TPaveText*)histo->GetListOfFunctions()->At(i);
292        p->Clear();
293        p->AddText(Form("%f events with vertex (%i out of %i processed events)",((Double_t)nevents_with_vertex)/((Double_t)nevents_physics),nevents_with_vertex,nevents_physics));
294        p->AddText(Form("%f events with good vertex (%i out of %i events with vertex)",histo->GetEntries()/((Double_t)nevents_with_vertex),(Int_t)histo->GetEntries(),nevents_with_vertex));
295       }
296      }
297      mv->WriteVertices(OUTPUT_FILE);
298
299      timeStamp->Set();
300      t2 = timeStamp->GetSec();
301      
302      
303     
304 #ifdef ALI_AMORE
305
306      if (TMath::Abs((t2-t1)-timeWindowExport) < timeErrWindowExport){
307        t1=t2;
308        // send the histos to AMORE pool
309        amore_status=vtxAmore.Send(mv->GetVertexXY()->GetName(),mv->GetVertexXY());
310        if(amore_status) printf("AMORE XY send status: %d\n",amore_status);
311        TH1D *hProj = mv->GetVertexXY()->ProjectionX();
312        amore_status=vtxAmore.Send(Form("%s_ProjX",mv->GetVertexXY()->GetName()),hProj); 
313        if(amore_status) printf("AMORE X send status: %d\n",amore_status);
314        if(hProj) delete hProj;
315        hProj = mv->GetVertexXY()->ProjectionY();
316        amore_status=vtxAmore.Send(Form("%s_ProjY",mv->GetVertexXY()->GetName()),hProj); 
317        if(amore_status) printf("AMORE Y send status: %d\n",amore_status);
318        if(hProj) hProj->Delete();
319        amore_status=vtxAmore.Send(mv->GetVertexZ()->GetName(),mv->GetVertexZ());
320        if(amore_status) printf("AMORE Z  send status: %d\n",amore_status);
321      }
322      if (TMath::Abs(t2-t1)> timeWindowExport) t1=t2;
323
324 #endif
325     }
326
327     delete rawReader;
328   }
329
330   /* free resources */
331   free(event);
332
333   /* exit when last event received, no need to wait for TERM signal */
334   if (eventT==END_OF_RUN) {
335    printf("EOR event detected\n");
336    break;
337   }
338  }
339
340  TH2F *histo = ((TH2F*)mv->GetVertexXY());
341  histo->SetStats(kFALSE);
342  histo->SetTitle("");
343  histo->GetListOfFunctions()->SetOwner(kTRUE);
344  if(histo->GetListOfFunctions()->GetEntries()<1) histo->GetListOfFunctions()->Add(new TPaveText(-5,4.5,5.,6.2,"br"));
345  for(Int_t i=0; i<histo->GetListOfFunctions()->GetEntries(); i++){
346   TString funcName = histo->GetListOfFunctions()->At(i)->ClassName();
347   if(funcName.Contains("TPaveText")){
348    TPaveText *p = (TPaveText*)histo->GetListOfFunctions()->At(i);
349    p->Clear(); 
350    p->AddText(Form("%f events with vertex (%i out of %i processed events)",((Double_t)nevents_with_vertex)/((Double_t)nevents_physics),nevents_with_vertex,nevents_physics));
351    p->AddText(Form("%f events with good vertex (%i out of %i events with vertex)",histo->GetEntries()/((Double_t)nevents_with_vertex),(Int_t)histo->GetEntries(),nevents_with_vertex));
352   }
353  }
354  mv->WriteVertices(OUTPUT_FILE);
355
356 #ifdef ALI_AMORE
357  if (TMath::Abs((t2-t1)-timeWindowExport) < timeErrWindowExport){
358        t1=t2;
359        // send the histos to AMORE pool
360        amore_status=vtxAmore.Send(mv->GetVertexXY()->GetName(),mv->GetVertexXY());  
361        if(amore_status) printf("AMORE XY send status: %d\n",amore_status);
362        TH1D *hProj = mv->GetVertexXY()->ProjectionX();
363        amore_status=vtxAmore.Send(Form("%s_ProjX",mv->GetVertexXY()->GetName()),hProj);
364        if(amore_status) printf("AMORE X send status: %d\n",amore_status);
365        if(hProj) delete hProj;
366        hProj = mv->GetVertexXY()->ProjectionY();
367        amore_status=vtxAmore.Send(Form("%s_ProjY",mv->GetVertexXY()->GetName()),hProj);
368        if(amore_status) printf("AMORE Y send status: %d\n",amore_status);
369        if(hProj) hProj->Delete();
370        amore_status=vtxAmore.Send(mv->GetVertexZ()->GetName(),mv->GetVertexZ());
371        if(amore_status) printf("AMORE Z  send status: %d\n",amore_status);
372  }
373  if (TMath::Abs(t2-t1)> timeWindowExport) t1=t2;
374
375 #endif
376        
377  delete mv;
378
379  /* write report */
380  printf("Run #%s, received %d events with vertex, out of %d physics and out of %d total events\n",getenv("DATE_RUN_NUMBER"),nevents_with_vertex,nevents_physics,nevents_total);
381
382  status=0;
383
384  /* export file to FXS */
385  if (daqDA_FES_storeFile(OUTPUT_FILE, "VertexDiamond")) {
386   status=-2;
387  }
388
389  return status;
390 }