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
5 void PMD_dig_reco (Int_t evNumber=1)
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.
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
17 /////////////////////////////////////////////////////////////////////////
18 // =========== Variables for Reconstruction==============================
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;
28 Int_t iVol4,iVol5 , ixpad,iypad;
39 long double r,eta,px,py,pz;
40 Float_t phi,Eparticle,theta1,theta;
44 fp1=fopen("kine_vertex.out","w");
45 fp2=fopen("kine_pmdhit.out","w");
46 fp4=fopen("det1.out","w");
48 // Dynamical link some shared libs
49 if (gClassTable->GetID("AliRun") < 0) {
50 gROOT->LoadMacro("loadlibs.C");
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");
58 // Get AliRun object from file or create it if not on file
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");
65 //Loop over the requested number of events
66 for (Int_t iEvt = 0; iEvt < evNumber; iEvt++) {
68 printf("evt no %d \n",iEvt);
73 // Import the Kine and Hits Trees for the event evNumber in the file
74 Int_t nparticles = gAlice->GetEvent(iEvt);
76 printf ("N of particles %d\n",nparticles);
78 if (nparticles <= 0) return;
80 Int_t nevt = iEvt + 1;
82 // ----********---- vertex track search BEGIN
84 //Loop over the requested number of events
85 // Import the Kine and Hits Trees for the event evNumber in the file
87 TObjArray *Particles = gAlice->Particles();
89 Float_t pl,cutmax,etamin,etamax,pmom,smin,smax;
94 // for (Int_t it = 0; it < 100000; it++)trk_Mult[it]=0.;
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);
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;
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);
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);
131 // 1999 file-- Get pointers to Alice detectors and Hits containers
132 //------********----- vertex track ENDs
134 Float_t x,y,z,mass,e;
143 const Int_t Npad2 = 72;
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();
158 TTree *TH = gAlice->TreeH();
159 // TTree *TK = gAlice->TreeK();
161 Int_t ntracks = TH->GetEntries();
162 printf("No.of Tracksin the TreeH %d \n", ntracks);
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 ;
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.);
182 // ================================================================
184 // =======================================================
185 // Start loop on tracks in the hits containers
186 for (Int_t track=0; track<ntracks;track++) {
189 //changes introduced ******
191 nbytes += TH->GetEvent(track);
192 // mbytes += TK->GetEvent(track);
193 // printf("%d ******** %d \n" ,track,nbytes);
194 // ======>Histogram 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();
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;
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);
215 Eparticle = particle->Energy();
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;
226 // if(TMath::Abs(eta) > 2.3 && TMath::Abs(eta) < 3.5){
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];
238 ypad = Vol2s/Npad2 + 1;
240 xpad = Vol2-ypad1*Npad2;
241 Edep = pmdHit->fEnergy;
246 // module [Vol4] [Vol5] [xpad] [ypad] = Edep ;
247 if(Edep>0)printf(" VOL4 %d %d %d %d %f \n",iVol4,iVol5,ixpad,iypad,Edep) ;
249 Modnum ->Fill(iVol5) ;
250 xycells->Fill(ixpad,iypad);
251 xyhits->Fill(xPos,yPos);
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);
261 // Closing Eta loop of the incident particles.
271 // TCanvas *c1 = new TCanvas("c1","PMD Simulated digits(Edep) ",400,10,600,700);
274 // gPad->SetFillColor(33);
275 // edp->SetFillColor(42);
278 // gPad->SetFillColor(33);
279 // Modnum->SetFillColor(41);
281 // c1->Print("baba.ps");
282 // TCanvas *c2 = new TCanvas("c2","PMD Simulated digits(Edep) ");
283 // gPad->SetFillColor(33);
284 // xycells->SetFillColor(42);
286 // c2->Print("baba1.ps");
287 // TCanvas *c3 = new TCanvas("c3","PMD Simulated digits(Edep) ");
288 // gPad->SetFillColor(33);
289 // xyhits->SetFillColor(42);
291 // c3->Print("baba2.ps");
294 // =========== RECONSTRUCTION ======================================
295 void pmdreco_(PMD_dig_reco* PPP) {
297 sqrth = sqrt(3.) / 2.;
299 printf(" I am in MAIN calling RanMAR & Intpts_ \n") ;
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;
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 ) ;
324 // Ending Reading pads in the Box .
326 // Ending Reading pads for 1-CPV,2-PMD.
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");
333 printf(" I am in MAIN after order for CPV(det=1) \n") ;
334 printf(" I am in MAIN calling ==crclust== \n") ;
336 printf(" I am in MAIN after ==crclust== & now call prints \n") ;
338 printf(" I am in MAIN Doing Clustering for PMD. \n") ;
339 printf (" I am in MAIN calling ==order== for PMD(det=2) \n");
341 printf(" I am in MAIN after order for PMD(det=2) \n") ;
342 printf(" I am in MAIN calling ==crclust== \n") ;
344 printf(" I am in MAIN after ==crclust== & now call prints \n") ;
347 for (Int_t i__ = 1; i__ <= 72; ++i__) {
348 for (j = 1; j <= 72; ++j) {
349 if ((*PPP).infocl[1] [i__] [j] == 1) {
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;
357 for (i__ = 1; i__ <= 72; ++i__) {
358 for (j = 1; j <= 72; ++j) {
359 if ((*PPP).infocl[1] [i__] [j] == 2) {
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;
368 /* better clustering ( Gaussian fit ) */
373 /* >>>>>>>>>>>> */void rmarin_(PMD_dig_reco* PPP)
375 Int_t i__, j, k, l, m;
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 */
381 /* To get the standard values in Marsaglia's paper, IJ=1802, KL=9373 */
384 if (ij <= 0 || ij >= 31328) {
387 if (kl <= 0 || kl >= 30081) {
390 printf("I am in RanMAR \n");
391 i__ = ij / 177 % 177 + 2;
393 k = kl / 169 % 178 + 1;
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) {
399 for (Int_t jj = 1; jj <= 24; ++jj) {
400 m = (((i__ * j) % 179) * k) % 179;
404 l = l * 53 + 1 % 169;
405 if (l * m % 64 >= 32) {
410 (*PPP).u[ii - 1] = s;
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) ;
418 printf( " +++++++ stop in rmarin: wrong params! ");
421 /*>>>>>>>>>>>>>>>> */void order_(Int_t jdet,PMD_dig_reco* PPP)
423 Int_t nmx,i__1, i__2, i__3;
426 Int_t idum, i__, j, k, l;
429 Int_t itest, i1, i2, id ;
431 printf("I am in order for DET=%d \n",jdet);
433 for (i1 = 1; i1 <= 72; ++i1) {
434 for (i2 = 1; i2 <= 72; ++i2) {
435 i__ = (i2 - 1) * 72 + i1;
438 if(jdet==1) d1[i__] = (*PPP).det1 [i1] [i2] ;
439 if(jdet==2) d1[i__] = (*PPP).det2 [i1] [i2] ;
443 for (j = 2; j <= nmx; ++j) {
448 for (k = 1; k <= i__2; ++k) {
449 if (adum > d1[k ] && itest == 0) {
452 for (l = j - 1; l >= i__3; --l) {
454 iord1 [l+1] = iord1 [l];
462 for (i__ = 1; i__ <= i__1; ++i__) {
465 i1 = id - (i2 - 1) * 72;
466 (*PPP).iorder[1] [i__] = i1;
467 (*PPP).iorder[2] [i__] = i2;
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__]);
473 /* crclust_ */ void crclust_(Int_t ndet, PMD_dig_reco* PPP)
475 Float_t d__ [73] [73] ;
477 Int_t neibx[7] = { 0,1,0,-1,-1,0,1 };
478 Int_t neiby[7] = { 0,0,1,1,0,-1,-1 };
483 Int_t icld, idum, i__, j, k;
484 Int_t nb, id1, id2, jd1, jd2, icle, icn, itr, nmx ;
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]) ;
489 for (j = 1; j <= 72; ++j) {
490 for (k = 1; k <= 72; ++k) {
492 if(ndet==1) d__ [j] [k] = (*PPP).det1 [j] [k] ;
493 if(ndet==2) d__ [j] [k] = (*PPP).det2 [j] [k] ;
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;
507 /* compute number of pads with nonzero count ( nmx1 ) */
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) {
514 printf(" nmx1 -> %d %d %f %d \n",id1,id2, d__ [id1] [id2],(*PPP).nmx1 ) ;
517 printf(" nmx1= %d \n",(*PPP).nmx1) ;
518 /* crude clustering. start with pad with largest count ( cluster 1 )
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 */
529 /* infocl(1,i1,i2) =1 for primary ( peak ) =2 for secondary ( */
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;
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] ) ;
544 /* start checking with next largest count and so on */
546 for (i__ = 2; i__ <= i__1; ++i__) {
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
553 /* cluster ( itr = 1 ) Otherwise itr = 0. */
554 if ((*PPP).infocl [1] [id1] [id2] > 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
563 /* icn = 1. Otherwise icn = 0. */
565 if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
566 if ((*PPP).infocl [1] [jd1] [jd2] != 0) {
568 icld = (*PPP).infocl [2] [jd1] [jd2] ;
572 /* itr = 0 and icn = 0 means new cluster. Store the cluster in
575 if (TMath::Abs(icn) <= 0.1) {
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]);
593 /* itr = 0 and icn = 1 means the pad belongs to cluster ic
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;
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;
625 printf(" ndet=%d,nmx1= %d , icle= %d icld= %d \n",ndet,(*PPP).nmx1,icle,icld) ;
632 /* Subroutine */ void intpts_(PMD_dig_reco* PPP)
634 Int_t i__1, i__2, i__3;
637 Float_t x, y, dx, dy;
648 for (ix = -ixpt; ix <= i__1; ++ix) {
652 wtx = (3. - idx - (TMath::Abs(i__2) ))/ 3. * dx;
653 yl = 1. / TMath::Sqrt(3.) * (1. - TMath::Abs(x));
657 for (iy = -ixpt; iy <= i__2; ++iy) {
661 wty = (3. - idy - (TMath::Abs(i__3))) / 3. * dy;
663 (*PPP).pts[1] [ipt] = x;
664 (*PPP).pts[2] [ipt] = y;
665 (*PPP).wts[ipt] = wtx * wty;
668 printf(" Exiting from intpts_ %f %f \n", ipt,(*PPP).wts[ipt] );
672 /* prints_ */ void prints_(Int_t ndet, PMD_dig_reco* PPP)
679 fp1=fopen("CPV_prime.out","w");
680 fp2=fopen("CPV_secondary.out","w");
685 fp1=fopen("PMD_prime.out","w");
686 fp2=fopen("PMD_secondary.out","w");
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] ) ;
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] ) ;
701 /* Subroutine */void arrange_(PMD_dig_reco* PPP)
703 Int_t i__1, i__2, i__3;
705 Int_t ihld1, ihld2, ihld3, i__, j, k, itest;
708 for (i__ = 2; i__ <= i__1; ++i__) {
710 ihld1 = (*PPP).inford [1] [i__];
711 ihld2 = (*PPP).inford [2] [i__];
712 ihld3 = (*PPP).inford [3] [i__];
714 for (j = 1; j <= i__2; ++j) {
715 if (itest == 0 && (*PPP).inford [1] [j] > ihld1) {
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] ;
723 (*PPP).inford [1] [j] = ihld1 ;
724 (*PPP).inford [2] [j] = ihld2 ;
725 (*PPP).inford [3] [j] = ihld3 ;
729 printf(" Exiting from arrange \n") ;
733 /* refclust_ */void refclust_(PMD_dig_reco* PPP)
737 Float_t d__ [73] [73] ;
738 Int_t npad[201], i__, j, k;
739 Int_t id, id1, id2, iclb;
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] ;
751 /* The No.of pads in each cluster = 200 */
752 for (i__ = 1; i__ <= 200; ++i__) {
755 for (i__ = 1; i__ <= 5000; ++i__) {
758 /* number of pads in each cluster determined. */
760 for (i__ = 1; i__ <= i__1; ++i__) {
761 id = (*PPP).inford[1] [i__];
763 Int_t nclid = (*PPP).ncl[id] ;
764 printf("id,ncl= %d %d\n",id,(*PPP).ncl[id]) ;
767 /* number of clusters determined. */
768 for (i__ = 1; i__ <= 5000; ++i__) {
769 id = (*PPP).ncl[i__];
775 printf("iclb= %d \n",iclb) ;
777 /* define coordinates of each pad in a cluster and call gausspara */
778 /* to determine Gaussian parameters. ( should it be eliminated ? ) */
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) {
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];
791 printf(" I am in refclust calling ==gausspara_== \n") ;
792 gausspara_(i__2,PPP );
793 printf(" I am in refclust after ==gausspara_== \n") ;
795 printf("clno= %d \n",(*PPP).clno);
799 /* gausspara_ */ void gausspara_(Int_t n,PMD_dig_reco* PPP)
801 Int_t i__1, i__2, i__3;
804 Int_t ihld, iord[501], idum, i__, j, k;
807 Float_t pi, xi, yi, xj, yj, rrr;
810 /* single Gaussian of the form exp(-((x-x_0)^2+(y-y_0)^2)/r_0^2) */
812 printf(" I am in ==gausspara_== %d \n" , n) ;
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];
823 /* Computing 2nd power */
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) {
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]) ;
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]);
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] ) ;
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 ). */
851 for (i__ = 1; i__ <= i__1; ++i__) {
855 for (i__ = 2; i__ <= i__1; ++i__) {
859 for (j = 1; j <= i__2; ++j) {
860 if (itest == 0 && (*PPP).z__[iord[j]] < (*PPP).z__[ihld]) {
863 for (k = i__ - 1; k >= i__3; --k) {
870 /* Gaussian centers determined. */
872 (*PPP).icl[ig] = iord[ig];
874 for (i__ = 2; i__ <= i__1; ++i__) {
875 xi = (*PPP).x[iord[i__]];
876 yi = (*PPP).y[iord[i__]];
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 */
884 /* Computing 2nd power */
886 rrr = TMath::Sqrt(d__1 * d__1 + d__2 * d__2);
893 (*PPP).icl[ig] = iord[i__];
897 printf("No. of pads = %d , No. of Gaussians = %f \n ",ig,(*PPP).clno ) ;
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__]];
904 sum += (*PPP).z__[(*PPP).icl[i__]];
907 for (i__ = 1; i__ <= i__1; ++i__) {
909 (*PPP).h0[i__] = (*PPP).z__[(*PPP).icl[i__]] / sum;
912 printf(" I am in ==gausspara_== Calling gausfit_ \n") ;
913 gausfit_(ig, nn, PPP );
914 printf(" I am in ==gausspara_ after gausfit_ \n") ;
919 /* gausfit_ */ void gausfit_(Int_t ng, Int_t nn, PMD_dig_reco* PPP )
931 printf(" ng=%d , nn=%d \n",ng,nn);
934 for (i__ = 1; i__ <= i__1; ++i__) {
938 printf(" s=%d \n",s);
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);
948 for (i__ = 1; i__ <= i__1; ++i__) {
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;
955 // printf(" st=%d ",st);
956 for (iran = 1; iran <= 5000; ++iran) {
959 for (i__ = 1; i__ <= i__1; ++i__) {
961 (*PPP).xx[i__] = (*PPP).x0[i__] + (rand - .5) * 1.;
963 (*PPP).yy[i__] = (*PPP).y0[i__] + (rand - .5) * 1.;
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);
972 for (i__ = 1; i__ <= i__1; ++i__) {
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;
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__];
988 printf(" st=%d \n",st);
991 /* print 1,((*PPP).x0(i),(*PPP).(*PPP).y0(i),(*PPP).r0(i),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] ) ;
997 printf(" I am ending gausfit_ \n") ;
1001 /* aintg_ */void aintg_(Int_t ng, Int_t n, Int_t bab, Float_t aint, PMD_dig_reco* PPP )
1004 Float_t d__1, d__2, d__3;
1007 Float_t xp, yp, fun;
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) ;
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];
1028 aint += fun * (*PPP).wts[i__];
1029 // printf("FUN,aint= %f %f \n",fun,aint) ;
1033 /* ranmar() >*/void ranmar_(Float_t random , PMD_dig_reco* PPP)
1035 /* Initialized data */
1041 /* ===========================================C==========================C
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];
1049 (*PPP).u[i__ ] = uni;
1051 if (i__ == 0) { i__ = 97; }
1053 if (j == 0) { j = 97; }
1054 (*PPP).c__ -= (*PPP).cd;
1055 if ((*PPP).c__ < 0.) { (*PPP).c__ += (*PPP).cm; }
1057 if (uni <= 0.) { uni += 1.; }