source: Sophya/trunk/Poubelle/archTOI.old/pisteetoile.cc@ 394

Last change on this file since 394 was 394, checked in by ansari, 26 years ago

Integration detecteur d'etoiles DY
File size: 10.4 KB
Line 
1#include "math.h"
2#include <new.h>
3#include "iostream.h"
4#include "archparam.h"
5#include "pisteetoile.h"
6
7TransFuncElec PisteEtoile::Trapani;
8FormePulse PisteEtoile::FPulse(Trapani);
9
10
11
12int SSTnbPasFit= (int) (2*TFExcursion/TPasTFit);
13
14PisteEtoile::PisteEtoile(short NoPhotDiod) {
15 PasTempsEch=0.005836818076;
16 //BUGG doit changer !!!!! Doit etre extrait d'un serveur de param d'acquisition
17 NoPhDiod=NoPhotDiod;
18 DiodVivante=true;
19 AmplNorm=0.;
20// On echantillonne dix fois plus fin pour avoir une bonne estimation de AmplNorm
21// Eviter les systematiques
22 for(int i=0; i<SSTLongIndexEvt()*10; i++){
23 AmplNorm+=FPulse.PulseShape(i/10.*TPasTFit)*FPulse.PulseShape(i/10.*TPasTFit);
24 }
25 AmplNorm/=10.;
26 indexDebutPiste=0;
27
28 PhDiodArray=new(nothrow) int[PhDiodTabLong];
29 PhDiodPiste=new(nothrow) double[PhDiodTabLong];
30 Ampl=new(nothrow) double[SSTnbPasFit];
31 X2=new(nothrow) double[SSTnbPasFit];
32
33 if((X2==0)||(Ampl==0)||(PhDiodPiste==0)||(PhDiodArray==0)) {
34 cerr<<"Erreur a la reservation de memoire dans constructeur pisteetoile"<<endl;
35 exit(-1);
36 }
37}
38
39PisteEtoile& PisteEtoile::operator =(PisteEtoile const & x) {
40 PasTempsEch=x.PasTempsEch;
41 //BUGG doit changer !!!!! Doit etre extrait d'un serveur de param d'acquisition
42 NoPhDiod=x.NoPhDiod;
43 DiodVivante=x.DiodVivante;
44 indexDebutPiste=x.indexDebutPiste;
45 DiodVivante=x.DiodVivante;
46 AmplNorm=x.AmplNorm;
47
48
49 for(int j=0; j<PhDiodTabLong;j++) PhDiodArray[j]=x.PhDiodArray[j];
50 for(int j=0; j<PhDiodTabLong;j++) PhDiodPiste[j]=x.PhDiodPiste[j];
51
52 LastEtoile=x.LastEtoile;
53
54 gachette=x.gachette;
55 Vrms=x.Vrms;
56 MaxIndex=x.MaxIndex; // Index de l'amplitude Max en analyse Trig()
57 AmplMax=x.AmplMax;
58
59 AmplRes=x.AmplRes;
60 TIndex=x.TIndex;
61 X2Min=x.X2Min;
62 SuccesAFine=x.SuccesAFine;
63 AmplNorm=x.AmplNorm;
64
65 return *this;
66}
67
68PisteEtoile::PisteEtoile (PisteEtoile const & x) {
69 PasTempsEch=x.PasTempsEch;
70 NoPhDiod=x.NoPhDiod;
71 DiodVivante=x.DiodVivante;
72 indexDebutPiste=x.indexDebutPiste;
73 DiodVivante=x.DiodVivante;
74 AmplNorm=x.AmplNorm;
75
76 PhDiodArray=new(nothrow) int[PhDiodTabLong];
77 PhDiodPiste=new(nothrow) double[PhDiodTabLong];
78 Ampl=new(nothrow) double[SSTnbPasFit];
79 X2=new(nothrow) double[SSTnbPasFit];
80
81
82 if((X2==0)||(Ampl==0)||(PhDiodPiste==0)||(PhDiodArray==0)) {
83 cerr<<"Erreur a la reservation de memoire dans constructeur pisteetoile"<<endl;
84 exit(-1);
85 }
86
87 for(int j=0; j<PhDiodTabLong;j++) PhDiodArray[j]=x.PhDiodArray[j];
88 for(int j=0; j<PhDiodTabLong;j++) PhDiodPiste[j]=x.PhDiodPiste[j];
89
90 LastEtoile=x.LastEtoile;
91
92 gachette=x.gachette;
93 Vrms=x.Vrms;
94 MaxIndex=x.MaxIndex; // Index de l'amplitude Max en analyse Trig()
95 AmplMax=x.AmplMax;
96
97 AmplRes=x.AmplRes;
98 TIndex=x.TIndex;
99 X2Min=x.X2Min;
100 SuccesAFine=x.SuccesAFine;
101 AmplNorm=x.AmplNorm;
102
103}
104
105
106
107PisteEtoile::~PisteEtoile() {
108 delete[] X2;
109 delete[] Ampl;
110 delete[] PhDiodPiste;
111 delete[] PhDiodArray;
112
113}
114
115void PisteEtoile::fill(int* pData, int InDebutPiste, int nData) {
116// if(!isAlive()) return;
117 indexDebutPiste=InDebutPiste;
118 int* PhDiod=PhDiodArray;
119 int* pDataC=pData;
120 for(int k=0;k<nData;k++,pDataC++, PhDiod++) (*PhDiod)=(*pDataC);
121 return;
122}
123
124void PisteEtoile::push(int* pData, int nbdata) {
125 int* PPhDiod=PhDiodArray;
126 int* PPhDiod2=PhDiodArray+nbdata;
127 int* pDataC=pData;
128
129 for(int i=0; i<PhDiodTabLong-nbdata; i++, PPhDiod++, PPhDiod2++) *PPhDiod=(*PPhDiod2);
130 for(int i=0; i<nbdata; i++, PPhDiod++,pDataC++) *PPhDiod=(*pDataC);
131 indexDebutPiste+=nbdata;
132 return;
133}
134
135bool PisteEtoile::traque() {
136 if(!isAlive()) return false; //La diode est-elle morte?
137
138 bool testEtoile=trig();
139 if(testEtoile) {
140 //On a declenché sur une Etoile
141 if(archParam.sst.analFine) {
142 demasque();
143// Un second Filtre est possible
144 // if(!(fabs(AmplRes)>SeuilGachette*Vrms)&&SuccesAFine) testEtoile=false;
145 // Test rajouté pour filtrer le bruit en escalier des voies SST
146 // Le fit de forme tasse le bruit et permet un meilleur filtrage.
147 }
148
149// Troisième filtre possible:
150 // On tue les evenements trop larges
151 if (widthHalfMax>(FPulse.FWHM()/PasTempsEch+2.)) testEtoile=false;
152
153 if(testEtoile) {
154 // On emballe une etoile.
155 LastEtoile.FitForme=SuccesAFine;
156 LastEtoile.X2Calc=X2Min;
157 LastEtoile.InpCurrNoisRMS=Vrms/FPulse.GainElectrique();
158 LastEtoile.TEchan=TIndex;
159 LastEtoile.InpCurrent= AmplRes/FPulse.GainElectrique();
160 LastEtoile.NoDiode=NoPhDiod;
161 LastEtoile.LargMiHauteur=widthHalfMax;
162 }
163 }
164 // On expedie
165 return testEtoile;
166}
167
168bool PisteEtoile::trig() {
169 MaxIndex=-1;
170 AmplMax=0.;
171 gachette=false;
172
173// Copie les données brutes sur la piste. Les met en Volt ADC 12 Bits-10 Volts
174 for(int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]=PhDiodArray[i]*SSTLsbVal;
175
176// Calcule Valeur Moy du prepulse
177 double Moy=0.;
178 for(int i=0; i<Prepulselong; i++) Moy+=PhDiodPiste[i];
179 Moy/=Prepulselong;
180
181// On soustrait la valeur Moy de la piste
182 for (int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]-=Moy;
183
184 double pente=0.;
185 if (archParam.sst.soustPente==true) {
186 // On calcule la pente du prepulse
187 double milprepulse=(Prepulselong-1.)/2.;
188 double SomCarre=0.;
189 for(int i=0; i<Prepulselong; i++) {
190 SomCarre+= (i-milprepulse)*(i-milprepulse);
191 pente+= (i-milprepulse)*PhDiodPiste[i];
192 }
193 pente=pente/SomCarre;
194
195 // On soustrait la pente sur la piste
196 for (int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]-=(i-milprepulse)*pente;
197 }
198
199// On calcule le Vrms sur le prepulse
200 Vrms=0.;
201 for (int i=0; i<Prepulselong; i++) Vrms+=PhDiodPiste[i]*PhDiodPiste[i];
202 Vrms=sqrt(Vrms)/Prepulselong;
203 if(Vrms<SSTLsbVal) Vrms=SSTLsbVal; // Pour eviter les VRMS=0.
204
205// Declenche-t-on?
206 gachette=false;
207 short iSeuil=0;
208 double* pPiste =PhDiodPiste+Prepulselong;
209
210// On demande deux échantillons à 2 Vrms
211 for (int i=0; i<(PhDiodTabLong-Prepulselong); i++, pPiste++) {
212 if (fabs(*pPiste)>(SeuilGachette*Vrms))
213 if (fabs(*(pPiste+1))>SeuilGachette*Vrms) {
214 gachette=true;
215 iSeuil=Prepulselong+i;
216 break;
217 }
218 }
219
220 if (gachette==true) {
221// On trouve le maximum
222 MaxIndex=iSeuil;
223 AmplMax=fabs(PhDiodPiste[iSeuil]);
224 for (int i=iSeuil; i<PhDiodTabLong; i++) {
225 if (fabs(PhDiodPiste[i])>AmplMax) {
226 MaxIndex=i;
227 AmplMax=fabs(PhDiodPiste[i]);
228 }
229 }
230 if ((MaxIndex>=PhDiodTabLong-4)||(MaxIndex<Prepulselong+4)) gachette=false;
231// Cela sera traite a l'iteration suivante
232// Ou ca a du etre traite a l'iteration precedente
233 }
234 if (gachette==true) {
235// On calcule la largeur a mi hauteur
236 int HalfMaxInf;
237 for(HalfMaxInf=MaxIndex;HalfMaxInf>Prepulselong;HalfMaxInf--)
238 if(fabs(PhDiodPiste[HalfMaxInf])<(AmplMax/2.)) break;
239 int HalfMaxSup;
240 for(HalfMaxSup=MaxIndex;HalfMaxSup<PhDiodTabLong;HalfMaxSup++)
241 if(fabs(PhDiodPiste[HalfMaxSup])<(AmplMax/2.)) break;
242 widthHalfMax=HalfMaxSup-HalfMaxInf;
243
244/*
245// On integre sur la larg a mi hauteur
246 integ=0.;
247 if((HalfMaxSup<PhDiodTabLong)&&(HalfMaxInf>Prepulselong))
248 for(int i=HalfMaxInf; i<=HalfMaxSup; i++) integ+=PhDiodPiste[i];
249*/
250
251
252// On fitte une parabole autour du maximum
253 double tmaxfit;
254 if (widthHalfMax>1) {
255 double ymax=PhDiodPiste[MaxIndex];
256 double yplus=PhDiodPiste[MaxIndex+1];
257 double ymoins=PhDiodPiste[MaxIndex-1];
258 // Estimation du temps
259 if(!(ymax==yplus)) {
260 double Y=(yplus-ymax)/(ymax-ymoins);
261 tmaxfit=0.5*(1.+Y)/(1-Y);
262 }
263 else tmaxfit=0.5;
264
265 // Estimation de l'amplitude
266 double dtplus2=(1.+tmaxfit)*(1.+tmaxfit);
267 double dtmoins2=(1.-tmaxfit)*(1.-tmaxfit);
268 if(dtplus2==dtmoins2) AmplRes=PhDiodPiste[MaxIndex];
269 else AmplRes= (yplus*dtplus2-ymoins*dtmoins2) / (dtplus2-dtmoins2);
270 }
271
272// On stocke les resultats
273 TIndex=indexDebutPiste+Prepulselong+MaxIndex+tmaxfit;
274
275 X2Min=0.;
276 SuccesAFine=false;
277 }
278 return gachette;
279}
280
281void PisteEtoile::demasque() {
282// Maintenant que l'on se doute que il y a une etoile
283// On veut connaitre son flux et son temps de passage
284// Nécessite d'avoir un ordre d'idee de la vitesse de rotation du ballon.
285
286 // BUGG devra etre extrait des headers acquisition. Ca peut changer!!!!!!
287
288
289
290 int IndexDebutEvt=MaxIndex-SSTLargEvtIndex();
291
292 pFitArray0= PhDiodPiste+IndexDebutEvt;
293 // Pointeur origine du tableau d'echantillon sur lequel sera calculé le fit
294 double DeltaTemps=0.;
295// Pour chaque decalage en temps, on estime un amplitude et calcule le X2
296 for(int k=0; k<SSTnbPasFit; k++) {
297 DeltaTemps=TFExcursion-k*TPasTFit;
298 pFitArray=pFitArray0;
299 double Som=0.;
300 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++) {
301// cout<<(*pFitArray)<<'\t'<<FPulse.PulseShape(i*PasTempsEch-DeltaTemps)<<endl;
302 Som=Som+FPulse.PulseShape(i*PasTempsEch-DeltaTemps)*(*pFitArray);
303 }
304 Ampl[k]=Som/AmplNorm;
305 X2[k]=0.;
306 double X2Som=0.; double dum=0.;pFitArray=pFitArray0;
307 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++) {
308 dum=(*pFitArray)-Ampl[k]*FPulse.PulseShape(i*PasTempsEch-DeltaTemps); // Différence
309 X2Som+=dum*dum; // Au carré
310 }
311 X2[k]=X2Som/(SSTLongIndexEvt()*Vrms*Vrms);
312 }
313
314// On recherche le minimum du X2
315 X2Min=1.e99;
316 int X2MinIndex=-1;
317 for(int k=0; k<SSTnbPasFit; k++)
318 if (X2[k]<X2Min) {
319 X2Min=X2[k];
320 X2MinIndex=k;
321 }
322
323 if((X2MinIndex>0)&&(X2MinIndex<(SSTnbPasFit-1))) {
324
325// On fitte une parabole sur les trois X2 qui encadrent le min
326 double Di; // Index (non entier) du minimum du X2
327 double X=(X2[X2MinIndex+1]-X2Min)/(X2Min-X2[X2MinIndex-1]);
328 Di=0.5*(X+1)/(X-1);
329
330// On calcule l'amplitude estimee a ce temps
331 double Som1=0.; pFitArray=pFitArray0;
332 DeltaTemps=TFExcursion-(X2MinIndex-Di)*TPasTFit;
333 for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++)
334 Som1+=(*pFitArray)*FPulse.PulseShape(i*PasTempsEch-DeltaTemps);
335 AmplRes=Som1/AmplNorm;
336
337//Le X2 est mimimum pour le n° de pas de fit X2MinIndex+Di
338
339#ifdef PEtoileDebug
340 double* PulseTrak=new double [SSTLongIndexEvt()];
341#endif
342
343 double X2Som1=0.; double dummy=0.;pFitArray=pFitArray0;
344 for(int i=0;i<SSTLongIndexEvt();i++,pFitArray++) {
345#ifdef PEtoileDebug
346 PulseTrak[i]=FPulse.PulseShape(i*PasTempsEch-DeltaTemps);
347#endif
348 dummy=(*pFitArray)-AmplRes*FPulse.PulseShape(i*PasTempsEch-DeltaTemps); // Différence
349 X2Som1+=dummy*dummy; // Au carré
350 }
351
352 X2Min=X2Som1/(SSTLongIndexEvt()*Vrms*Vrms);
353
354 TIndex=indexDebutPiste+IndexDebutEvt+DeltaTemps/PasTempsEch;
355 //Ouff!!!!
356 SuccesAFine=true;
357 }
358 else { // Echec, c'est tres bizarre! On doit etre dans le bruit
359 X2Min=0.;
360 }
361 return;
362}
363
364int PisteEtoile::SSTLargEvtIndex() {
365// Demi-largeur a 3 sigma de l'impulsion
366 int dum=(int) (15.e-3/PasTempsEch+1.);
367 return dum;
368}
369
370int PisteEtoile::SSTLongIndexEvt() {
371// Nombre d'echantillon sur lequel est calculé le fit de forme
372 int temp=2*SSTLargEvtIndex()+1;
373 return temp;
374}
375
Note: See TracBrowser for help on using the repository browser.