]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PMD/PMD_dig_reco.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Version 3 of the PMD. (Tapan K. Nayak)
[u/mrichter/AliRoot.git] / PMD / PMD_dig_reco.C
CommitLineData
9c02cab7 1#include "PMD_dig_reco.h"
2#include <iostream.h> // define cout stream
3#include <fstream.h> // defines ofstream class
4#include <stdlib.h> // defines exit() functions
5void PMD_dig_reco (Int_t evNumber=1)
6{
7/////////////////////////////////////////////////////////////////////////
8// This macro is a small example of a ROOT macro
9// illustrating how to read the output of GALICE
10// and fill some histograms.
11//
12//
13// Created by : Baba Potukuchi.V.K.S. ( JU / CERN ) 17/10/2000.
14// Modified by : Baba Potukuchi.V.K.S. ( JU / CERN ) 6/02/2001.
15// e-mail : potukuchi.baba@cern.ch
16//
17/////////////////////////////////////////////////////////////////////////
18// =========== Variables for Reconstruction==============================
19 FILE *fp4;
20 PMD_dig_reco* PPP = new PMD_dig_reco;
21// ======================================================================
22 Int_t padmult,cpvmult,padmult_gam,tot_trk,tot_trk_cpv;
23 Int_t padmult_cut, cpvmult_cut;
24 Int_t trk_Mult[100000];
25 Int_t trk_Mult_cpv[100000];
26 Float_t xPos,yPos,zPos,Edep,Vol1,Vol2,Vol2s,Vol3,Vol4,Vol5;
27 Float_t xpad,ypad;
28 Int_t iVol4,iVol5 , ixpad,iypad;
29 Int_t nbytes = 0;
30 // Float_t px,py,pz;
31 Float_t px1,py1,pz1;
32 Float_t mass;
33 Float_t edep_tot;
34 Int_t id_part;
35 Int_t gam_in;
36 Int_t jjj;
37 Int_t pip_in,pim_in;
38 TObjArray *points;
39 long double r,eta,px,py,pz;
40 Float_t phi,Eparticle,theta1,theta;
41 FILE *fp1;
42 FILE *fp2;
43
44 fp1=fopen("kine_vertex.out","w");
45 fp2=fopen("kine_pmdhit.out","w");
46 fp4=fopen("det1.out","w");
47
48// Dynamical link some shared libs
49 if (gClassTable->GetID("AliRun") < 0) {
50 gROOT->LoadMacro("loadlibs.C");
51 loadlibs();
52 }
53
54// Connect the Root Galice file containing Geometry, Kine and Hits
55 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
56 if (!file) file = new TFile("galice.root");
57
58// Get AliRun object from file or create it if not on file
59// if (!gAlice) {
60 gAlice = (AliRun*)file->Get("gAlice");
61 if (gAlice) printf("AliRun object found on file\n");
62 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
63// }
64 // Just for test
65//Loop over the requested number of events
66 for (Int_t iEvt = 0; iEvt < evNumber; iEvt++) {
67
68 printf("evt no %d \n",iEvt);
69 // End Just for test
70
71
72
73// Import the Kine and Hits Trees for the event evNumber in the file
74 Int_t nparticles = gAlice->GetEvent(iEvt);
75
76 printf ("N of particles %d\n",nparticles);
77
78 if (nparticles <= 0) return;
79
80 Int_t nevt = iEvt + 1;
81
82// ----********---- vertex track search BEGIN
83
84//Loop over the requested number of events
85// Import the Kine and Hits Trees for the event evNumber in the file
86
87 TObjArray *Particles = gAlice->Particles();
88 Float_t phi;
89 Float_t pl,cutmax,etamin,etamax,pmom,smin,smax;
90 Float_t *pxyz;
91 TObjArray *points;
92 TParticle *particle;
93 AliPoints *pm;
94 // for (Int_t it = 0; it < 100000; it++)trk_Mult[it]=0.;
95
96 for (Int_t in_track = 0; in_track < nparticles; in_track++) {
97 // for (Int_t in_track = 0; in_track < 10; in_track++) {
98 particle = (TParticle*)gAlice->Particle(in_track);
99 if (!particle) continue;
100 id_part = particle->GetPdgCode();
101 // mass = TDatabasePDG::Instance()->GetParticle(id_part)->Mass();
102 // printf("_____________--------%d %f\n",id_part,mass);
103
104 px = particle->Px();
105 py = particle->Py();
106 pz = particle->Pz();
107 if(px !=0 || py!=0 || pz!=0) {
108 if( (px!=0 || py!=0) && pz < 0 ){
109// if(id_part == 22)gam_in++;
110 // printf("%d\n",gam_in);
111 r = TMath::Sqrt(px*px + py*py);
112 theta = TMath::ATan2(r,TMath::Abs(pz));
113 phi = TMath::ATan2(py,px);
114 phi = (phi*180.)/3.14;
115// if(theta < -1.5707961)printf( " THETA=%f\n",theta) ;
116 if(theta > -1.5707961)eta = -TMath::Log(TMath::Tan(0.5*theta));
117 theta1=theta*180./3.14159;
118
119// if(TMath::Abs(eta) > 2.3 && TMath::Abs(eta) < 3.5){
120 if(id_part == 22)gam_in++;
121 if(id_part == 211)pip_in++;
122 if(id_part == -211)pim_in++;
123 fprintf(fp1," %d %d %f %f %f %f %f %f \n",nevt,id_part,px,py,pz,eta,theta1,phi);
124// }
125// printf("%d %d %d\n",gam_in,pip_in,pim_in);
126// printf("%f %f %f %f %d\n",ptg,pxg,pyg,pzg,id_part);
127 }
128 }
129 }
130
131// 1999 file-- Get pointers to Alice detectors and Hits containers
132 //------********----- vertex track ENDs
133
134 Float_t x,y,z,mass,e;
135 Int_t nbytes = 0;
136 Int_t j,hit,ipart;
137 Int_t npmd;
138 Int_t mbytes = 0;
139 Int_t nhits,gam_in;
140 Int_t ipmd,id_part;
141 Int_t ypad1;
142 Int_t mpmd;
143 const Int_t Npad2 = 72;
144
145 Int_t dVol;
146 // Int_t iEvt;
147
148 Int_t sector,plane;
149 TParticle *particle;
150 AliPMDhit *pmdHit;
151
152// Get pointers to Alice detectors and Hits containers
153 AliDetector *PMD = gAlice->GetDetector("PMD");
154// TClonesArray *fParticles = gAlice->Particles();
155 if (PMD) TClonesArray *PMDhits = PMD->Hits();
156 if (PMD) TObjArray *points = PMD->Points();
157
158 TTree *TH = gAlice->TreeH();
159 // TTree *TK = gAlice->TreeK();
160
161 Int_t ntracks = TH->GetEntries();
162 printf("No.of Tracksin the TreeH %d \n", ntracks);
163
164// =========== Initializing the array ============================
165 for (Int_t iDet = 1; iDet <= 2 ; ++iDet) {
166 printf(" Initializing for Detector %d \n",iDet) ;
167 for (Int_t iBox = 1; iBox <= 30 ; ++iBox) {
168 for (Int_t iCol = 1; iCol <= 72 ; ++iCol) {
169 for (Int_t iRow = 1; iRow <= 72 ; ++iRow) {
170 if(iDet==1)(*PPP).module1 [iBox] [iCol] [iRow] = 0.0 ;
171 if(iDet==2)(*PPP).module2 [iBox] [iCol] [iRow] = 0.0 ;
172 }
173 }
174 }
175 }
176// Creating histograms .
177 TH1F *edp = new TH1F("edp","Energy Deposited in each cell" ,100,100.,100000.);
178 TH1F *Modnum = new TH1F("Modnum","Module number" ,100,1.,100.);
179 TH2F *xycells = new TH2F("xycells","XY-of Cells in an Unit Modlule " ,100,0.,80.,100,0.,80.);
180 TH2F *xyhits = new TH2F("xyhits","Hits of the cells " ,100,-150.,150.,100,-150.,150.);
181
182// ================================================================
183
184// =======================================================
185// Start loop on tracks in the hits containers
186 for (Int_t track=0; track<ntracks;track++) {
187 gAlice->ResetHits();
188
189 //changes introduced ******
190
191 nbytes += TH->GetEvent(track);
192 // mbytes += TK->GetEvent(track);
193 // printf("%d ******** %d \n" ,track,nbytes);
194 // ======>Histogram PMD
195
196 if (PMD) {
197 npmd = PMDhits->GetEntriesFast();
198 //printf("npmd=%d, for track=%d %d \n" ,npmd,track,nbytes);
199 for (ipmd = 0; ipmd < npmd; ipmd++) {
200 pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
201 ipart = pmdHit->GetTrack();
202
203 // printf("ID part and hits on PMD --------%d %d \n" ,ipart,track);
204// get kinematics of the particles
205// particle = (TParticle*)fParticles->UncheckedAt(ipart);
206 particle = gAlice->Particle(ipart);
207// if (!particle) continue;
208 //
209 id_part = particle->GetPdgCode();
210 mass = TDatabasePDG::Instance()->GetParticle(id_part)->Mass();
211 // printf("ID part and hit on PMD %d %d \n" ,id_part,ipmd);
212 px = particle->Px();
213 py = particle->Py();
214 pz = particle->Pz();
215 Eparticle = particle->Energy();
216
217 r = TMath::Sqrt(px*px + py*py);
218 theta = TMath::ATan2(r,TMath::Abs(pz));
219 //phi = TMath::ATan2(TMath::Abs(py),TMath::Abs(px));
220 phi = TMath::ATan2(py,px);
221 phi = (phi*180.)/3.14159;
222 eta = -TMath::Log(TMath::Tan(0.5*theta));
223 if (pz < 0) eta = -eta;
224 theta=theta*180./3.14159;
225
226// if(TMath::Abs(eta) > 2.3 && TMath::Abs(eta) < 3.5){
227 //
228 xPos = pmdHit->X();
229 yPos = pmdHit->Y();
230 zPos = pmdHit->Z();
231 Vol1 = pmdHit->fVolume[0];
232 Vol2 = pmdHit->fVolume[1];
233 Vol3 = pmdHit->fVolume[2];
234 Vol4 = pmdHit->fVolume[3];
235 Vol5 = pmdHit->fVolume[4];
236
237 Vol2s = Vol2-1;
238 ypad = Vol2s/Npad2 + 1;
239 ypad1 = ypad-1;
240 xpad = Vol2-ypad1*Npad2;
241 Edep = pmdHit->fEnergy;
242 ixpad = xpad+0.2 ;
243 iypad = ypad+0.2 ;
244 iVol4 = Vol4+0.2 ;
245 iVol5 = Vol5+0.2 ;
246// module [Vol4] [Vol5] [xpad] [ypad] = Edep ;
247 if(Edep>0)printf(" VOL4 %d %d %d %d %f \n",iVol4,iVol5,ixpad,iypad,Edep) ;
248 edp ->Fill(Edep) ;
249 Modnum ->Fill(iVol5) ;
250 xycells->Fill(ixpad,iypad);
251 xyhits->Fill(xPos,yPos);
252
253 if(iVol4==1)(*PPP).module1 [iVol5] [ixpad] [iypad] = Edep ;
254 if(iVol4==2)(*PPP).module2 [iVol5] [ixpad] [iypad] = Edep ;
255 Float_t XXP = (*PPP).module2 [iVol5] [ixpad] [iypad] ;
256 if(Edep>0)printf(" VOL4 %d %d %d %d %f \n",Vol4,Vol5,xpad,ypad,XXP) ;
257 // Write to the kine.out file
258 fprintf(fp2," %d %d %d %f %f %f %d %d %d %d %f %f %f %f %f \n",
259 nevt,ipart,id_part,xPos,yPos,zPos,xpad,ypad,Vol5,Vol4,Edep,px,py,pz,Eparticle);
260
261// Closing Eta loop of the incident particles.
262// }
263
264 }
265 }
266 }
267 }
268
269 pmdreco_(PPP);
270
271// TCanvas *c1 = new TCanvas("c1","PMD Simulated digits(Edep) ",400,10,600,700);
272// c1->Divide(1,2);
273// c1->cd(1);
274// gPad->SetFillColor(33);
275// edp->SetFillColor(42);
276// edp->Draw();
277// c1->cd(2);
278// gPad->SetFillColor(33);
279// Modnum->SetFillColor(41);
280// Modnum->Draw();
281// c1->Print("baba.ps");
282// TCanvas *c2 = new TCanvas("c2","PMD Simulated digits(Edep) ");
283// gPad->SetFillColor(33);
284// xycells->SetFillColor(42);
285// xycells->Draw();
286// c2->Print("baba1.ps");
287// TCanvas *c3 = new TCanvas("c3","PMD Simulated digits(Edep) ");
288// gPad->SetFillColor(33);
289// xyhits->SetFillColor(42);
290// xyhits->Draw();
291// c3->Print("baba2.ps");
292}
293
294// =========== RECONSTRUCTION ======================================
295void pmdreco_(PMD_dig_reco* PPP) {
296 Float_t sqrth;
297 sqrth = sqrt(3.) / 2.;
298
299 printf(" I am in MAIN calling RanMAR & Intpts_ \n") ;
300 rmarin_(PPP);
301 intpts_(PPP);
302// Generate COORDINATES .
303 for (Int_t i__ = 1; i__ <= 72; ++i__) {
304 for (Int_t j = 1; j <= 72; ++j) {
305 (*PPP).coord [1] [i__] [j] = i__ + j / 2.;
306 (*PPP).coord [2] [i__] [j] = j * sqrth;
307 }
308 }
309//-------------------------
310 for (Int_t iDet = 1; iDet <= 2 ; ++iDet) {
311 printf("idet= %d \n",iDet) ;
312 for (Int_t iBox = 1; iBox <= 30 ; ++iBox) {
313 for (Int_t iCol = 1; iCol <= 72 ; ++iCol) {
314 for (Int_t iRow = 1; iRow <= 72 ; ++iRow) {
315 Float_t PVAL = (*PPP).module1 [iBox] [iCol] [iRow];
316// printf("Det1 = %f %f %f %f %f \n",iDet, iBox, iCol, iRow,PVAL ) ;
317 if(iDet==1)(*PPP).det1 [iCol] [iRow] = (*PPP).module1 [iBox] [iCol] [iRow] ;
318 if(iDet==2)(*PPP).det2 [iCol] [iRow] = (*PPP).module2 [iBox] [iCol] [iRow] ;
319 if(PVAL>0.0)printf("Det1 = %d %d %d %d %f \n",iDet, iBox, iCol, iRow,PVAL ) ;
320// fprintf(fp4," Det1 = %d %d %d %d %f \n",iDet, iBox, iCol, iRow,VAL ) ;
321
322 }
323 }
324// Ending Reading pads in the Box .
325 }
326// Ending Reading pads for 1-CPV,2-PMD.
327 }
328/* analyse detector 1 */
329/* CPV clustering. Ordering data */
330 printf(" I am in MAIN Doing Clustering for CPV. \n") ;
331 printf (" I am in MAIN calling ==order== for CPV(det=1) \n");
332 order_(1,PPP);
333 printf(" I am in MAIN after order for CPV(det=1) \n") ;
334 printf(" I am in MAIN calling ==crclust== \n") ;
335 crclust_(1, PPP) ;
336 printf(" I am in MAIN after ==crclust== & now call prints \n") ;
337 prints_(1,PPP) ;
338 printf(" I am in MAIN Doing Clustering for PMD. \n") ;
339printf (" I am in MAIN calling ==order== for PMD(det=2) \n");
340 order_(2,PPP);
341 printf(" I am in MAIN after order for PMD(det=2) \n") ;
342 printf(" I am in MAIN calling ==crclust== \n") ;
343 crclust_(2, PPP) ;
344 printf(" I am in MAIN after ==crclust== & now call prints \n") ;
345 prints_(2,PPP) ;
346 (*PPP).incr = 0;
347 for (Int_t i__ = 1; i__ <= 72; ++i__) {
348 for (j = 1; j <= 72; ++j) {
349 if ((*PPP).infocl[1] [i__] [j] == 1) {
350 ++(*PPP).incr;
351 (*PPP).inford[1] [(*PPP).incr] = (*PPP).infocl[2] [i__] [j] ;
352 (*PPP).inford[2] [(*PPP).incr] = i__;
353 (*PPP).inford[3] [(*PPP).incr] = j;
354 }
355 }
356 }
357 for (i__ = 1; i__ <= 72; ++i__) {
358 for (j = 1; j <= 72; ++j) {
359 if ((*PPP).infocl[1] [i__] [j] == 2) {
360 ++(*PPP).incr;
361 (*PPP).inford[1] [(*PPP).incr] = (*PPP).infocl[2] [i__] [j];
362 (*PPP).inford[2] [(*PPP).incr] = i__;
363 (*PPP).inford[3] [(*PPP).incr] = j;
364 }
365 }
366 }
367 arrange_(PPP);
368/* better clustering ( Gaussian fit ) */
369 refclust_(PPP);
370 return ;
371} /* LOOP pmdreco */
372//
373 /* >>>>>>>>>>>> */void rmarin_(PMD_dig_reco* PPP)
374{
375 Int_t i__, j, k, l, m;
376 Float_t s, t;
377/* Initializing routine for RANMAR, must be called before */
378/* generating any pseudorandom numbers with RANMAR. The input */
379/* values should be in the ranges: 0<=IJ<=31328 */
380/* 0<=KL<=30081 */
381/* To get the standard values in Marsaglia's paper, IJ=1802, KL=9373 */
382 Int_t ij = 8775;
383 Int_t kl = 5926;
384 if (ij <= 0 || ij >= 31328) {
385 goto L900;
386 }
387 if (kl <= 0 || kl >= 30081) {
388 goto L900;
389 }
390 printf("I am in RanMAR \n");
391 i__ = ij / 177 % 177 + 2;
392 j = ij % 177 + 2;
393 k = kl / 169 % 178 + 1;
394 l = kl % 169;
395 printf(" ranmar initialized: %d %d %d %d %d %d \n",ij,kl,i__,j,k,l);
396 for (Int_t ii = 1; ii <= 97; ++ii) {
397 s = 0.;
398 t = .5;
399 for (Int_t jj = 1; jj <= 24; ++jj) {
400 m = (((i__ * j) % 179) * k) % 179;
401 i__ = j;
402 j = k;
403 k = m ;
404 l = l * 53 + 1 % 169;
405 if (l * m % 64 >= 32) {
406 s += t;
407 }
408 t *= .5;
409 }
410 (*PPP).u[ii - 1] = s;
411 }
412 (*PPP).c__ = .021602869033813477;
413 (*PPP).cd = .45623308420181274;
414 (*PPP).cm = .99999982118606567;
415 printf(" Exiting from rmarin %f %f %f \n",(*PPP).c__,(*PPP).cd,(*PPP).cm) ;
416 return ;
417L900:
418 printf( " +++++++ stop in rmarin: wrong params! ");
419} /* rmarin_ */
420
421/*>>>>>>>>>>>>>>>> */void order_(Int_t jdet,PMD_dig_reco* PPP)
422{
423 Int_t nmx,i__1, i__2, i__3;
424
425 Float_t adum;
426 Int_t idum, i__, j, k, l;
427 Float_t d1[5185];
428 Float_t iord1[5185];
429 Int_t itest, i1, i2, id ;
430
431 printf("I am in order for DET=%d \n",jdet);
432 nmx = 5184;
433 for (i1 = 1; i1 <= 72; ++i1) {
434 for (i2 = 1; i2 <= 72; ++i2) {
435 i__ = (i2 - 1) * 72 + i1;
436 d1[i__] = 0.0 ;
437 iord1 [i__] = 0;
438 if(jdet==1) d1[i__] = (*PPP).det1 [i1] [i2] ;
439 if(jdet==2) d1[i__] = (*PPP).det2 [i1] [i2] ;
440 iord1 [i__] = i__;
441 }
442 }
443 for (j = 2; j <= nmx; ++j) {
444 itest = 0;
445 adum = d1[j];
446 idum = iord1 [j];
447 i__2 = j - 1;
448 for (k = 1; k <= i__2; ++k) {
449 if (adum > d1[k ] && itest == 0) {
450 itest = 1 ;
451 i__3 = k;
452 for (l = j - 1; l >= i__3; --l) {
453 d1[l+1] = d1[l];
454 iord1 [l+1] = iord1 [l];
455 }
456 d1[k] = adum;
457 iord1 [k] = idum;
458 }
459 }
460 }
461 i__1 = nmx;
462 for (i__ = 1; i__ <= i__1; ++i__) {
463 id = iord1 [i__];
464 i2 = id / 72 + 1;
465 i1 = id - (i2 - 1) * 72;
466 (*PPP).iorder[1] [i__] = i1;
467 (*PPP).iorder[2] [i__] = i2;
468
469 printf("I am in order for DET=%d %d %d %d %d %d \n",jdet,i__,i1,i2,(*PPP).iorder[1] [i__],(*PPP).iorder[2] [i__]);
470 }
471} /* order_ */
472
473/* crclust_ */ void crclust_(Int_t ndet, PMD_dig_reco* PPP)
474 {
475 Float_t d__ [73] [73] ;
476
477 Int_t neibx[7] = { 0,1,0,-1,-1,0,1 };
478 Int_t neiby[7] = { 0,0,1,1,0,-1,-1 };
479
480 Int_t i__1;
481 Float_t sqrth ;
482
483 Int_t icld, idum, i__, j, k;
484 Int_t nb, id1, id2, jd1, jd2, icle, icn, itr, nmx ;
485
486 printf(" neibx %d %d %d %d %d %d \n",neibx[1],neibx[2],neibx[3],neibx[4],neibx[5],neibx[6]) ;
487 printf(" neiby %d %d %d %d %d %d \n",neiby[1],neiby[2],neiby[3],neiby[4],neiby[5],neiby[6]) ;
488
489 for (j = 1; j <= 72; ++j) {
490 for (k = 1; k <= 72; ++k) {
491 d__ [j] [k] = 0.0 ;
492 if(ndet==1) d__ [j] [k] = (*PPP).det1 [j] [k] ;
493 if(ndet==2) d__ [j] [k] = (*PPP).det2 [j] [k] ;
494
495 }
496 }
497 nmx = 5148;
498 sqrth=TMath::Sqrt(3.)/2. ;
499/* initialise infocl */
500 for (i__ = 1; i__ <= 2; ++i__) {
501 for (j = 1; j <= 72; ++j) {
502 for (k = 1; k <= 72; ++k) {
503 (*PPP).infocl[i__] [j] [k] = 0;
504 }
505 }
506 }
507/* compute number of pads with nonzero count ( nmx1 ) */
508 i__1 = nmx;
509 for (i__ = 1; i__ <= i__1; ++i__) {
510 id1 = (*PPP).iorder [1] [i__];
511 id2 = (*PPP).iorder [2] [i__];
512 if (d__[id1] [id2] > 0.0000000000001) {
513 (*PPP).nmx1 = i__;
514 printf(" nmx1 -> %d %d %f %d \n",id1,id2, d__ [id1] [id2],(*PPP).nmx1 ) ;
515 }
516 }
517 printf(" nmx1= %d \n",(*PPP).nmx1) ;
518/* crude clustering. start with pad with largest count ( cluster 1 )
519*/
520 icle = 1;
521 idum = 0;
522 id1 = (*PPP).iorder [1] [1];
523 id2 = (*PPP).iorder [2] [1];
524 (*PPP).infocl[1] [id1] [id2] = 1;
525 (*PPP).infocl[2] [id1] [id2] = icle;
526 printf(" infocl %d %d %d %d %d \n",id1,id2,icle,(*PPP).infocl[1] [id1] [id2],(*PPP).infocl[2] [id1] [id2] ) ;
527/* check if neighbouring count in nonzero. If yes, include in the */
528/* cluster */
529/* infocl(1,i1,i2) =1 for primary ( peak ) =2 for secondary ( */
530/* neighbour ) */
531/* infocl(2,i1,i2) stores cluster no. */
532 for (nb = 1; nb <= 6; ++nb) {
533 jd1 = id1 + neibx[nb];
534 jd2 = id2 + neiby[nb];
535 printf(" infocl %d %d %d %d %d %f %d %d\n",nb,jd1,id1,jd2,id2,d__[jd1] [jd2],icle,(*PPP).infocl[1] [jd1] [jd2] ) ;
536 if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
537 if (d__[jd1] [jd2] > 0.0) {
538 (*PPP).infocl[1] [jd1] [jd2] = 2;
539 (*PPP).infocl[2] [jd1] [jd2] = icle;
540 }
541 }
542 printf(" infocl %d %d %d %d %d %f %d %d\n",nb,jd1,id1,jd2,id2,d__[jd1] [jd2],icle,(*PPP).infocl[1] [jd1] [jd2] ) ;
543 }
544/* start checking with next largest count and so on */
545 i__1 = (*PPP).nmx1;
546 for (i__ = 2; i__ <= i__1; ++i__) {
547 itr = 0;
548 id1 = (*PPP).iorder[1] [i__];
549 id2 = (*PPP).iorder[2] [i__];
550 printf("iorder %d %d %d %d %d %d %d\n",i__,id1,id2,(*PPP).infocl [1] [id1] [id2],jd1,jd2,(*PPP).infocl [1] [jd1] [jd2] ) ;
551/* nonzero infocl(1,id1,id2) means the pad is included in earlier
552*/
553/* cluster ( itr = 1 ) Otherwise itr = 0. */
554 if ((*PPP).infocl [1] [id1] [id2] > 0) {
555 itr = 1;
556 }
557 if (itr == 0) {
558 for (nb = 1; nb <= 6; ++nb) {
559 jd1 = id1 + neibx[nb ];
560 jd2 = id2 + neiby[nb ];
561/* check if a neighbouring pad belongs to previous cluster
562. If yes, */
563/* icn = 1. Otherwise icn = 0. */
564 icn = 0;
565 if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
566 if ((*PPP).infocl [1] [jd1] [jd2] != 0) {
567 icn = 1;
568 icld = (*PPP).infocl [2] [jd1] [jd2] ;
569 }
570 }
571 }
572/* itr = 0 and icn = 0 means new cluster. Store the cluster in
573fo. in */
574/* infocl. */
575 if (TMath::Abs(icn) <= 0.1) {
576 ++icle;
577 (*PPP).infocl[1] [id1] [id2] = 1;
578 (*PPP).infocl[2] [id1] [id2] = icle;
579 printf("ID1 icn %d %d %d %d %d\n",icn,icle,id1,id2,(*PPP).infocl[2] [id1] [id2]);
580 for (nb = 1; nb <= 6; ++nb) {
581/* neighbours of the peak put with the peak. */
582 jd1 = id1 + neibx[nb ];
583 jd2 = id2 + neiby[nb ];
584 if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
585 if (d__[jd1] [jd2] != 0 && (*PPP).infocl [1] [jd1] [jd2] == 0) {
586 (*PPP).infocl[1] [jd1] [jd2] = 2;
587 (*PPP).infocl[2] [jd1] [jd2] = icle;
588 printf("JD1 icn %d %d %d %d %d\n",icn,icle,jd1,jd2,(*PPP).infocl[2] [jd1] [jd2]);
589 }
590 }
591 }
592 } else {
593/* itr = 0 and icn = 1 means the pad belongs to cluster ic
594ld */
595 (*PPP).infocl[1] [id1] [id2] = 2;
596 (*PPP).infocl[2] [id1] [id2] = icld;
597 for (nb = 1; nb <= 6; ++nb) {
598/* and also the neighbours. */
599 jd1 = id1 + neibx[nb ];
600 jd2 = id2 + neiby[nb ];
601 if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
602 if (d__[jd1] [jd2] != 0. && (*PPP).infocl[1] [jd1] [jd2] == 0) {
603 (*PPP).infocl[1] [jd1] [jd2] = 2;
604 (*PPP).infocl[2] [jd1] [jd2] = icld;
605 }
606 }
607 }
608 } else {
609/* itr = 0, so not a new cluster */
610 icld = (*PPP).infocl[2] [id1] [id2] ;
611 for (nb = 1; nb <= 6; ++nb) {
612/* neighbours put with the previous cluster. */
613 jd1 = id1 + neibx[nb ];
614 jd2 = id2 + neiby[nb ];
615 if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
616 if (d__[jd1] [jd2] != 0. && (*PPP).infocl[1] [jd1] [jd2] == 0) {
617 (*PPP).infocl[1] [jd1] [jd2] = 2;
618 (*PPP).infocl[2] [jd1] [jd2] = icld;
619 }
620 }
621 }
622 }
623 }
624}
625 printf(" ndet=%d,nmx1= %d , icle= %d icld= %d \n",ndet,(*PPP).nmx1,icle,icld) ;
626
627/* L2: */
628/* L3: */
629 return ;
630} /* crclust_ */
631
632/* Subroutine */ void intpts_(PMD_dig_reco* PPP)
633{
634 Int_t i__1, i__2, i__3;
635
636 Int_t ixpt;
637 Float_t x, y, dx, dy;
638 Int_t ix, iy;
639 Float_t yl;
640 Int_t idx, idy, ipt;
641 Float_t wtx, wty;
642
643 ixpt = 2;
644 idx = -1;
645 ipt = 0;
646 dx = .25;
647 i__1 = ixpt;
648 for (ix = -ixpt; ix <= i__1; ++ix) {
649 x = ix * dx;
650 idx = -idx;
651 i__2 = ix / ixpt ;
652 wtx = (3. - idx - (TMath::Abs(i__2) ))/ 3. * dx;
653 yl = 1. / TMath::Sqrt(3.) * (1. - TMath::Abs(x));
654 dy = yl * 2. / 4.;
655 idy = -1;
656 i__2 = ixpt;
657 for (iy = -ixpt; iy <= i__2; ++iy) {
658 idy = -idy;
659 y = iy * dy;
660 i__3 = iy / ixpt ;
661 wty = (3. - idy - (TMath::Abs(i__3))) / 3. * dy;
662 ++ipt;
663 (*PPP).pts[1] [ipt] = x;
664 (*PPP).pts[2] [ipt] = y;
665 (*PPP).wts[ipt] = wtx * wty;
666 }
667 }
668 printf(" Exiting from intpts_ %f %f \n", ipt,(*PPP).wts[ipt] );
669 return ;
670} /* intpts_ */
671
672/* prints_ */ void prints_(Int_t ndet, PMD_dig_reco* PPP)
673{
674 Int_t i__, j;
675
676 if (ndet == 1) {
677 FILE *fp1;
678 FILE *fp2;
679 fp1=fopen("CPV_prime.out","w");
680 fp2=fopen("CPV_secondary.out","w");
681 }
682 if (ndet == 2) {
683 FILE *fp1;
684 FILE *fp2;
685 fp1=fopen("PMD_prime.out","w");
686 fp2=fopen("PMD_secondary.out","w");
687 }
688 for (i__ = 1; i__ <= 72; ++i__) {
689 for (j = 1; j <= 72; ++j) {
690 if( (*PPP).infocl [1] [i__] [j] == 1) {
691 fprintf(fp1," %f %f %d %d %d \n",(*PPP).coord [1] [i__] [j],(*PPP).coord [2] [i__] [j],i__,j,(*PPP).infocl [2] [i__] [j] ) ;
692 }
693 if( (*PPP).infocl [1] [i__] [j] == 2) {
694 fprintf(fp2," %f %f %d %d %d \n",(*PPP).coord [1] [i__] [j],(*PPP).coord [2] [i__] [j],i__,j,(*PPP).infocl [2] [i__] [j] ) ;
695 }
696 }
697 }
698 return ;
699} /* prints_ */
700
701/* Subroutine */void arrange_(PMD_dig_reco* PPP)
702{
703 Int_t i__1, i__2, i__3;
704
705 Int_t ihld1, ihld2, ihld3, i__, j, k, itest;
706
707 i__1 = (*PPP).incr;
708 for (i__ = 2; i__ <= i__1; ++i__) {
709 itest = 0;
710 ihld1 = (*PPP).inford [1] [i__];
711 ihld2 = (*PPP).inford [2] [i__];
712 ihld3 = (*PPP).inford [3] [i__];
713 i__2 = i__ - 1;
714 for (j = 1; j <= i__2; ++j) {
715 if (itest == 0 && (*PPP).inford [1] [j] > ihld1) {
716 itest = 1;
717 i__3 = j;
718 for (k = i__ -1; k >= i__3; --k) {
719 (*PPP).inford [1] [k+1] = (*PPP).inford [1] [k] ;
720 (*PPP).inford [2] [k+1] = (*PPP).inford [2] [k] ;
721 (*PPP).inford [3] [k+1] = (*PPP).inford [3] [k] ;
722 }
723 (*PPP).inford [1] [j] = ihld1 ;
724 (*PPP).inford [2] [j] = ihld2 ;
725 (*PPP).inford [3] [j] = ihld3 ;
726 }
727 }
728 }
729 printf(" Exiting from arrange \n") ;
730 return ;
731} /* arrange_ */
732
733/* refclust_ */void refclust_(PMD_dig_reco* PPP)
734{
735
736 Int_t i__1, i__2;
737 Float_t d__ [73] [73] ;
738 Int_t npad[201], i__, j, k;
739 Int_t id, id1, id2, iclb;
740
741/* using information from crude clustering, a refined algorithm for */
742/* clustering tried. At the moment using Gaussian fit with a number */
743/* of Gaussian functions. */
744 for (j = 1; j <= 72; ++j) {
745 for (k = 1; k <= 72; ++k) {
746 d__ [j] [k] = (*PPP).det2 [j] [k] ;
747 }
748 }
749
750 (*PPP).clno = 0;
751/* The No.of pads in each cluster = 200 */
752 for (i__ = 1; i__ <= 200; ++i__) {
753 npad[i__] = 0;
754 }
755 for (i__ = 1; i__ <= 5000; ++i__) {
756 (*PPP).ncl[i__] = 0;
757 }
758/* number of pads in each cluster determined. */
759 i__1 = (*PPP).incr;
760 for (i__ = 1; i__ <= i__1; ++i__) {
761 id = (*PPP).inford[1] [i__];
762 ++(*PPP).ncl[id];
763 Int_t nclid = (*PPP).ncl[id] ;
764 printf("id,ncl= %d %d\n",id,(*PPP).ncl[id]) ;
765 }
766 iclb = 0;
767/* number of clusters determined. */
768 for (i__ = 1; i__ <= 5000; ++i__) {
769 id = (*PPP).ncl[i__];
770 if (id != 0 ) {
771 ++iclb;
772 ++npad[id];
773 }
774 }
775 printf("iclb= %d \n",iclb) ;
776 id = 0;
777/* define coordinates of each pad in a cluster and call gausspara */
778/* to determine Gaussian parameters. ( should it be eliminated ? ) */
779 i__1 = iclb;
780 for (i__ = 1; i__ <= i__1; ++i__) {
781 i__2 = (*PPP).ncl[i__];
782 printf("(*PPP).ncl(i)= %d %d\n",i__,(*PPP).ncl[i__]) ;
783 for (j = 1; j <= i__2; ++j) {
784 ++id;
785 id1 = (*PPP).inford[2] [id];
786 id2 = (*PPP).inford[3] [id];
787 (*PPP).x[j] = (*PPP).coord[1] [id1] [id2];
788 (*PPP).y[j] = (*PPP).coord[2] [id1] [id2];
789 (*PPP).z__[j] = d__[id1] [id2];
790 }
791 printf(" I am in refclust calling ==gausspara_== \n") ;
792 gausspara_(i__2,PPP );
793 printf(" I am in refclust after ==gausspara_== \n") ;
794 }
795 printf("clno= %d \n",(*PPP).clno);
796 return ;
797} /* refclust_ */
798
799/* gausspara_ */ void gausspara_(Int_t n,PMD_dig_reco* PPP)
800{
801 Int_t i__1, i__2, i__3;
802 Float_t d__1, d__2;
803
804 Int_t ihld, iord[501], idum, i__, j, k;
805 Int_t itest;
806 Int_t ig;
807 Float_t pi, xi, yi, xj, yj, rrr;
808 Float_t sum;
809
810/* single Gaussian of the form exp(-((x-x_0)^2+(y-y_0)^2)/r_0^2) */
811
812 printf(" I am in ==gausspara_== %d \n" , n) ;
813 pi = 3.141593;
814 Int_t nn = n ;
815 if (nn == 1) {
816 ++(*PPP).clno;
817/* single pad cluster. Assumed fit to single Gaussian centered at the */
818/* center of the pad with Gaussian parameter r0 about 0.5 since */
819/* neighbouring pads don't fire. */
820 (*PPP).x0[1] = (*PPP).x[1];
821 (*PPP).y0[1] = (*PPP).y[1];
822 (*PPP).r0[1] = 0.5;
823/* Computing 2nd power */
824 d__1 = (*PPP).r0[1];
825 (*PPP).h0[1] = (*PPP).z__[1] / pi / (d__1 * d__1);
826 printf("single pad ,clust no. = %d %d %f %f \n",n, (*PPP).clno, (*PPP).x0[1],(*PPP).y0[1] ) ;
827 } else if (nn == 2) {
828 ++((*PPP).clno);
829/* two pad cluster. Assumed fit to single Gaussian centered at the */
830/* weighted average of the centers of two pads with r_0 about 0.5. */
831/* Could be a pair of clusters but no way of differentiating the two. */
832/* print *,'two pads ','clust no. ', clno */
833 printf("two pad clusters = %d \n",(*PPP).clno ) ;
834// printf("z__(1),(2) = %f %f \n",(*PPP).z__[1],(*PPP).z__[2] ) ;
835 Int_t btest = ((*PPP).z__[1])+((*PPP).z__[2]) ;
836 if ( btest != 0 ) {
837 (*PPP).x0[1] = ((*PPP).z__[1] * (*PPP).x[1] + (*PPP).z__[2] * (*PPP).x[2]) / ((*PPP).z__[1] + (*PPP).z__[2]);
838 (*PPP).y0[1] = ((*PPP).z__[1] * (*PPP).y[1] + (*PPP).z__[2] * (*PPP).y[2]) / ((*PPP).z__[1] + (*PPP).z__[2]);
839 (*PPP).r0[1] = 0.5;
840 }
841/* Computing 2nd power */
842 d__1 = (*PPP).r0[1] = 0.5;
843 (*PPP).h0[1] = ((*PPP).z__[1] + (*PPP).z__[2]) / pi / (d__1 * d__1);
844// printf("Two pad clusters = %d %d %f %f \n",n, (*PPP).clno, (*PPP).x0[1],(*PPP).y0[1] ) ;
845 } else {
846/* guessing possible no. of Gaussians required. */
847/* assuming that if pads are separated by more than one unit, they */
848/* are possible candidates for Gaussian centers. */
849/* do this from larger values of energy ( z here ). */
850 i__1 = nn;
851 for (i__ = 1; i__ <= i__1; ++i__) {
852 iord[i__] = i__;
853 }
854 i__1 = nn;
855 for (i__ = 2; i__ <= i__1; ++i__) {
856 itest = 0;
857 ihld = iord[i__];
858 i__2 = i__ - 1;
859 for (j = 1; j <= i__2; ++j) {
860 if (itest == 0 && (*PPP).z__[iord[j]] < (*PPP).z__[ihld]) {
861 itest = 1;
862 i__3 = j;
863 for (k = i__ - 1; k >= i__3; --k) {
864 iord[k+1] = iord[k];
865 }
866 iord[j] = ihld;
867 }
868 }
869 }
870/* Gaussian centers determined. */
871 ig = 1;
872 (*PPP).icl[ig] = iord[ig];
873 i__1 = nn;
874 for (i__ = 2; i__ <= i__1; ++i__) {
875 xi = (*PPP).x[iord[i__]];
876 yi = (*PPP).y[iord[i__]];
877 idum = 0;
878 i__2 = ig;
879 for (j = 1; j <= i__2; ++j) {
880 xj = (*PPP).x[(*PPP).icl[j]];
881 yj = (*PPP).y[(*PPP).icl[j]];
882/* Computing 2nd power */
883 d__1 = xi - xj;
884/* Computing 2nd power */
885 d__2 = yi - yj;
886 rrr = TMath::Sqrt(d__1 * d__1 + d__2 * d__2);
887 if (1.01 < rrr) {
888 ++idum;
889 }
890 }
891 if (idum == ig) {
892 ++ig;
893 (*PPP).icl[ig] = iord[i__];
894 }
895 }
896 (*PPP).clno += ig;
897 printf("No. of pads = %d , No. of Gaussians = %f \n ",ig,(*PPP).clno ) ;
898 sum = 0.;
899 i__1 = ig;
900 for (i__ = 1; i__ <= i__1; ++i__) {
901 (*PPP).x0[i__] = (*PPP).x[(*PPP).icl[i__]];
902 (*PPP).y0[i__] = (*PPP).y[(*PPP).icl[i__]];
903 (*PPP).r0[i__] = .5;
904 sum += (*PPP).z__[(*PPP).icl[i__]];
905 }
906 i__1 = ig;
907 for (i__ = 1; i__ <= i__1; ++i__) {
908 if ( sum != 0 ) {
909 (*PPP).h0[i__] = (*PPP).z__[(*PPP).icl[i__]] / sum;
910 }
911 }
912 printf(" I am in ==gausspara_== Calling gausfit_ \n") ;
913 gausfit_(ig, nn, PPP );
914 printf(" I am in ==gausspara_ after gausfit_ \n") ;
915 }
916 return ;
917} /* gausspara_ */
918
919/* gausfit_ */ void gausfit_(Int_t ng, Int_t nn, PMD_dig_reco* PPP )
920{
921 Int_t i__1;
922 Float_t d__1;
923
924 Int_t iran;
925 Float_t aint;
926 Int_t i__;
927 Float_t s;
928 Float_t pi, ss, st;
929
930 pi = 3.141593;
931 printf(" ng=%d , nn=%d \n",ng,nn);
932 s = 0.;
933 i__1 = nn;
934 for (i__ = 1; i__ <= i__1; ++i__) {
935 s +=(*PPP).z__[i__];
936 }
937/* print *,s */
938 printf(" s=%d \n",s);
939 i__1 = ng;
940 for (i__ = 1; i__ <= i__1; ++i__) {
941/* Computing 2nd power */
942// printf("i__, X0,y0,r0 = %d %f %f %f \n",i__, (*PPP).x0[i__],(*PPP).y0[i__],(*PPP).r0[i__]) ;
943 d__1 = (*PPP).r0[i__];
944 (*PPP).h0[i__] = (*PPP).h0[i__] * s / pi / (d__1 * d__1);
945 }
946 st = 0.;
947 i__1 = nn;
948 for (i__ = 1; i__ <= i__1; ++i__) {
949 Int_t pbab = 1 ;
950 aintg_(ng, i__, pbab, aint, PPP);
951/* Computing 2nd power */
952 d__1 = aint - (*PPP).z__[i__];
953 if( s > 0. ) st += d__1 * d__1 / s;
954 }
955// printf(" st=%d ",st);
956 for (iran = 1; iran <= 5000; ++iran) {
957 Float_t rand;
958 i__1 = ng;
959 for (i__ = 1; i__ <= i__1; ++i__) {
960 ranmar_(rand ,PPP) ;
961 (*PPP).xx[i__] = (*PPP).x0[i__] + (rand - .5) * 1.;
962 ranmar_(rand ,PPP) ;
963 (*PPP).yy[i__] = (*PPP).y0[i__] + (rand - .5) * 1.;
964 ranmar_(rand ,PPP) ;
965 (*PPP).rr[i__] = (*PPP).r0[i__] * (((rand -.5) * .2 )+ 1.);
966/* Computing 2nd power */
967 d__1 = (*PPP).r0[i__] / (*PPP).rr[i__];
968 (*PPP).hh[i__] = (*PPP).h0[i__] * (d__1 * d__1);
969 }
970 ss = 0.;
971 i__1 = nn;
972 for (i__ = 1; i__ <= i__1; ++i__) {
973 Int_t pbab = 2 ;
974 aintg_(ng, i__, pbab, aint, PPP);
975/* Computing 2nd power */
976 d__1 = aint - (*PPP).z__[i__];
977 if(s>0.0)ss += d__1 * d__1 / s;
978 }
979 if (ss < st) {
980 i__1 = ng;
981 for (i__ = 1; i__ <= i__1; ++i__) {
982 (*PPP).x0[i__] = (*PPP).xx[i__];
983 (*PPP).y0[i__] = (*PPP).yy[i__];
984 (*PPP).r0[i__] = (*PPP).rr[i__];
985 (*PPP).h0[i__] = (*PPP).hh[i__];
986 }
987 st = ss;
988 printf(" st=%d \n",st);
989 }
990 }
991/* print 1,((*PPP).x0(i),(*PPP).(*PPP).y0(i),(*PPP).r0(i),i=1,ng) */
992/* print *,st */
993 i__1 = ng;
994 for (i__ = 1; i__ <= i__1; ++i__) {
995// printf("x0(i),y0(i),r0(i) =%f %f %f \n",(*PPP).x0[i__1],(*PPP).y0[i__1],(*PPP).r0[i__1] ) ;
996 }
997 printf(" I am ending gausfit_ \n") ;
998 return ;
999} /* gausfit_ */
1000
1001/* aintg_ */void aintg_(Int_t ng, Int_t n, Int_t bab, Float_t aint, PMD_dig_reco* PPP )
1002{
1003 Int_t i__1;
1004 Float_t d__1, d__2, d__3;
1005
1006 Int_t i__, j;
1007 Float_t xp, yp, fun;
1008
1009 aint = 0.;
1010 for (i__ = 1; i__ <= 25; ++i__) {
1011 xp = (*PPP).pts[1] [i__] + (*PPP).x[n] ;
1012 yp = (*PPP).pts[2] [i__] + (*PPP).y[n] ;
1013// printf("x[n],y[n],pts[1] [i__],[2] [i__], xp,yp =%f %f %f %f %f %f\n",(*PPP).x[n],(*PPP).y[n],(*PPP).pts[1] [i__],(*PPP).pts[2] [i__],xp,yp) ;
1014 fun = 0.;
1015 i__1 = ng;
1016 for (j = 1; j <= i__1; ++j) {
1017//printf("J,D__1,D__2,D__3= %d %f %f %f %f\n",j,(*PPP).xx[j],xp,(*PPP).yy[j],(*PPP).rr[j]) ;
1018/* Computing 2nd power */
1019 d__1 = xp - (*PPP).xx[j];
1020/* Computing 2nd power */
1021 d__2 = yp - (*PPP).yy[j];
1022/* Computing 2nd power */
1023 if( bab == 1) d__3 = (*PPP).r0[j];
1024 if( bab == 2) d__3 = (*PPP).rr[j];
1025// printf("J,D__1,D__2,D__3= %d %f %f %f \n",j,d__1,d__2,d__3) ;
1026 fun += TMath::exp(-(d__1 * d__1 + d__2 * d__2) / (d__3 * d__3)) * (*PPP).hh[j];
1027 }
1028 aint += fun * (*PPP).wts[i__];
1029// printf("FUN,aint= %f %f \n",fun,aint) ;
1030 }
1031 return ;
1032} /* aintg_ */
1033/* ranmar() >*/void ranmar_(Float_t random , PMD_dig_reco* PPP)
1034{
1035 /* Initialized data */
1036
1037 Int_t i__ = 97;
1038 Int_t j = 33;
1039 Float_t uni;
1040
1041/* ===========================================C==========================C
1042 */
1043/* Universal random number generator proposed by Marsaglia and Zaman */
1044/* in report FSU-SCRI-87-50 */
1045 uni = (*PPP).u[i__ - 1] - (*PPP).u[j - 1];
1046 if (uni < 0.) {
1047 uni += 1.;
1048 }
1049 (*PPP).u[i__ ] = uni;
1050 --i__;
1051 if (i__ == 0) { i__ = 97; }
1052 --j;
1053 if (j == 0) { j = 97; }
1054 (*PPP).c__ -= (*PPP).cd;
1055 if ((*PPP).c__ < 0.) { (*PPP).c__ += (*PPP).cm; }
1056 uni -= (*PPP).c__;
1057 if (uni <= 0.) { uni += 1.; }
1058 random = uni;
1059} /* ranmar_ */
1060