]>
Commit | Line | Data |
---|---|---|
eb35e591 | 1 | /* |
2 | Contact: cvetan.cheshkov@cern.ch | |
3 | Link: missing | |
4 | Run Type: PHYSICS | |
5 | DA Type: MON | |
6 | Number of events needed: 10000 | |
7 | Input Files: | |
8 | Output Files: | |
9 | Trigger types used: PHYSICS | |
10 | */ | |
11 | ||
12 | #define X_LIMIT 2.0 | |
13 | #define Y_LIMIT 2.0 | |
14 | #define Z_LIMIT 50.0 | |
15 | #define X_DELTA 0.02 | |
16 | #define Y_DELTA 0.02 | |
17 | #define Z_DELTA 0.2 | |
18 | #define OUTPUT_FILE "SPDVertexDiamondDA.root" | |
19 | #define CDB_STORAGE "local://$ALICE_ROOT" | |
20 | #define N_EVENTS_AUTOSAVE 50 | |
21 | #define NFITPARAMS 5 | |
22 | ||
23 | extern "C" { | |
24 | #include "daqDA.h" | |
25 | } | |
26 | ||
27 | #include "event.h" | |
28 | #include "monitor.h" | |
29 | ||
55c5e86d | 30 | #ifdef ALI_AMORE |
31 | #include <AmoreDA.h> | |
308c2f7c | 32 | //int amore::da::Updated(char const*) {} |
55c5e86d | 33 | #endif |
34 | ||
eb35e591 | 35 | #include <TTree.h> |
36 | #include <TH1.h> | |
37 | #include <TH2.h> | |
38 | #include <TF2.h> | |
39 | #include <TFile.h> | |
40 | #include <TPluginManager.h> | |
41 | #include <TROOT.h> | |
42 | #include <TFitter.h> | |
43 | ||
44 | #include "AliRawReaderDate.h" | |
45 | #include "AliGeomManager.h" | |
46 | #include "AliCDBManager.h" | |
47 | #include "AliESDVertex.h" | |
48 | #include "AliITSDetTypeRec.h" | |
49 | #include "AliITSInitGeometry.h" | |
50 | #include "AliITSVertexer3DTapan.h" | |
51 | ||
52 | TF2 *fitFcn = 0x0; | |
53 | ||
54 | AliESDVertex* FitVertexDiamond(TH2F *hXY, TH1F *hZ); | |
55 | ||
56 | int main(int argc, char **argv) { | |
57 | ||
58 | gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo", | |
59 | "*", | |
60 | "TStreamerInfo", | |
61 | "RIO", | |
62 | "TStreamerInfo()"); | |
63 | TFitter *minuitFit = new TFitter(NFITPARAMS); | |
64 | TVirtualFitter::SetFitter(minuitFit); | |
65 | ||
66 | int status; | |
67 | if (argc<2) { | |
68 | printf("Wrong number of arguments\n"); | |
69 | return -1; | |
70 | } | |
71 | ||
72 | /* define data source : this is argument 1 */ | |
73 | status=monitorSetDataSource( argv[1] ); | |
74 | if (status!=0) { | |
75 | printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status)); | |
76 | return -1; | |
77 | } | |
78 | ||
79 | /* declare monitoring program */ | |
80 | status=monitorDeclareMp( __FILE__ ); | |
81 | if (status!=0) { | |
82 | printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status)); | |
83 | return -1; | |
84 | } | |
85 | ||
86 | /* define wait event timeout - 1s max */ | |
87 | monitorSetNowait(); | |
88 | monitorSetNoWaitNetworkTimeout(1000); | |
89 | ||
90 | /* log start of process */ | |
91 | printf("Vertex-Diamond SPD DA started\n"); | |
92 | ||
93 | /* init some counters */ | |
94 | int nevents_with_vertex = 0; | |
95 | int nevents_physics=0; | |
96 | int nevents_total=0; | |
97 | ||
98 | struct eventHeaderStruct *event; | |
99 | eventTypeType eventT; | |
100 | ||
101 | // Get run number | |
102 | if (getenv("DATE_RUN_NUMBER")==0) { | |
103 | printf("DATE_RUN_NUMBER not properly set.\n"); | |
104 | return -1; | |
105 | } | |
106 | int runNr = atoi(getenv("DATE_RUN_NUMBER")); | |
107 | ||
108 | // Histograms initialization | |
109 | TH2F *hXY = new TH2F("hXY","Vertex Diamond (Y vs X)", | |
110 | 2*(Int_t)(X_LIMIT/X_DELTA),-X_LIMIT,X_LIMIT, | |
111 | 2*(Int_t)(Y_LIMIT/Y_DELTA),-Y_LIMIT,Y_LIMIT); | |
112 | TH1F *hZ = new TH1F("hZ"," Longitudinal Vertex Profile", | |
113 | 2*(Int_t)(Z_LIMIT/Z_DELTA),-Z_LIMIT,Z_LIMIT); | |
114 | ||
115 | // Global initializations | |
116 | AliCDBManager *man = AliCDBManager::Instance(); | |
117 | man->SetDefaultStorage(CDB_STORAGE); | |
118 | man->SetRun(runNr); | |
119 | AliGeomManager::LoadGeometry("geometry.root"); | |
120 | AliGeomManager::ApplyAlignObjsFromCDB("ITS"); | |
121 | ||
122 | // ITS initializations | |
123 | AliITSInitGeometry initgeom; | |
124 | AliITSgeom *geom = initgeom.CreateAliITSgeom(); | |
125 | printf("Geometry name: %s\n",(initgeom.GetGeometryName()).Data()); | |
126 | ||
127 | AliITSDetTypeRec *detTypeRec = new AliITSDetTypeRec(); | |
128 | detTypeRec->SetITSgeom(geom); | |
129 | detTypeRec->SetDefaults(); | |
130 | detTypeRec->SetDefaultClusterFindersV2(kTRUE); | |
131 | ||
55c5e86d | 132 | // Initialization of AMORE sender |
133 | #ifdef ALI_AMORE | |
134 | amore::da::AmoreDA vtxAmore(amore::da::AmoreDA::kSender); | |
135 | #endif | |
eb35e591 | 136 | /* main loop (infinite) */ |
137 | for(;;) { | |
138 | ||
139 | /* check shutdown condition */ | |
140 | if (daqDA_checkShutdown()) {break;} | |
141 | ||
142 | /* get next event (blocking call until timeout) */ | |
143 | status=monitorGetEventDynamic((void **)&event); | |
144 | if (status==MON_ERR_EOF) { | |
145 | printf ("End of File detected\n"); | |
146 | break; /* end of monitoring file has been reached */ | |
147 | } | |
148 | ||
149 | if (status!=0) { | |
150 | printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status)); | |
151 | break; | |
152 | } | |
153 | ||
154 | /* retry if got no event */ | |
155 | if (event==NULL) { | |
156 | continue; | |
157 | } | |
158 | ||
159 | nevents_total++; | |
160 | eventT=event->eventType; | |
161 | switch (event->eventType){ | |
162 | ||
163 | /* START OF RUN */ | |
164 | case START_OF_RUN: | |
165 | break; | |
166 | /* END START OF RUN */ | |
167 | ||
168 | /* END OF RUN */ | |
169 | case END_OF_RUN: | |
170 | break; | |
171 | ||
172 | case PHYSICS_EVENT: | |
173 | nevents_physics++; | |
174 | AliRawReader *rawReader = new AliRawReaderDate((void*)event); | |
175 | ||
176 | // Run SPD cluster finder | |
177 | TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree | |
178 | detTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD"); | |
179 | ||
180 | // Run vertex-finder | |
181 | AliITSVertexer3DTapan *vertexer = new AliITSVertexer3DTapan(geom,1000); | |
308c2f7c | 182 | AliESDVertex *vtx = vertexer->FindVertexForCurrentEvent(clustersTree); |
eb35e591 | 183 | |
184 | if (TMath::Abs(vtx->GetChi2()) < 0.1) { | |
185 | // Fill the vertex into the histos | |
186 | nevents_with_vertex++; | |
187 | hXY->Fill(vtx->GetXv(),vtx->GetYv()); | |
188 | hZ->Fill(vtx->GetZv()); | |
189 | ||
190 | // Auto save | |
191 | if ((nevents_with_vertex%N_EVENTS_AUTOSAVE) == 0) { | |
192 | TFile outFile(OUTPUT_FILE, "update"); | |
193 | AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ); | |
194 | if (fitVtx) { | |
195 | fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite); | |
196 | TH1 *fithXY = fitFcn->CreateHistogram(); | |
197 | fithXY->Write(fithXY->GetName(),TObject::kOverwrite); | |
198 | delete fithXY; | |
199 | } | |
200 | hXY->Write(hXY->GetName(),TObject::kOverwrite); | |
201 | hZ->Write(hZ->GetName(),TObject::kOverwrite); | |
202 | outFile.Close(); | |
203 | delete fitVtx; | |
204 | } | |
205 | } | |
206 | ||
207 | delete vtx; | |
208 | delete vertexer; | |
209 | ||
210 | delete clustersTree; | |
211 | delete rawReader; | |
212 | ||
213 | } | |
214 | ||
215 | /* free resources */ | |
216 | free(event); | |
217 | ||
218 | /* exit when last event received, no need to wait for TERM signal */ | |
219 | if (eventT==END_OF_RUN) { | |
220 | printf("EOR event detected\n"); | |
221 | break; | |
222 | } | |
223 | } | |
224 | ||
225 | if (detTypeRec) delete detTypeRec; | |
226 | ||
227 | // Store the final histograms | |
228 | TFile outFile(OUTPUT_FILE, "update"); | |
229 | if (nevents_with_vertex > N_EVENTS_AUTOSAVE) { | |
230 | // Fit XY & Z histograms | |
231 | AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ); | |
232 | if (fitVtx) { | |
233 | fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite); | |
234 | TH1 *fithXY = fitFcn->CreateHistogram(); | |
235 | fithXY->Write(fithXY->GetName(),TObject::kOverwrite); | |
236 | delete fithXY; | |
237 | } | |
238 | delete fitVtx; | |
239 | } | |
240 | hXY->Write(hXY->GetName(),TObject::kOverwrite); | |
241 | hZ->Write(hZ->GetName(),TObject::kOverwrite); | |
242 | outFile.Close(); | |
243 | ||
55c5e86d | 244 | #ifdef ALI_AMORE |
245 | // send the histos to AMORE pool | |
246 | printf("AMORE send status: %d",vtxAmore.Send(hXY->GetName(),hXY)); | |
247 | #endif | |
248 | ||
eb35e591 | 249 | delete minuitFit; |
250 | TVirtualFitter::SetFitter(0); | |
251 | ||
252 | /* write report */ | |
253 | 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); | |
254 | ||
255 | status=0; | |
256 | ||
257 | /* export file to FXS */ | |
258 | if (daqDA_FES_storeFile(OUTPUT_FILE, "VertexDiamond")) { | |
259 | status=-2; | |
260 | } | |
261 | ||
262 | return status; | |
263 | } | |
264 | ||
265 | Double_t fitFunction(Double_t *x, Double_t *par) { | |
266 | ||
267 | Double_t t1 = x[0] - par[1]; | |
268 | Double_t t2 = x[1] - par[2]; | |
269 | ||
270 | return par[0]*TMath::Exp(-0.5*(t1*t1/(par[3]*par[3])+t2*t2/(par[4]*par[4]))); | |
271 | } | |
272 | ||
273 | AliESDVertex* FitVertexDiamond(TH2F *hXY, TH1F *hZ) | |
274 | { | |
275 | ||
276 | if (!fitFcn) { | |
277 | fitFcn = new TF2("fitFcn",fitFunction, | |
278 | -X_LIMIT,X_LIMIT, | |
279 | -Y_LIMIT,Y_LIMIT,NFITPARAMS); | |
280 | fitFcn->SetNpx(2*(Int_t)(X_LIMIT/X_DELTA)); | |
281 | fitFcn->SetNpy(2*(Int_t)(Y_LIMIT/Y_DELTA)); | |
282 | fitFcn->SetParameters(hXY->GetMaximum(),0,0,hXY->GetRMS(1),hXY->GetRMS(2)); | |
283 | fitFcn->Update(); | |
284 | } | |
285 | ||
286 | if (hXY->Fit("fitFcn","L0V+") != 0) { | |
287 | printf("XY fit failed!"); | |
288 | return 0x0; | |
289 | } | |
290 | ||
291 | Double_t pos[3],poserr[3]; | |
292 | pos[0] = fitFcn->GetParameter(1); | |
293 | pos[1] = fitFcn->GetParameter(2); | |
294 | poserr[0] = fitFcn->GetParameter(3); | |
295 | poserr[1] = fitFcn->GetParameter(4); | |
296 | ||
297 | // Could be replaced by something more robust... | |
298 | pos[2] = hZ->GetMean(); | |
299 | poserr[2] = hZ->GetRMS(); | |
300 | ||
301 | return new AliESDVertex(pos,poserr); | |
302 | } |