MUON/AliMUONClusterFinderv0.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 */
  19
  20 #include "AliMUONClusterFinderv0.h"
  21 #include "AliMUONSegResV1.h"
  22 //#include "TTree.h"
  23 //#include "AliRun.h"
  24 //#include <TCanvas.h>
  25 //#include <TH1.h>
  26 //#include <TPad.h>
  27 //#include <TGraph.h>
  28
  29 //----------------------------------------------------------
  30 ClassImp(AliMUONClusterFinderv0)
  31
  32     AliMUONClusterFinderv0::AliMUONClusterFinderv0
  33 (AliMUONsegmentation *segmentation, AliMUONresponse *response,
  34  TClonesArray *digits, Int_t chamber) : AliMUONClusterFinder(segmentation,response,digits,chamber)
  35 {;}
  36
  37     AliMUONClusterFinderv0::AliMUONClusterFinderv0():AliMUONClusterFinder()
  38 {;}
  39
  40 /*
  41 void AliMUONClusterFinder::AddRawCluster(const AliMUONRawCluster c)
  42 {
  43   //
  44   // Add a raw cluster copy to the list
  45   //
  46     AliMUON *MUON=(AliMUON*)gAlice->GetModule("MUON");
  47     MUON->AddRawCluster(fChamber,c);
  48     fNRawClusters++;
  49 }
  50 */
  51
  52
  53
  54 void AliMUONClusterFinderv0::Decluster(AliMUONRawCluster *cluster)
  55 {
  56 //    AliMUONdigit *dig;
  57 //    Int_t q;
  58     static int done=0;
  59     if (!done) {
  60         printf("Calling decluster\n");
  61         done=1;
  62     }
  63
  64
  65
  66     Int_t mul = cluster->fMultiplicity;
  67 //    printf("Decluster - multiplicity   %d \n",mul);
  68
  69     if (mul == 1) {
  70 //      printf("\n Nothing special for 1-clusters \n");
  71 //
  72 // Nothing special for 1-clusters
  73 //
  74         AddRawCluster(*cluster);
  75     } else if (mul ==2) {
  76 //
  77 // 2-cluster, compute offset
  78 //
  79         SetOffset(cluster);
  80         FillCluster(cluster);
  81         AddRawCluster(*cluster);
  82     } else if (mul ==3) {
  83 //
  84 // 3-cluster, check topology
  85 //      printf("\n 3-cluster, check topology \n");
  86 //
  87         if (Centered(cluster)) {
  88 //
  89 // ok, cluster is centered
  90 //          printf("\n ok, cluster is centered \n");
  91         } else {
  92 //
  93 // cluster is not centered, split into 2+1
  94 //          printf("\n cluster is not centered, split into 2+1 \n");
  95         }
  96
  97     } else {
  98         if (mul >(50-5)) printf("Decluster - multiplicity %d approaching 50\n",mul);
  99 //
 100 // 4-and more-pad clusters
 101 //
 102         SplitByLocalMaxima(cluster);
 103     } // multiplicity
 104 }
 105
 106 Int_t AliMUONClusterFinderv0::PeakOffsetAndCoordinates(Int_t DigitIndex, Float_t *X, Float_t *Y)
 107 //
 108 // Computes for which allowed offsets the digit has the highest neighbouring charge
 109 // Returns the value of the offset, and sets the pyisical coordinates of that pad
 110 // Loop on physical neighbours is specific to AliMUONsegmentationV1
 111 {
 112 Int_t nPara, offset, returnOffset=0 ;
 113 AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(DigitIndex);
 114 AliMUONsegmentationV1* seg = (AliMUONsegmentationV1*) fSegmentation;
 115 seg->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);
 116 if (nPara>1)
 117   {
 118   Float_t qMax=0;
 119   for (Int_t i=0;i<nPara; i++)
 120     {
 121     // Compute the charge on the 9 neighbouring pads
 122     // We assume that there are no pads connected in parallel in the neighbourhood
 123     //
 124     Float_t q=0;
 125     for (Int_t dx=-1;dx<2;dx++)
 126       for (Int_t dy=-1;dy<2;dy++)
 127         {
 128         if (dx==dy && dy==0)
 129           continue;
 130         Int_t padY=dig->fPadY+dy;
 131         Int_t padX=seg->Ix((Int_t) (dig->fPadX+dx+i*offset) , padY);
 132         if (fHitMap->TestHit(padX, padY)==empty)
 133            continue;
 134         AliMUONdigit* digt = (AliMUONdigit*) fHitMap->GetHit(padX,padY);
 135         q += digt->fSignal;
 136         }
 137     if (q>qMax)
 138       {
 139       returnOffset=i*offset;
 140       qMax=q;
 141       }
 142     }
 143   }
 144 fSegmentation->GetPadCxy(dig->fPadX+returnOffset,dig->fPadY,*X,*Y);
 145 return returnOffset;
 146 }
 147
 148
 149 void AliMUONClusterFinderv0::SetOffset(AliMUONRawCluster *cluster)
 150 // compute the offsets assuming that there is only one peak !
 151 {
 152 //DumpCluster(cluster);
 153 Float_t X,Y;
 154 cluster->fOffsetMap[0]=PeakOffsetAndCoordinates(cluster->fIndexMap[0],&X,&Y);
 155 for (Int_t i=1;i<cluster->fMultiplicity;i++) {
 156   AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[i]);
 157   fSegmentation->Distance2AndOffset(dig->fPadX,dig->fPadY,X,Y,&(cluster->fOffsetMap[i]));
 158   }
 159 }
 160
 161 void AliMUONClusterFinderv0::DumpCluster(AliMUONRawCluster *cluster)
 162 {
 163 printf ("other cluster\n");
 164 for (Int_t i=0; i<cluster->fMultiplicity; i++)
 165     {
 166         AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[i]);
 167         Int_t nPara, offset;
 168         fSegmentation->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);
 169
 170         printf("X %d Y %d Q %d NPara %d \n",dig->fPadX, dig->fPadY,dig->fSignal, nPara);
 171     }
 172 }
 173
 174 Bool_t AliMUONClusterFinderv0::Centered(AliMUONRawCluster *cluster)
 175 {
 176     AliMUONdigit* dig;
 177     dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
 178     Int_t ix=dig->fPadX;
 179     Int_t iy=dig->fPadY;
 180     Int_t nn;
 181     Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
 182     fSegmentation->Neighbours(ix,iy,&nn,X,Y);
 183
 184     Int_t nd=0;
 185     for (Int_t i=0; i<nn; i++) {
 186         if (fHitMap->TestHit(X[i],Y[i]) == used) {
 187             XN[nd]=X[i];
 188             YN[nd]=Y[i];
 189             nd++;
 190         }
 191     }
 192     if (nd==2) {
 193 //
 194 // cluster is centered !
 195         SetOffset(cluster);
 196         FillCluster(cluster);
 197         AddRawCluster(*cluster);
 198         return kTRUE;
 199     } else if (nd ==1) {
 200 //
 201 // Highest signal on an edge, split cluster into 2+1
 202 //
 203 // who is the neighbour ?
 204         Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
 205         Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
 206         Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
 207 //
 208 // 2-cluster
 209         AliMUONRawCluster cnew;
 210         cnew.fMultiplicity=2;
 211         cnew.fIndexMap[0]=cluster->fIndexMap[0];
 212         cnew.fIndexMap[1]=cluster->fIndexMap[i1];
 213         SetOffset(&cnew);
 214         FillCluster(&cnew);
 215         AddRawCluster(cnew);
 216 //
 217 // 1-cluster
 218         cluster->fMultiplicity=1;
 219         cluster->fIndexMap[0]=cluster->fIndexMap[i2];
 220         cluster->fIndexMap[1]=0;
 221         cluster->fIndexMap[2]=0;
 222         FillCluster(cluster);
 223         AddRawCluster(*cluster);
 224         return kFALSE;
 225     } else {
 226         printf("\n Completely screwed up %d !! \n",nd);
 227
 228     }
 229
 230         return kFALSE;
 231 }
 232
 233
 234 void AliMUONClusterFinderv0::SplitByLocalMaxima(AliMUONRawCluster *c)
 235 {
 236     AliMUONdigit* dig[50], *digt;
 237     Int_t ix[50], iy[50], q[50];
 238     Float_t x[50], y[50];
 239     Int_t i; // loops over digits
 240     Int_t j; // loops over local maxima
 241
 242     Int_t mul=c->fMultiplicity;
 243 //
 244 //  dump digit information into arrays
 245 //
 246     for (i=0; i<mul; i++)
 247     {
 248         dig[i]= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
 249         ix[i]= dig[i]->fPadX;
 250         iy[i]= dig[i]->fPadY;
 251         q[i] = dig[i]->fSignal;
 252         fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
 253     }
 254 //
 255 //  Find local maxima
 256 //
 257     Bool_t IsLocal[50];
 258     Int_t NLocal=0;
 259     Int_t AssocPeak[50];
 260     Int_t IndLocal[50];
 261     Int_t nn;
 262     Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
 263     for (i=0; i<mul; i++) {
 264         fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
 265         IsLocal[i]=kTRUE;
 266         for (j=0; j<nn; j++) {
 267             if (fHitMap->TestHit(X[j], Y[j])==empty) continue;
 268             digt=(AliMUONdigit*) fHitMap->GetHit(X[j], Y[j]);
 269             if (digt->fSignal > q[i]) {
 270                 IsLocal[i]=kFALSE;
 271                 break;
 272 //
 273 // handle special case of neighbouring pads with equal signal
 274             } else if (digt->fSignal == q[i]) {
 275                 if (NLocal >0) {
 276                     for (Int_t k=0; k<NLocal; k++) {
 277                         if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
 278                             IsLocal[i]=kFALSE;
 279                         }
 280                     }
 281                 }
 282             }
 283         } // loop over next neighbours
 284         if (IsLocal[i]) {
 285             IndLocal[NLocal]=i;
 286             // New for LYON : we guess which is the actual position of the pad hit
 287             // But this would run like that for normal chamber !
 288             c->fOffsetMap[i]=PeakOffsetAndCoordinates(c->fIndexMap[i], &(x[i]), &(y[i]));
 289             NLocal++;
 290         }
 291     } // loop over all digits
 292 //    printf("Found %d local Maxima",NLocal);
 293 //
 294 // Check if enough local clusters have been found,
 295 // if not add global maxima to the list
 296 //
 297 // But what the hell is that ? (Manu)
 298 //
 299 //    Int_t nPerMax=mul/NLocal;
 300 //    if (nPerMax > 5) {
 301 //      Int_t nGlob=mul/5-NLocal+1;
 302 //      if (nGlob > 0) {
 303 //          Int_t nnew=0;
 304 //          for (i=0; i<mul; i++) {
 305 //              if (!IsLocal[i]) {
 306 //                  IndLocal[NLocal]=i;
 307 //                  IsLocal[i]=kTRUE;
 308 //                  NLocal++;
 309 //                  nnew++;
 310 //              }
 311 //              if (nnew==nGlob) break;
 312 //          }
 313 //      }
 314 //   }
 315 //
 316 //
 317 // Associate hits to peaks
 318 //
 319     for (i=0; i<mul; i++) {
 320         //
 321         // loop on digits
 322         //
 323         Float_t dmin=1.E10;
 324         Float_t qmax=0;
 325         Int_t offset;
 326         if (IsLocal[i]) continue;
 327         for (j=0; j<NLocal; j++) {
 328             //
 329             // Loop on peaks
 330             //
 331             Int_t il=IndLocal[j];
 332 //          Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
 333 //                                +(y[i]-y[il])*(y[i]-y[il]));
 334             // Can run like that for non-Lyon chambers
 335             Float_t d = fSegmentation->Distance2AndOffset(ix[i],iy[i],x[il],y[il], &offset);
 336             Float_t ql=q[il];
 337 //
 338 // Select nearest peak
 339 //
 340             if (d<dmin) {
 341                 dmin=d;
 342                 qmax=ql;
 343                 AssocPeak[i]=j;
 344                 c->fOffsetMap[i]=offset;
 345             } else if (d==dmin) {
 346 //
 347 // If more than one take highest peak
 348 //
 349                 if (ql>qmax) {
 350                     dmin=d;
 351                     qmax=ql;
 352                     AssocPeak[i]=j;
 353                     c->fOffsetMap[i]=offset;
 354                 }
 355             } // end if
 356         } // End loop on peaks
 357     } // end loop on digits
 358 //
 359 // One cluster for each maximum
 360 //
 361     for (j=0; j<NLocal; j++) {
 362         AliMUONRawCluster cnew;
 363         cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
 364         cnew.fOffsetMap[0]=c->fOffsetMap[IndLocal[j]];
 365         cnew.fMultiplicity=1;
 366         for (i=0; i<mul; i++) {
 367             if (IsLocal[i]) continue;
 368             if (AssocPeak[i]==j) {
 369                 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
 370                 cnew.fOffsetMap[cnew.fMultiplicity]=c->fOffsetMap[i];
 371                 cnew.fMultiplicity++;
 372             }
 373         }
 374         FillCluster(&cnew);
 375         AddRawCluster(cnew);
 376     }
 377 }
 378
 379 /*
 380 void  AliMUONClusterFinderv0::FillCluster(AliMUONRawCluster* c)
 381 {
 382 //
 383 //  Completes cluster information starting from list of digits
 384 //
 385     AliMUONdigit* dig;
 386     Float_t x, y;
 387     Int_t  ix, iy;
 388
 389     c->fPeakSignal=0;
 390     c->fX=0;
 391     c->fY=0;
 392     c->fQ=0;
 393     for (Int_t i=0; i<c->fMultiplicity; i++)
 394     {
 395         dig= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
 396         ix=dig->fPadX + c.fOffsetMap[i]; // should be 0 for non-LYON
 397         iy=dig->fPadY;
 398         Int_t q=dig->fSignal;
 399 //
 400 //
 401 // peak signal and track list
 402         if (q>c->fPeakSignal) {
 403             c->fPeakSignal=0;
 404             c->fTracks[0]=dig->fTracks[0];
 405             c->fTracks[1]=dig->fTracks[1];
 406             c->fTracks[2]=dig->fTracks[2];
 407         }
 408 //
 409 // centre of gravity
 410         fSegmentation->GetPadCxy(ix, iy, x, y);
 411         c->fX += q*x;
 412         c->fY += q*y;
 413         c->fQ += q;
 414     }
 415     c->fX/=c->fQ;
 416 // Not valid for inclined tracks in X !!! (Manu)
 417 //    c->fX=fSegmentation->GetAnod(c->fX);
 418     c->fY/=c->fQ;
 419 //
 420 //  apply correction to the coordinate along the anode wire
 421 //
 422     if (fCogCorr) {
 423         x=c->fX;
 424         y=c->fY;
 425         fSegmentation->GetPadIxy(x, y, ix, iy);
 426         fSegmentation->GetPadCxy(ix, iy, x, y);
 427         Float_t YonPad=(c->fY-y)/fSegmentation->Dpy();
 428         c->fY=y-fCogCorr->Eval(YonPad, 0, 0);
 429     }
 430
 431 }
 432 */
 433 /*
 434 void  AliMUONClusterFinder::FindCluster(Int_t i, Int_t j, AliMUONRawCluster &c){
 435 //
 436 //  Find clusters
 437 //
 438 //
 439 //  Add i,j as element of the cluster
 440 //
 441
 442     Int_t idx = fHitMap->GetHitIndex(i,j);
 443     AliMUONdigit* dig = (AliMUONdigit*) fHitMap->GetHit(i,j);
 444     Int_t q=dig->fSignal;
 445     if (q > TMath::Abs(c.fPeakSignal)) {
 446         c.fPeakSignal=q;
 447         c.fTracks[0]=dig->fTracks[0];
 448         c.fTracks[1]=dig->fTracks[1];
 449         c.fTracks[2]=dig->fTracks[2];
 450     }
 451 //
 452 //  Make sure that list of digits is ordered
 453 //
 454     Int_t mu=c.fMultiplicity;
 455     c.fIndexMap[mu]=idx;
 456
 457     if (mu > 0) {
 458         for (Int_t ind=mu-1; ind>=0; ind--) {
 459             Int_t ist=(c.fIndexMap)[ind];
 460             Int_t ql=((AliMUONdigit*)fDigits
 461                       ->UncheckedAt(ist))->fSignal;
 462             if (q>ql) {
 463                 c.fIndexMap[ind]=idx;
 464                 c.fIndexMap[ind+1]=ist;
 465             } else {
 466                 break;
 467             }
 468         }
 469     }
 470
 471     c.fMultiplicity++;
 472
 473     if (c.fMultiplicity >= 50 ) {
 474         printf("FindCluster - multiplicity >50  %d \n",c.fMultiplicity);
 475         c.fMultiplicity=50-1;
 476     }
 477
 478 // Prepare center of gravity calculation
 479     Float_t x, y;
 480     fSegmentation->GetPadCxy(i, j, x, y);
 481     c.fX += q*x;
 482     c.fY += q*y;
 483     c.fQ += q;
 484 // Flag hit as taken
 485     fHitMap->FlagHit(i,j);
 486 //
 487 //  Now look recursively for all neighbours
 488 //
 489     Int_t nn;
 490     Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
 491     fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
 492     for (Int_t in=0; in<nn; in++) {
 493         Int_t ix=Xlist[in];
 494         Int_t iy=Ylist[in];
 495         if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
 496     }
 497 }
 498 */
 499
 500 //_____________________________________________________________________________
 501
 502 void AliMUONClusterFinderv0::FindRawClusters()
 503 {
 504   //
 505   // simple MUON cluster finder from digits -- finds neighbours and
 506   // fill the tree with raw clusters
 507   //
 508     if (!fNdigits) return;
 509
 510     fHitMap = new AliMUONHitMapA1(fSegmentation, fDigits);
 511
 512     AliMUONdigit *dig;
 513
 514     int ndig;
 515     int nskip=0;
 516
 517     fHitMap->FillHits();
 518     for (ndig=0; ndig<fNdigits; ndig++) {
 519         dig = (AliMUONdigit*)fDigits->UncheckedAt(ndig);
 520         Int_t i=dig->fPadX;
 521         Int_t j=dig->fPadY;
 522         if (fHitMap->TestHit(i,j)==used ||fHitMap->TestHit(i,j)==empty) {
 523             nskip++;
 524             continue;
 525         }
 526         AliMUONRawCluster c;
 527         c.fMultiplicity=0;
 528 //      c.fPeakSignal=dig->fSignal;
 529 //      c.fTracks[0]=dig->fTracks[0];
 530 //      c.fTracks[1]=dig->fTracks[1];
 531 //      c.fTracks[2]=dig->fTracks[2];
 532         c.fPeakSignal=0;
 533         FindCluster(i,j, c);
 534         // center of gravity
 535         c.fX /= c.fQ;
 536         c.fX=fSegmentation->GetAnod(c.fX);
 537         c.fY /= c.fQ;
 538 //
 539 //  apply correction to the coordinate along the anode wire
 540 //
 541
 542
 543     if (fCogCorr) {
 544         Int_t ix,iy;
 545         Float_t x=c.fX;
 546         Float_t y=c.fY;
 547         fSegmentation->GetPadIxy(x, y, ix, iy);
 548         fSegmentation->GetPadCxy(ix, iy, x, y);
 549         Float_t YonPad=(c.fY-y)/fSegmentation->Dpy();
 550         c.fY=y-fCogCorr->Eval(YonPad,0,0);
 551     }
 552 //
 553 //      Analyse cluster and decluster if necessary
 554 //
 555         Decluster(&c);
 556 //
 557 //
 558 //
 559 //      reset Cluster object
 560         for (int k=0;k<c.fMultiplicity;k++) {
 561             c.fIndexMap[k]=0;
 562             c.fOffsetMap[k]=0;
 563         }
 564         c.fMultiplicity=0;
 565     } // end loop ndig
 566     delete fHitMap;
 567 }
 568
 569 /*
 570
 571 void AliMUONClusterFinder::
 572 CalibrateCOG()
 573 {
 574     Float_t x[5];
 575     Float_t y[5];
 576     Int_t n, i;
 577     TF1 func;
 578
 579     if (fSegmentation) {
 580         fSegmentation->GiveTestPoints(n, x, y);
 581         for (i=0; i<n; i++) {
 582             Float_t xtest=x[i];
 583             Float_t ytest=y[i];
 584             SinoidalFit(xtest, ytest, func);
 585         }
 586         fCogCorr = new TF1(func);
 587     }
 588 }
 589 */
 590 /*
 591
 592 void AliMUONClusterFinder::
 593 SinoidalFit(Float_t x, Float_t y, TF1 &func)
 594 {
 595 //
 596     static Int_t count=0;
 597     char canvasname[3];
 598     count++;
 599     sprintf(canvasname,"c%d",count);
 600
 601 // MANU : without const, error on HP
 602     const Int_t ns=101;
 603     Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
 604     Float_t xsig[ns], ysig[ns];
 605
 606     AliMUONsegmentation *segmentation=fSegmentation;
 607
 608     Int_t ix,iy;
 609     segmentation->GetPadIxy(x,y,ix,iy);
 610     segmentation->GetPadCxy(ix,iy,x,y);
 611     Int_t isec=segmentation->Sector(ix,iy);
 612 // Pad Limits
 613     Float_t xmin = x-segmentation->GetRealDpx(isec)/2;
 614     Float_t ymin = y-segmentation->Dpy()/2;
 615 //
 616 //      Integration Limits
 617     Float_t dxI=fResponse->Nsigma()*fResponse->ChwX();
 618     Float_t dyI=fResponse->Nsigma()*fResponse->ChwY();
 619
 620 //
 621 //  Scanning
 622 //
 623     Int_t i;
 624     Float_t qp;
 625 //
 626 //  y-position
 627     Float_t yscan=ymin;
 628     Float_t dy=segmentation->Dpy()/(ns-1);
 629
 630     for (i=0; i<ns; i++) {
 631 //
 632 //      Pad Loop
 633 //
 634         Float_t sum=0;
 635         Float_t qcheck=0;
 636         segmentation->SigGenInit(x, yscan, 0);
 637
 638         for (segmentation->FirstPad(x, yscan, dxI, dyI);
 639              segmentation->MorePads();
 640              segmentation->NextPad())
 641         {
 642             qp=fResponse->IntXY(segmentation);
 643             qp=TMath::Abs(qp);
 644 //
 645 //
 646             if (qp > 1.e-4) {
 647                 qcheck+=qp;
 648                 Int_t ixs=segmentation->Ix();
 649                 Int_t iys=segmentation->Iy();
 650                 Float_t xs,ys;
 651                 segmentation->GetPadCxy(ixs,iys,xs,ys);
 652                 sum+=qp*ys;
 653             }
 654         } // Pad loop
 655         Float_t ycog=sum/qcheck;
 656         yg[i]=(yscan-y)/segmentation->Dpy();
 657         yrg[i]=(ycog-y)/segmentation->Dpy();
 658         ysig[i]=ycog-yscan;
 659         yscan+=dy;
 660     } // scan loop
 661 //
 662 //  x-position
 663     Float_t xscan=xmin;
 664     Float_t dx=segmentation->GetRealDpx(isec)/(ns-1);
 665
 666     for (i=0; i<ns; i++) {
 667 //
 668 //      Pad Loop
 669 //
 670         Float_t sum=0;
 671         Float_t qcheck=0;
 672         segmentation->SigGenInit(xscan, y, 0);
 673
 674         for (segmentation->FirstPad(xscan, y, dxI, dyI);
 675              segmentation->MorePads();
 676              segmentation->NextPad())
 677         {
 678             qp=fResponse->IntXY(segmentation);
 679             qp=TMath::Abs(qp);
 680 //
 681 //
 682             if (qp > 1.e-2) {
 683                 qcheck+=qp;
 684                 Int_t ixs=segmentation->Ix();
 685                 Int_t iys=segmentation->Iy();
 686                 Float_t xs,ys;
 687                 segmentation->GetPadCxy(ixs,iys,xs,ys);
 688                 sum+=qp*xs;
 689             }
 690         } // Pad loop
 691         Float_t xcog=sum/qcheck;
 692         xcog=segmentation->GetAnod(xcog);
 693
 694         xg[i]=(xscan-x)/segmentation->GetRealDpx(isec);
 695         xrg[i]=(xcog-x)/segmentation->GetRealDpx(isec);
 696         xsig[i]=xcog-xscan;
 697         xscan+=dx;
 698     }
 699
 700    TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
 701    TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 702    TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 703    TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 704    TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 705    pad11->SetFillColor(11);
 706    pad12->SetFillColor(11);
 707    pad13->SetFillColor(11);
 708    pad14->SetFillColor(11);
 709    pad11->Draw();
 710    pad12->Draw();
 711    pad13->Draw();
 712    pad14->Draw();
 713    TGraph *graphx = new TGraph(ns,xg ,xsig);
 714    TGraph *graphxr= new TGraph(ns,xrg,xsig);
 715    TGraph *graphy = new TGraph(ns,yg ,ysig);
 716    TGraph *graphyr= new TGraph(ns,yrg,ysig);
 717 //
 718 // Creates a Root function based on function sinoid above
 719 // and perform the fit
 720 //
 721    Double_t sinoid(Double_t *x, Double_t *par);
 722    TF1 *sinoidf = new TF1("sinoidf",sinoid,0.5,0.5,5);
 723    graphyr->Fit("sinoidf","V");
 724    sinoidf->Copy(func);
 725    func.Eval(0,0,0);
 726 //
 727    pad11->cd();
 728    graphx->SetFillColor(42);
 729    graphx->SetMarkerColor(4);
 730    graphx->SetMarkerStyle(21);
 731    graphx->Draw("AC");
 732    graphx->GetHistogram()->SetXTitle("x on pad");
 733    graphx->GetHistogram()->SetYTitle("xcog-x");
 734
 735
 736    pad12->cd();
 737    graphxr->SetFillColor(42);
 738    graphxr->SetMarkerColor(4);
 739    graphxr->SetMarkerStyle(21);
 740    graphxr->Draw("AP");
 741    graphxr->GetHistogram()->SetXTitle("xcog on pad");
 742    graphxr->GetHistogram()->SetYTitle("xcog-x");
 743
 744
 745    pad13->cd();
 746    graphy->SetFillColor(42);
 747    graphy->SetMarkerColor(4);
 748    graphy->SetMarkerStyle(21);
 749    graphy->Draw("AF");
 750    graphy->GetHistogram()->SetXTitle("y on pad");
 751    graphy->GetHistogram()->SetYTitle("ycog-y");
 752
 753
 754
 755    pad14->cd();
 756    graphyr->SetFillColor(42);
 757    graphyr->SetMarkerColor(4);
 758    graphyr->SetMarkerStyle(21);
 759    graphyr->Draw("AF");
 760    graphyr->GetHistogram()->SetXTitle("ycog on pad");
 761    graphyr->GetHistogram()->SetYTitle("ycog-y");
 762
 763    c1->Update();
 764
 765 }
 766 */
 767 /*
 768 Double_t sinoid(Double_t *x, Double_t *par)
 769 {
 770     Double_t arg = -2*TMath::Pi()*x[0];
 771     Double_t fitval= par[0]*TMath::Sin(arg)+
 772         par[1]*TMath::Sin(2*arg)+
 773         par[2]*TMath::Sin(3*arg)+
 774         par[3]*TMath::Sin(4*arg)+
 775         par[4]*TMath::Sin(5*arg);
 776     return fitval;
 777  }
 778 */
 779
 780
 781
 782
 783
 784
 785