MUON/AliMUONCluster.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 #include <Riostream.h>
  19 #include <TMath.h>
  20 #include <TObjArray.h>
  21 #include <TVirtualPad.h>
  22 #include <TVirtualX.h>
  23
  24 #include "AliMUONCluster.h"
  25
  26 #include "AliLog.h"
  27 #include "AliMUONPad.h"
  28
  29 /// \class AliMUONCluster
  30 ///
  31 /// A group of adjacent pads
  32 ///
  33 /// Besides holding an internal array of AliMUONPads, this object
  34 /// also computes some global characteristics for that pad sets.
  35 ///
  36 /// \author Laurent Aphecetche
  37 ///
  38
  39 /// \cond CLASSIMP
  40 ClassImp(AliMUONCluster)
  41 /// \endcond
  42
  43 namespace
  44 {
  45   //___________________________________________________________________________
  46   Bool_t
  47   ShouldUsePad(const AliMUONPad& pad,
  48                Int_t cathode, Int_t statusMask, Bool_t matchMask)
  49   {
  50
  51       // FIXME : we should only use >=0 status, so we can fully
  52       // use masking possibility ?
  53     if ( pad.Status() < 0 ) return kFALSE;
  54
  55     if ( pad.Cathode() == cathode && pad.IsReal() && !pad.IsSaturated() )
  56     {
  57       Bool_t test = ( ( pad.Status() & statusMask ) != 0 );
  58       if ( !statusMask )
  59       {
  60         test = ( pad.Status() == 0 );
  61       }
  62       if ( ( test && matchMask ) || ( !test && !matchMask ) )
  63       {
  64         return kTRUE;
  65       }
  66     }
  67     return kFALSE;
  68   }
  69
  70   //___________________________________________________________________________
  71   Int_t Unique(Int_t n, Double_t* array, Double_t precision)
  72   {
  73     /// Return the number of *different* elements in array
  74     /// where different is up to precision
  75     /// Note that we assume that n is >= 1
  76
  77     Int_t count(1);
  78
  79     Int_t* index = new Int_t[n];
  80
  81     TMath::Sort(n,array,index);
  82
  83     for ( Int_t i = 1; i < n; ++i )
  84     {
  85       if ( array[index[i]] - array[index[i-1]] < -precision ) ++count;
  86     }
  87
  88     delete[] index;
  89
  90     return count;
  91   }
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliMUONCluster::AliMUONCluster()
  96 : TObject(),
  97 fPads(0x0),
  98 fHasPosition(kFALSE),
  99 fPosition(1E9,1E9),
 100 fPositionError(1E9,1E9),
 101 fHasCharge(kFALSE),
 102 fChi2(0),
 103 fIsSorted(kFALSE)
 104 {
 105   /// ctor
 106   fMultiplicity[0]=fMultiplicity[1]=0;
 107   fRawCharge[0]=fRawCharge[1]=0;
 108   fCharge[0]=fCharge[1]=0;
 109   fIsSaturated[0]=fIsSaturated[1]=kFALSE;
 110 }
 111
 112 //_____________________________________________________________________________
 113 AliMUONCluster::AliMUONCluster(const AliMUONCluster& src)
 114 : TObject(src),
 115 fPads(0x0),
 116 fHasPosition(kFALSE),
 117 fPosition(1E9,1E9),
 118 fPositionError(1E9,1E9),
 119 fHasCharge(kFALSE),
 120 fChi2(0),
 121 fIsSorted(kFALSE)
 122 {
 123   /// copy ctor
 124   src.Copy(*this);
 125 }
 126
 127 //_____________________________________________________________________________
 128 AliMUONCluster&
 129 AliMUONCluster::operator=(const AliMUONCluster& src)
 130 {
 131   /// assignement operator
 132   AliMUONCluster c(src);
 133   c.Copy(*this);
 134   return *this;
 135 }
 136
 137 //_____________________________________________________________________________
 138 AliMUONCluster::~AliMUONCluster()
 139 {
 140   /// dtor : note that we're owner of our pads
 141   delete fPads;
 142 }
 143
 144 //_____________________________________________________________________________
 145 void
 146 AliMUONCluster::AddPad(const AliMUONPad& pad)
 147 {
 148   /// Add a pad to our pad array, and update some internal information
 149   /// accordingly.
 150   /// If pad array was sorted prior to this call, we re-sort it after
 151   /// actual addition.
 152
 153   if (!fPads)
 154   {
 155     fPads = new TObjArray(10);
 156     fPads->SetOwner(kTRUE);
 157   }
 158   AliMUONPad* p = new AliMUONPad(pad);
 159   fPads->AddLast(p);
 160   p->SetClusterId(GetUniqueID());
 161   Int_t cathode = p->Cathode();
 162   ++(fMultiplicity[cathode]);
 163   fRawCharge[cathode] += p->Charge();
 164   if ( p->IsSaturated() )
 165   {
 166     fIsSaturated[p->Cathode()]=kTRUE;
 167   }
 168   //AZ if ( fIsSorted ) Sort();
 169 }
 170
 171 //_____________________________________________________________________________
 172 AliMpArea
 173 AliMUONCluster::Area() const
 174 {
 175   /// Return the geometrical area covered by this cluster
 176
 177   // Start by finding the (x,y) limits of this cluster
 178   TVector2 lowerLeft(1E9,1E9);
 179   TVector2 upperRight(-1E9,-1E9);
 180
 181   for ( Int_t i = 0; i < Multiplicity(); ++i )
 182   {
 183     AliMUONPad* pad = Pad(i);
 184     TVector2 ll = pad->Position() - pad->Dimensions();
 185     TVector2 ur = pad->Position() + pad->Dimensions();
 186     lowerLeft.Set( TMath::Min(ll.X(),lowerLeft.X()),
 187                    TMath::Min(ll.Y(),lowerLeft.Y()) );
 188     upperRight.Set( TMath::Max(ur.X(),upperRight.X()),
 189                     TMath::Max(ur.Y(),upperRight.Y()) );
 190   }
 191
 192   // then construct the area from those limits
 193   return AliMpArea((lowerLeft+upperRight)/2,
 194                    (upperRight-lowerLeft)/2);
 195 }
 196
 197 //_____________________________________________________________________________
 198 void
 199 AliMUONCluster::Copy(TObject& obj) const
 200 {
 201   ///
 202   /// Copy this cluster to (cluster&)obj
 203   ///
 204   TObject::Copy(obj);
 205   AliMUONCluster& dest = static_cast<AliMUONCluster&>(obj);
 206   dest.fPads = 0x0;
 207   if ( fPads )
 208   {
 209     dest.fPads = static_cast<TObjArray*>(fPads->Clone());
 210   }
 211   dest.fHasPosition = fHasPosition;
 212   dest.fPosition = fPosition;
 213   dest.fPositionError = fPositionError;
 214   dest.fHasCharge = fHasCharge;
 215   dest.fIsSorted = fIsSorted;
 216   dest.fChi2 = fChi2;
 217   for ( Int_t i = 0; i < 2; ++i )
 218   {
 219     dest.fRawCharge[i] = fRawCharge[i];
 220     dest.fCharge[i] = fCharge[i];
 221     dest.fMultiplicity[i] = fMultiplicity[i];
 222     dest.fIsSaturated[i] = fIsSaturated[i];
 223   }
 224 }
 225
 226 //_____________________________________________________________________________
 227 Float_t
 228 AliMUONCluster::Charge() const
 229 {
 230   /// Return the average charge over both cathodes
 231   return (Charge(0)+Charge(1))/2.0;
 232 }
 233
 234 //_____________________________________________________________________________
 235 Float_t
 236 AliMUONCluster::Charge(Int_t cathode) const
 237 {
 238   /// Returns the charge of a given cathode
 239   if ( !fHasCharge ) return RawCharge(cathode);
 240
 241   if ( cathode == 0 || cathode == 1 )
 242   {
 243     return fCharge[cathode];
 244   }
 245   return 0;
 246 }
 247
 248 //_____________________________________________________________________________
 249 Float_t
 250 AliMUONCluster::ChargeAsymmetry() const
 251 {
 252   /// Returns the charge asymmetry
 253   if ( Charge() > 0 )
 254   {
 255     return TMath::Abs(Charge(0)-Charge(1))/Charge();
 256   }
 257   return 0;
 258 }
 259
 260 //_____________________________________________________________________________
 261 TVector2
 262 AliMUONCluster::MinPadDimensions(Int_t statusMask, Bool_t matchMask) const
 263 {
 264   /// Returns the minimum pad dimensions (half sizes), only considering
 265   /// pads matching (or not, depending matchMask) a given mask
 266
 267   TVector2 cath0(MinPadDimensions(0,statusMask,matchMask));
 268   TVector2 cath1(MinPadDimensions(1,statusMask,matchMask));
 269
 270   return TVector2( TMath::Min(cath0.X(),cath1.X()),
 271                    TMath::Min(cath0.Y(),cath1.Y()) );
 272 }
 273
 274 //_____________________________________________________________________________
 275 TVector2
 276 AliMUONCluster::MinPadDimensions(Int_t cathode,
 277                                  Int_t statusMask, Bool_t matchMask) const
 278 {
 279   /// Returns the minimum pad dimensions (half sizes), only considering
 280   /// pads matching (or not, depending matchMask) a given mask, within a
 281   /// given cathode
 282
 283   Double_t xmin(1E9);
 284   Double_t ymin(1E9);
 285
 286   for ( Int_t i = 0; i < Multiplicity(); ++i )
 287   {
 288     AliMUONPad* pad = Pad(i);
 289     if ( ShouldUsePad(*pad,cathode,statusMask,matchMask) )
 290     {
 291       xmin = TMath::Min(xmin,pad->DX());
 292       ymin = TMath::Min(ymin,pad->DY());
 293     }
 294   }
 295   return TVector2(xmin,ymin);
 296 }
 297
 298 //_____________________________________________________________________________
 299 Int_t
 300 AliMUONCluster::Multiplicity() const
 301 {
 302   /// Returns the total number of pads in this cluster
 303   return Multiplicity(0)+Multiplicity(1);
 304 }
 305
 306 //_____________________________________________________________________________
 307 Int_t
 308 AliMUONCluster::Multiplicity(Int_t cathode) const
 309 {
 310   /// Returns the number of pads in this cluster, in the given cathode
 311   if ( cathode == 0 || cathode == 1 )
 312   {
 313     return fMultiplicity[cathode];
 314   }
 315   return 0;
 316 }
 317
 318 //_____________________________________________________________________________
 319 AliMpIntPair
 320 AliMUONCluster::NofPads(Int_t statusMask, Bool_t matchMask) const
 321 {
 322   /// Number of pads satisfying (or not, depending matchMask) a
 323   /// given mask
 324
 325   Int_t nx, ny;
 326
 327   TVector2 dim0(MinPadDimensions(0,statusMask,matchMask));
 328   TVector2 dim1(MinPadDimensions(1,statusMask,matchMask));
 329
 330   AliMpIntPair npad0(NofPads(0,statusMask,matchMask));
 331   AliMpIntPair npad1(NofPads(1,statusMask,matchMask));
 332
 333   if ( TMath::Abs( (dim0-dim1).X() ) < 1E-3 )
 334   {
 335     nx = TMath::Max( npad0.GetFirst(), npad1.GetFirst() );
 336   }
 337   else
 338   {
 339     nx = dim0.X() < dim1.X() ? npad0.GetFirst() : npad1.GetFirst();
 340   }
 341
 342   if ( TMath::Abs( (dim0-dim1).Y() ) < 1E-3 )
 343   {
 344     ny = TMath::Max( npad0.GetSecond(), npad1.GetSecond() );
 345   }
 346   else
 347   {
 348     ny = dim0.Y() < dim1.Y() ? npad0.GetSecond() : npad1.GetSecond();
 349   }
 350
 351   return AliMpIntPair(nx,ny);
 352 }
 353
 354 //_____________________________________________________________________________
 355 AliMpIntPair
 356 AliMUONCluster::NofPads(Int_t cathode,
 357                         Int_t statusMask, Bool_t matchMask) const
 358 {
 359   /// Number of pads of a given cathode, satisfying (or not,
 360   /// depending matchMask) a given mask
 361
 362   Int_t n = Multiplicity(cathode);
 363   if (!n)
 364   {
 365     return AliMpIntPair(0,0);
 366   }
 367   Double_t* x = new Double_t[n];
 368   Double_t* y = new Double_t[n];
 369   Int_t np(0);
 370
 371   for ( Int_t i = 0; i < Multiplicity(); ++i )
 372   {
 373     AliMUONPad* pad = Pad(i);
 374     if ( ShouldUsePad(*pad,cathode,statusMask,matchMask) )
 375     {
 376       x[np] = pad->X();
 377       y[np] = pad->Y();
 378       ++np;
 379     }
 380   }
 381
 382   Int_t cx = Unique(np,x,0.01);
 383   Int_t cy = Unique(np,y,0.01);
 384
 385   delete[] x;
 386   delete[] y;
 387
 388   return AliMpIntPair(cx,cy);
 389 }
 390
 391 //_____________________________________________________________________________
 392 AliMUONPad*
 393 AliMUONCluster::Pad(Int_t index) const
 394 {
 395   /// Returns the index-th pad
 396
 397   if (!fPads) return 0x0;
 398   if ( index < fPads->GetLast()+1 )
 399   {
 400     return static_cast<AliMUONPad*>(fPads->At(index));
 401   }
 402   else
 403   {
 404     AliError(Form("Requesting index %d out of bounds (%d)",index,fPads->GetLast()));
 405   }
 406   return 0x0;
 407 }
 408
 409
 410 //_____________________________________________________________________________
 411 void
 412 AliMUONCluster::Paint(Option_t*)
 413 {
 414   /// Paint this cluster
 415   if (!Multiplicity()) return;
 416
 417   AliMpArea area(Area());
 418
 419   gPad->Range(area.LeftBorder(),area.DownBorder(),area.RightBorder(),area.UpBorder());
 420
 421   gVirtualX->SetFillStyle(0);
 422
 423   gVirtualX->SetLineColor(2);
 424   gVirtualX->SetLineWidth(4);
 425   for ( Int_t i = 0; i < Multiplicity(); ++i)
 426   {
 427     AliMUONPad* pad = Pad(i);
 428     if ( pad->Cathode() == 0 ) pad->Paint();
 429   }
 430
 431   gVirtualX->SetLineColor(4);
 432   gVirtualX->SetLineWidth(2);
 433   for ( Int_t i = 0; i < Multiplicity(); ++i)
 434   {
 435     AliMUONPad* pad = Pad(i);
 436     if ( pad->Cathode() == 1 ) pad->Paint();
 437   }
 438
 439 }
 440
 441 //_____________________________________________________________________________
 442 void
 443 AliMUONCluster::Print(Option_t* opt) const
 444 {
 445   /// printout
 446   cout << "Cluster Id " << GetUniqueID() << " npads=" << Multiplicity()
 447   << "(" << Multiplicity(0) << "," << Multiplicity(1) << ") RawCharge="
 448   << RawCharge() << " (" << RawCharge(0) << "," << RawCharge(1)
 449   << ") Charge=(" << Charge(0) << "," << Charge(1) <<")";
 450   if ( HasPosition() )
 451   {
 452     cout << " (x,y)=(" << Position().X() << "," << Position().Y() << ")";
 453     cout << " (errX,errY)=(" << PositionError().X() << "," << PositionError().Y() << ")";
 454   }
 455   cout << " " << Area() << endl;
 456   if (fPads)
 457   {
 458     fPads->Print("",opt);
 459   }
 460 }
 461
 462 //_____________________________________________________________________________
 463 void
 464 AliMUONCluster::Sort()
 465 {
 466   /// Sort the pad array
 467   fPads->Sort();
 468   fIsSorted = kTRUE;
 469 }
 470
 471 //_____________________________________________________________________________
 472 void
 473 AliMUONCluster::RemovePad(AliMUONPad* pad)
 474 {
 475   /// Remove a pad.
 476   /// As a consequence, some internal information must be updated
 477
 478   fPads->Remove(pad);
 479   fPads->Compress();
 480   // update cluster's data
 481   fIsSaturated[0]=fIsSaturated[1]=kFALSE;
 482   fMultiplicity[0]=fMultiplicity[1]=0;
 483   fRawCharge[0]=fRawCharge[1]=0;
 484   for ( Int_t i = 0; i <= fPads->GetLast(); ++i )
 485   {
 486     AliMUONPad* p = Pad(i);
 487     if ( p->IsSaturated() )
 488     {
 489       fIsSaturated[p->Cathode()] = kTRUE;
 490     }
 491     ++fMultiplicity[p->Cathode()];
 492     fRawCharge[p->Cathode()] += p->Charge();
 493   }
 494   if (fIsSorted) Sort();
 495 }
 496
 497 //_____________________________________________________________________________
 498 Float_t
 499 AliMUONCluster::RawCharge() const
 500 {
 501   /// Returns the raw average charge
 502   return (RawCharge(0)+RawCharge(1))/2.0;
 503 }
 504
 505 //_____________________________________________________________________________
 506 Float_t
 507 AliMUONCluster::RawCharge(Int_t cathode) const
 508 {
 509   /// Returns the average charge of a given cathode
 510   if ( cathode == 0 || cathode == 1 )
 511   {
 512     return fRawCharge[cathode];
 513   }
 514   return 0;
 515 }
 516
 517 //_____________________________________________________________________________
 518 Float_t
 519 AliMUONCluster::RawChargeAsymmetry() const
 520 {
 521   /// Returns the raw charge asymmetry
 522    if ( RawCharge() > 0 )
 523    {
 524      return TMath::Abs(RawCharge(0)-RawCharge(1))/RawCharge();
 525    }
 526   return 0;
 527 }