MUON/AliMUONContourMaker.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 // $Id$
  17
  18 ///
  19 /// Maker/merger of contours. Can create contour from a set polygons or
  20 /// merger a set of contours into a single one.
  21 ///
  22 /// This is based on (one of the) algorithm found in
  23 /// Diane L. Souvaine and Iliana Bjorling-Sachs,
  24 /// Proceedings of the IEEE, Vol. 80, No. 9, September 1992, p. 1449
  25 ///
  26 /// Note that besides the AliMUON prefix, nothing is really MUON specific
  27 /// in this class...
  28 ///
  29 /// \author Laurent Aphecetche, Subatech
  30 ///
  31
  32 #include "AliMUONContourMaker.h"
  33
  34 #include "AliCodeTimer.h"
  35 #include "AliLog.h"
  36 #include "AliMUONContour.h"
  37 #include "AliMUONPointWithRef.h"
  38 #include "AliMUONPolygon.h"
  39 #include "AliMUONSegment.h"
  40 #include "AliMUONSegmentTree.h"
  41 #include "Riostream.h"
  42 #include "TArrayD.h"
  43 #include "TMath.h"
  44 #include <cassert>
  45 #include "TArrayI.h"
  46
  47 /// \cond CLASSIMP
  48 ClassImp(AliMUONContourMaker)
  49 /// \endcond
  50
  51 namespace
  52 {
  53   void PrintSegments(const TObjArray& array)
  54   {
  55     TIter next(&array);
  56     AliMUONSegment* s;
  57     Int_t i(0);
  58     while ( ( s = static_cast<AliMUONSegment*>(next()) ) )
  59     {
  60       cout << Form("i=%d %s",i,s->AsString()) << endl;
  61       ++i;
  62     }
  63   }
  64 }
  65
  66 //_____________________________________________________________________________
  67 AliMUONContourMaker::AliMUONContourMaker()
  68 {
  69 /// Default constructor
  70 }
  71
  72
  73 //_____________________________________________________________________________
  74 AliMUONContourMaker::~AliMUONContourMaker()
  75 {
  76 /// Destructor
  77 }
  78
  79 //_____________________________________________________________________________
  80 AliMUONContour*
  81 AliMUONContourMaker::CreateContour(const TObjArray& polygons, const char* name) const
  82 {
  83   /// Create the contour of the polygon array
  84   /// and get back the intermediate verticals and horizontal segments
  85   /// both arrays are arrays of AliMUONSegment objects.
  86
  87   AliCodeTimerAuto("",0);
  88
  89   if ( polygons.IsEmpty() ) return 0x0; // protection against user error...
  90
  91   // Sanity check : insure that all polygons are oriented counter-clockwise
  92   TIter next(&polygons);
  93   AliMUONPolygon* pol;
  94   while ( ( pol = static_cast<AliMUONPolygon*>(next()) ) )
  95   {
  96     if ( !pol->IsCounterClockwiseOriented() )
  97     {
  98       AliError(Form("Got a clockwise oriented polygon in CreateContour(%s). That's not OK !",name));
  99       StdoutToAliError(polygons.Print());
 100       return 0x0;
 101     }
 102   }
 103
 104   AliMUONContour* contour(0x0);
 105
 106   if ( polygons.GetLast() == 0 )
 107   {
 108     AliCodeTimerAuto("Trivial case",1);
 109     contour = new AliMUONContour(name);
 110     pol = static_cast<AliMUONPolygon*>(polygons.First());
 111     contour->Add(*pol);
 112     contour->AssertOrientation();
 113     return contour;
 114   }
 115
 116   TObjArray polygonVerticalEdges; // arrray of AliMUONSegment objects
 117   polygonVerticalEdges.SetOwner(kTRUE);
 118   // get vertical edges of input polygons
 119   GetVerticalEdges(polygons,polygonVerticalEdges);
 120
 121   // sort them in ascending x order
 122   // if same x, insure that left edges are before right edges
 123   // within same x, order by increasing bottommost y (see AliMUONSegment::Compare method)
 124   polygonVerticalEdges.Sort();
 125
 126   if ( polygonVerticalEdges.GetLast()+1 < 2 )
 127   {
 128     AliError(Form("Got too few edges here for createContour %s",name));
 129     polygons.Print();
 130     TObject* o(0x0);
 131     o->Print();
 132   }
 133
 134   // Find the vertical edges of the merged contour. This is the meat of the algorithm...
 135   TObjArray contourVerticalEdges;
 136   contourVerticalEdges.SetOwner(kTRUE);
 137   Sweep(polygonVerticalEdges,contourVerticalEdges);
 138
 139   TObjArray horizontals;
 140   horizontals.SetOwner(kTRUE);
 141   VerticalToHorizontal(contourVerticalEdges,horizontals);
 142
 143   contour = FinalizeContour(contourVerticalEdges,horizontals);
 144
 145   if ( contour && name ) contour->SetName(name);
 146
 147   return contour;
 148 }
 149
 150 //_____________________________________________________________________________
 151 AliMUONContour*
 152 AliMUONContourMaker::FinalizeContour(const TObjArray& verticals,
 153                                      const TObjArray& horizontals) const
 154 {
 155   /// For a list of vertical and horizontal edges, we build the final
 156   /// contour object.
 157
 158   AliCodeTimerAuto("",0);
 159
 160   TObjArray all; // array of AliMUONSegment
 161   TObjArray inorder; // array of AliMUONSegment
 162
 163   all.SetOwner(kFALSE);
 164   inorder.SetOwner(kFALSE);
 165
 166   for ( Int_t i = 0; i <= verticals.GetLast(); ++i )
 167   {
 168     all.Add(verticals.UncheckedAt(i));
 169     all.Add(horizontals.UncheckedAt(i));
 170   }
 171
 172   TArrayI alreadyAdded(all.GetLast()+1);
 173   alreadyAdded.Reset();
 174
 175   Int_t i(0);
 176
 177   AliMUONContour* contour = new AliMUONContour;
 178
 179   int total(0);
 180
 181   while ( !all.IsEmpty() )
 182   {
 183     total++;
 184
 185     if ( total > 1000 )
 186     {
 187       AliError("Total 1000 reached !!!!");
 188       return 0x0;
 189     }
 190
 191     AliMUONSegment* si = static_cast<AliMUONSegment*>(all.UncheckedAt(i));
 192     inorder.Add(si);
 193     alreadyAdded[i] = 1;
 194     const AliMUONSegment* all0 = static_cast<const AliMUONSegment*>(all.First());
 195     if ( i != 0 && AliMUONSegment::AreEqual(si->EndX(),all0->StartX()) && AliMUONSegment::AreEqual(si->EndY(),all0->StartY()) )
 196     {
 197       Int_t n(-1);
 198
 199       AliMUONPolygon polygon(inorder.GetLast()+2);
 200
 201       // we got a cycle. Add it to the contour
 202       for ( Int_t j = 0; j <= inorder.GetLast(); ++j )
 203       {
 204         AliMUONSegment* s = static_cast<AliMUONSegment*>(inorder.UncheckedAt(j));
 205         polygon.SetVertex(++n,s->StartX(),s->StartY());
 206         all.Remove(s);
 207       }
 208
 209       all.Compress();
 210
 211       polygon.Close();
 212
 213       contour->Add(polygon);
 214
 215       if ( ! all.IsEmpty() )
 216       {
 217         i = 0;
 218         inorder.Clear();
 219         alreadyAdded.Set(all.GetLast()+1);
 220         alreadyAdded.Reset();
 221       }
 222       continue;
 223     }
 224
 225     for ( Int_t j = 0; j <= all.GetLast(); ++j)
 226     {
 227       if ( j != i && alreadyAdded[j] == 0 )
 228       {
 229         const AliMUONSegment* sj = static_cast<const AliMUONSegment*>(all.UncheckedAt(j));
 230         if ( AliMUONSegment::AreEqual(si->EndX(),sj->StartX()) && AliMUONSegment::AreEqual(si->EndY(),sj->StartY()))
 231         {
 232           i = j;
 233           break;
 234         }
 235       }
 236     }
 237   }
 238
 239   contour->AssertOrientation(kTRUE);
 240   return contour;
 241 }
 242
 243
 244 //_____________________________________________________________________________
 245 void
 246 AliMUONContourMaker::GetVerticalEdges(const TObjArray& polygons, TObjArray& polygonVerticalEdges) const
 247 {
 248   /// From an array of polygons, extract the list of vertical edges.
 249   /// Output array polygonVerticalEdges should be empty before calling.
 250
 251   AliCodeTimerAuto("",0);
 252
 253   for ( Int_t i = 0; i <= polygons.GetLast(); ++i )
 254   {
 255     const AliMUONPolygon* g = static_cast<const AliMUONPolygon*>(polygons.UncheckedAt(i));
 256     for ( Int_t j = 0; j < g->NumberOfVertices()-1; ++j )
 257     {
 258       if ( AliMUONSegment::AreEqual(g->X(j),g->X(j+1)) ) // segment is vertical
 259       {
 260         polygonVerticalEdges.Add(new AliMUONSegment(g->X(j),g->Y(j),g->X(j+1),g->Y(j+1)));
 261       }
 262     }
 263   }
 264 }
 265
 266
 267 //_____________________________________________________________________________
 268 void
 269 AliMUONContourMaker::GetYPositions(const TObjArray& polygonVerticalEdges,
 270                                    TArrayD& yPositions) const
 271 {
 272   /// Fill the array yPositions with the different y positions found in
 273   /// polygonVerticalEdges
 274
 275   AliCodeTimerAuto("",0);
 276
 277   Double_t* y = new Double_t[polygonVerticalEdges.GetSize()*2];
 278   Int_t n(0);
 279
 280   for ( Int_t i = 0; i < polygonVerticalEdges.GetLast(); ++i )
 281   {
 282     AliMUONSegment* s = static_cast<AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 283     y[n] = s->StartY();
 284     y[n+1] = s->EndY();
 285     n += 2;
 286   }
 287   Int_t* ix = new Int_t[n+1];
 288
 289   TMath::Sort(n,y,ix,kFALSE);
 290
 291   yPositions.Set(n+1);
 292
 293   Int_t u(0);
 294   Double_t x(FLT_MAX);
 295
 296   for ( Int_t i = 0; i < n; ++i )
 297   {
 298     if ( y[ix[i]] != x )
 299     {
 300       yPositions[u] = y[ix[i]];
 301       x = y[ix[i]];
 302       ++u;
 303     }
 304   }
 305
 306   yPositions.Set(u);
 307
 308   delete[] ix;
 309   delete[] y;
 310
 311 }
 312
 313 //_____________________________________________________________________________
 314 AliMUONContour*
 315 AliMUONContourMaker::MergeContour(const TObjArray& contours, const char* name) const
 316 {
 317   /// Merge all the polygons of all contours into a single contour
 318
 319   AliCodeTimerAuto("",0);
 320
 321   TObjArray polygons;
 322   polygons.SetOwner(kTRUE);
 323
 324   TIter next(&contours);
 325   AliMUONContour* contour;
 326   while ( ( contour = static_cast<AliMUONContour*>(next()) ) )
 327   {
 328     const TObjArray* contourPolygons = contour->Polygons();
 329     TIter nextPol(contourPolygons);
 330     AliMUONPolygon* pol;
 331     while ( ( pol = static_cast<AliMUONPolygon*>(nextPol()) ) )
 332     {
 333       polygons.Add(new AliMUONPolygon(*pol));
 334     }
 335   }
 336
 337   if ( polygons.IsEmpty() ) return 0x0;
 338
 339   contour = CreateContour(polygons,name);
 340
 341   return contour;
 342 }
 343
 344 //_____________________________________________________________________________
 345 void
 346 AliMUONContourMaker::SortPoints(const TObjArray& polygonVerticalEdges,
 347                                 TObjArray& sortedPoints) const
 348 {
 349   /// Sort the point of the vertical edges in ascending order, first on ordinate,
 350   /// then on abcissa, and put them in output vector sortedPoints.
 351   /// Output array sortedPoints should be empty before calling this method.
 352
 353   AliCodeTimerAuto("",0);
 354
 355   for ( Int_t i = 0; i <= polygonVerticalEdges.GetLast(); ++i )
 356   {
 357     const AliMUONSegment* e = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 358     sortedPoints.Add(new AliMUONPointWithRef(e->StartX(),e->StartY(),i));
 359     sortedPoints.Add(new AliMUONPointWithRef(e->EndX(),e->EndY(),i));
 360     // note that we keep track of the original edge, which is used
 361     // later on to deduce orientation of horizontal edges.
 362   }
 363
 364   sortedPoints.Sort();
 365 }
 366
 367 //_____________________________________________________________________________
 368 void
 369 AliMUONContourMaker::Sweep(const TObjArray& polygonVerticalEdges,
 370                            TObjArray& contourVerticalEdges) const
 371 {
 372   /// This is the meat of the algorithm of the contour merging...
 373
 374   AliCodeTimerAuto("",0);
 375
 376   TArrayD yPositions;
 377   GetYPositions(polygonVerticalEdges,yPositions);
 378
 379   AliMUONSegmentTree segmentTree(yPositions);
 380
 381   for ( Int_t i = 0; i <= polygonVerticalEdges.GetLast(); ++i )
 382   {
 383     const AliMUONSegment* edge = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i));
 384
 385     assert(edge!=0x0);
 386
 387     if ( edge->IsLeftEdge() )
 388     {
 389       segmentTree.Contribution(edge->Bottom(),edge->Top());
 390       segmentTree.InsertInterval(edge->Bottom(),edge->Top());
 391     }
 392     else
 393     {
 394       segmentTree.DeleteInterval(edge->Bottom(),edge->Top());
 395       segmentTree.Contribution(edge->Bottom(),edge->Top());
 396     }
 397
 398     AliMUONSegment e1(*edge);
 399
 400     if ( i < polygonVerticalEdges.GetLast() )
 401     {
 402       const AliMUONSegment* next = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(i+1));
 403       e1 = *next;
 404     }
 405
 406     if ( ( edge->IsLeftEdge() != e1.IsLeftEdge() ) ||
 407         ( !AliMUONSegment::AreEqual(edge->StartX(),e1.StartX() ) ) ||
 408         ( i == polygonVerticalEdges.GetLast() ) )
 409     {
 410       const TObjArray& stack = segmentTree.Stack();
 411
 412       double x = edge->StartX();
 413
 414       for ( Int_t j = 0; j <= stack.GetLast(); ++j )
 415       {
 416         AliMUONSegment* sj = static_cast<AliMUONSegment*>(stack.UncheckedAt(j));
 417         AliMUONSegment* s = new AliMUONSegment(x,sj->StartY(),x,sj->EndY());
 418
 419         if  (s->IsAPoint())
 420         {
 421           delete s;
 422           continue;
 423         }
 424
 425         if ( edge->IsLeftEdge() != s->IsLeftEdge() )
 426         {
 427           s->Set(x,sj->EndY(),x,sj->StartY());
 428         }
 429         contourVerticalEdges.Add(s);
 430       }
 431       segmentTree.ResetStack();
 432     }
 433   }
 434 }
 435
 436 //_____________________________________________________________________________
 437 void
 438 AliMUONContourMaker::VerticalToHorizontal(const TObjArray& polygonVerticalEdges,
 439                                           TObjArray& horizontalEdges) const
 440 {
 441   /// Deduce the set of horizontal edges from the vertical edges
 442   /// Output array horizontalEdges should be empty before calling this method
 443
 444   AliCodeTimerAuto("",0);
 445
 446   TObjArray points; // array of AliMUONPointWithRef
 447   points.SetOwner(kTRUE);
 448
 449   SortPoints(polygonVerticalEdges,points);
 450
 451   for ( Int_t k = 0; k < (points.GetLast()+1)/2; ++k )
 452   {
 453     const AliMUONPointWithRef* p1 = static_cast<AliMUONPointWithRef*>(points.UncheckedAt(k*2));
 454     const AliMUONPointWithRef* p2 = static_cast<AliMUONPointWithRef*>(points.UncheckedAt(k*2+1));
 455
 456     const AliMUONSegment* refEdge = static_cast<const AliMUONSegment*>(polygonVerticalEdges.UncheckedAt(p1->Ref()));
 457
 458     // (p1,p2) is the horizontal edge.
 459     // refEdge is used to deduce the orientation of (p1,p2)
 460
 461     if ( AliMUONSegment::AreEqual(p1->X(),refEdge->EndX()) && AliMUONSegment::AreEqual(p1->Y(),refEdge->EndY()) )
 462 //    if ( AreEqual(p1,refEdge->End()) )
 463     {
 464       horizontalEdges.Add(new AliMUONSegment(p1->X(),p1->Y(),p2->X(),p2->Y()));
 465     }
 466     else
 467     {
 468       horizontalEdges.Add(new AliMUONSegment(p2->X(),p2->Y(),p1->X(),p1->Y()));
 469     }
 470   }
 471 }
 472