MUON/AliMUONVTrackReconstructor.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 ///
  20 /// \class AliMUONVTrackReconstructor
  21 /// Virtual MUON track reconstructor in ALICE (class renamed from AliMUONEventReconstructor)
  22 ///
  23 /// This class contains as data:
  24 /// * a pointer to the array of hits to be reconstructed (the event)
  25 /// * a pointer to the array of segments made with these hits inside each station
  26 /// * a pointer to the array of reconstructed tracks
  27 ///
  28 /// It contains as methods, among others:
  29 /// * EventReconstruct to build the muon tracks
  30 /// * EventReconstructTrigger to build the trigger tracks
  31 ///
  32 ///  \author Philippe Pillot
  33 ///
  34 ////////////////////////////////////
  35
  36 #include "AliMUONVTrackReconstructor.h"
  37
  38 #include "AliMpDEManager.h"
  39 #include "AliMUONConstants.h"
  40 #include "AliMUONHitForRec.h"
  41 #include "AliMUONObjectPair.h"
  42 #include "AliMUONTriggerTrack.h"
  43 #include "AliMUONTriggerCircuit.h"
  44 #include "AliMUONLocalTrigger.h"
  45 #include "AliMUONGlobalTrigger.h"
  46 #include "AliMUONTrack.h"
  47 #include "AliMUONTrackParam.h"
  48 #include "AliMUONTrackExtrap.h"
  49 #include "AliMUONTrackHitPattern.h"
  50 #include "AliMUONVTrackStore.h"
  51 #include "AliMUONVClusterStore.h"
  52 #include "AliMUONRawCluster.h"
  53 #include "AliLog.h"
  54 #include "AliTracker.h"
  55 #include "AliMUONVTriggerStore.h"
  56 #include "AliMUONVTriggerTrackStore.h"
  57
  58 #include <Riostream.h>
  59 #include <TClonesArray.h>
  60 #include <TMath.h>
  61
  62 /// \cond CLASSIMP
  63 ClassImp(AliMUONVTrackReconstructor) // Class implementation in ROOT context
  64 /// \endcond
  65
  66 //************* Defaults parameters for reconstruction
  67 const Double_t AliMUONVTrackReconstructor::fgkDefaultMinBendingMomentum = 3.0;
  68 const Double_t AliMUONVTrackReconstructor::fgkDefaultMaxBendingMomentum = 3000.0;
  69 const Double_t AliMUONVTrackReconstructor::fgkDefaultBendingResolution = 0.01;
  70 const Double_t AliMUONVTrackReconstructor::fgkDefaultNonBendingResolution = 0.144;
  71 const Double_t AliMUONVTrackReconstructor::fgkDefaultBendingVertexDispersion = 10.;
  72 const Double_t AliMUONVTrackReconstructor::fgkDefaultNonBendingVertexDispersion = 10.;
  73 const Double_t AliMUONVTrackReconstructor::fgkDefaultMaxNormChi2MatchTrigger = 16.0;
  74
  75 //__________________________________________________________________________
  76 AliMUONVTrackReconstructor::AliMUONVTrackReconstructor()
  77   : TObject(),
  78     fMinBendingMomentum(fgkDefaultMinBendingMomentum),
  79     fMaxBendingMomentum(fgkDefaultMaxBendingMomentum),
  80     fBendingResolution(fgkDefaultBendingResolution),
  81     fNonBendingResolution(fgkDefaultNonBendingResolution),
  82     fBendingVertexDispersion(fgkDefaultBendingVertexDispersion),
  83     fNonBendingVertexDispersion(fgkDefaultNonBendingVertexDispersion),
  84     fMaxNormChi2MatchTrigger(fgkDefaultMaxNormChi2MatchTrigger),
  85     fHitsForRecPtr(0x0),
  86     fNHitsForRec(0),
  87     fNHitsForRecPerChamber(0x0),
  88     fIndexOfFirstHitForRecPerChamber(0x0),
  89     fRecTracksPtr(0x0),
  90     fNRecTracks(0)
  91 {
  92   /// Constructor for class AliMUONVTrackReconstructor
  93   fNHitsForRecPerChamber = new Int_t[AliMUONConstants::NTrackingCh()];
  94   fIndexOfFirstHitForRecPerChamber = new Int_t[AliMUONConstants::NTrackingCh()];
  95
  96   // Memory allocation for the TClonesArray of hits for reconstruction
  97   // Is 10000 the right size ????
  98   fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
  99   fHitsForRecPtr->SetOwner(kTRUE);
 100
 101   // set the magnetic field for track extrapolations
 102   const AliMagF* kField = AliTracker::GetFieldMap();
 103   if (!kField) AliFatal("No field available");
 104   AliMUONTrackExtrap::SetField(kField);
 105 }
 106
 107   //__________________________________________________________________________
 108 AliMUONVTrackReconstructor::~AliMUONVTrackReconstructor(void)
 109 {
 110   /// Destructor for class AliMUONVTrackReconstructor
 111   delete [] fNHitsForRecPerChamber;
 112   delete [] fIndexOfFirstHitForRecPerChamber;
 113   delete fHitsForRecPtr;
 114 }
 115
 116 //__________________________________________________________________________
 117 void
 118 AliMUONVTrackReconstructor::AddHitsForRecFromRawClusters(const AliMUONVClusterStore& clusterStore)
 119 {
 120   /// Build internal array of hit for rec from clusterStore
 121
 122   TIter next(clusterStore.CreateIterator());
 123   AliMUONRawCluster* clus(0x0);
 124   Int_t iclus(0);
 125
 126   while ( ( clus = static_cast<AliMUONRawCluster*>(next()) ) )
 127   {
 128     // new AliMUONHitForRec from raw cluster
 129     // and increment number of AliMUONHitForRec's (total and in chamber)
 130     AliMUONHitForRec* hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
 131     fNHitsForRec++;
 132     // more information into HitForRec
 133     hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
 134     hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
 135     //  original raw cluster
 136     Int_t ch = AliMpDEManager::GetChamberId(clus->GetDetElemId());
 137     hitForRec->SetChamberNumber(ch);
 138     hitForRec->SetHitNumber(iclus);
 139     // Z coordinate of the raw cluster (cm)
 140     hitForRec->SetZ(clus->GetZ(0));
 141     if (AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 3) {
 142       cout << "Chamber " << ch <<" raw cluster  " << iclus << " : " << endl;
 143       clus->Print("full");
 144       cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 145       hitForRec->Print("full");
 146     }
 147     ++iclus;
 148   } // end of chamber loop
 149
 150   SortHitsForRecWithIncreasingChamber();
 151
 152   AliDebug(1,"End of AddHitsForRecFromRawClusters");
 153
 154   if (AliLog::GetGlobalDebugLevel() > 0)
 155   {
 156     AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
 157     for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++)
 158     {
 159       AliDebug(1, Form("Chamber(0...): %d",ch));
 160       AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
 161       AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
 162       for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
 163            hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
 164            hit++) {
 165         AliDebug(1, Form("HitForRec index(0...): %d",hit));
 166         ((*fHitsForRecPtr)[hit])->Dump();
 167       }
 168     }
 169   }
 170 }
 171
 172 //__________________________________________________________________________
 173 void AliMUONVTrackReconstructor::EventReconstruct(const AliMUONVClusterStore& clusterStore,
 174                                                   AliMUONVTrackStore& trackStore)
 175 {
 176   /// To reconstruct one event
 177   AliDebug(1,"");
 178
 179   ResetTracks(); //AZ
 180   ResetHitsForRec(); //AZ
 181   AddHitsForRecFromRawClusters(clusterStore);
 182   MakeTracks();
 183
 184   // Add tracks to MUON data container
 185   for (Int_t i=0; i<fNRecTracks; ++i)
 186   {
 187     AliMUONTrack * track = (AliMUONTrack*) fRecTracksPtr->At(i);
 188     trackStore.Add(*track);
 189   }
 190 }
 191
 192 //__________________________________________________________________________
 193 void AliMUONVTrackReconstructor::ResetTracks(void)
 194 {
 195   /// To reset the TClonesArray of reconstructed tracks
 196   if (fRecTracksPtr) fRecTracksPtr->Clear("C");
 197   // Delete in order that the Track destructors are called,
 198   // hence the space for the TClonesArray of pointers to TrackHit's is freed
 199   fNRecTracks = 0;
 200   return;
 201 }
 202
 203   //__________________________________________________________________________
 204 void AliMUONVTrackReconstructor::ResetHitsForRec(void)
 205 {
 206   /// To reset the array and the number of HitsForRec,
 207   /// and also the number of HitsForRec
 208   /// and the index of the first HitForRec per chamber
 209   if (fHitsForRecPtr) fHitsForRecPtr->Clear("C");
 210   fNHitsForRec = 0;
 211   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++)
 212     fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
 213   return;
 214 }
 215
 216   //__________________________________________________________________________
 217 void AliMUONVTrackReconstructor::SortHitsForRecWithIncreasingChamber()
 218 {
 219   /// Sort HitsForRec's in increasing order with respect to chamber number.
 220   /// Uses the function "Compare".
 221   /// Update the information for HitsForRec per chamber too.
 222   Int_t ch, nhits, prevch;
 223   fHitsForRecPtr->Sort();
 224   for (ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 225     fNHitsForRecPerChamber[ch] = 0;
 226     fIndexOfFirstHitForRecPerChamber[ch] = 0;
 227   }
 228   prevch = 0; // previous chamber
 229   nhits = 0; // number of hits in current chamber
 230   // Loop over HitsForRec
 231   for (Int_t hit = 0; hit < fNHitsForRec; hit++) {
 232     // chamber number (0...)
 233     ch = ((AliMUONHitForRec*)  ((*fHitsForRecPtr)[hit]))->GetChamberNumber();
 234     // increment number of hits in current chamber
 235     (fNHitsForRecPerChamber[ch])++;
 236     // update index of first HitForRec in current chamber
 237     // if chamber number different from previous one
 238     if (ch != prevch) {
 239       fIndexOfFirstHitForRecPerChamber[ch] = hit;
 240       prevch = ch;
 241     }
 242   }
 243   return;
 244 }
 245
 246   //__________________________________________________________________________
 247 TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsInStation(Int_t station)
 248 {
 249   /// To make the list of segments in station(0..) "Station" from the list of hits to be reconstructed.
 250   /// Return a new TClonesArray of segments.
 251   /// It is the responsibility of the user to delete it afterward.
 252   AliDebug(1,Form("Enter MakeSegmentsPerStation (0...) %d",station));
 253
 254   AliMUONHitForRec *hit1Ptr, *hit2Ptr;
 255   AliMUONObjectPair *segment;
 256   Double_t bendingSlope = 0, impactParam = 0., bendingMomentum = 0.; // to avoid compilation warning
 257                                                                      // first and second chambers (0...) in the station
 258   Int_t ch1 = 2 * station;
 259   Int_t ch2 = ch1 + 1;
 260
 261   // list of segments
 262   TClonesArray *segments = new TClonesArray("AliMUONObjectPair", fNHitsForRecPerChamber[ch2]);
 263   segments->SetOwner(kTRUE);
 264
 265   // Loop over HitForRec's in the first chamber of the station
 266   for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
 267        hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
 268        hit1++)
 269   {
 270     // pointer to the HitForRec
 271     hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
 272     // Loop over HitsForRec's in the second chamber of the station
 273     for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
 274          hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
 275          hit2++)
 276     {
 277       // pointer to the HitForRec
 278       hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]);
 279       if ( hit1Ptr->GetZ() - hit2Ptr->GetZ() != 0. )
 280       {
 281         // bending slope
 282         bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) / (hit1Ptr->GetZ() - hit2Ptr->GetZ());
 283         // impact parameter
 284         impactParam = hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope;
 285         // absolute value of bending momentum
 286         bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(impactParam));
 287       } else
 288       {
 289         AliWarning("hit1Ptr->GetZ() = hit2Ptr->GetZ(): no segment created");
 290         continue;
 291       }
 292       // check for bending momentum within tolerances
 293       if ((bendingMomentum < fMaxBendingMomentum) && (bendingMomentum > fMinBendingMomentum))
 294       {
 295         // make new segment
 296         segment = new ((*segments)[segments->GetLast()+1]) AliMUONObjectPair(hit1Ptr, hit2Ptr, kFALSE, kFALSE);
 297         if (AliLog::GetGlobalDebugLevel() > 1) {
 298           cout << "segmentIndex(0...): " << segments->GetLast() << endl;
 299           segment->Dump();
 300           cout << "HitForRec in first chamber" << endl;
 301           hit1Ptr->Dump();
 302           cout << "HitForRec in second chamber" << endl;
 303           hit2Ptr->Dump();
 304         }
 305       }
 306     } //for (Int_t hit2
 307   } // for (Int_t hit1...
 308   AliDebug(1,Form("Station: %d  NSegments:  %d ", station, segments->GetEntriesFast()));
 309   return segments;
 310 }
 311
 312 //__________________________________________________________________________
 313 void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(AliMUONVTrackStore& trackStore,
 314                                                            const AliMUONVTriggerTrackStore& triggerTrackStore,
 315                                                            const AliMUONVTriggerStore& triggerStore,
 316                                                            const AliMUONTrackHitPattern& trackHitPattern)
 317 {
 318   /// Try to match track from tracking system with trigger track
 319   static const Double_t kDistSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY
 320
 321   Int_t matchTrigger;
 322   Int_t loTrgNum(-1);
 323   Double_t distTriggerTrack[3];
 324   Double_t xTrack, yTrack, ySlopeTrack, chi2MatchTrigger, minChi2MatchTrigger, chi2;
 325
 326   TIter itTrack(trackStore.CreateIterator());
 327   AliMUONTrack* track;
 328
 329   while ( ( track = static_cast<AliMUONTrack*>(itTrack()) ) )
 330   {
 331     matchTrigger = 0;
 332     chi2MatchTrigger = 0.;
 333     loTrgNum = -1;
 334     Int_t doubleMatch=-1; // Check if track matches 2 trigger tracks
 335     Double_t doubleChi2 = -1.;
 336
 337     AliMUONTrackParam trackParam(*((AliMUONTrackParam*) (track->GetTrackParamAtHit()->Last())));
 338     AliMUONTrackExtrap::ExtrapToZ(&trackParam, AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber
 339
 340     xTrack = trackParam.GetNonBendingCoor();
 341     yTrack = trackParam.GetBendingCoor();
 342     ySlopeTrack = trackParam.GetBendingSlope();
 343     minChi2MatchTrigger = 999.;
 344
 345     AliMUONTriggerTrack *triggerTrack;
 346     TIter itTriggerTrack(triggerTrackStore.CreateIterator());
 347     while ( ( triggerTrack = static_cast<AliMUONTriggerTrack*>(itTriggerTrack() ) ) )
 348     {
 349       distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/kDistSigma[0];
 350       distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/kDistSigma[1];
 351       distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/kDistSigma[2];
 352       chi2 = 0.;
 353       for (Int_t iVar = 0; iVar < 3; iVar++) chi2 += distTriggerTrack[iVar]*distTriggerTrack[iVar];
 354       chi2 /= 3.; // Normalized Chi2: 3 degrees of freedom (X,Y,slopeY)
 355       if (chi2 < fMaxNormChi2MatchTrigger)
 356       {
 357         Bool_t isDoubleTrack = (TMath::Abs(chi2 - minChi2MatchTrigger)<1.);
 358         if (chi2 < minChi2MatchTrigger && chi2 < fMaxNormChi2MatchTrigger)
 359         {
 360           if(isDoubleTrack)
 361           {
 362             doubleMatch = loTrgNum;
 363             doubleChi2 = chi2MatchTrigger;
 364           }
 365           minChi2MatchTrigger = chi2;
 366           chi2MatchTrigger = chi2;
 367           loTrgNum = triggerTrack->GetLoTrgNum();
 368           AliMUONLocalTrigger* locTrg = triggerStore.FindLocal(loTrgNum);
 369           matchTrigger=1;
 370           if(locTrg->LoLpt()>0)matchTrigger=2;
 371           if(locTrg->LoHpt()>0)matchTrigger=3;
 372         }
 373         else if(isDoubleTrack)
 374         {
 375           doubleMatch = triggerTrack->GetLoTrgNum();
 376           doubleChi2 = chi2;
 377         }
 378       }
 379     }
 380     if(doubleMatch>=0)
 381     { // If two trigger tracks match, select the one passing more trigger cuts
 382       AliDebug(1, Form("Two candidates found: %i and %i",loTrgNum,doubleMatch));
 383       AliMUONLocalTrigger* locTrg1 = triggerStore.FindLocal(doubleMatch);
 384       if((locTrg1->LoLpt()>0 && matchTrigger<2) || (locTrg1->LoHpt() && matchTrigger<3))
 385       {
 386         if(locTrg1->LoHpt()>0)matchTrigger=3;
 387         else matchTrigger=2;
 388         loTrgNum = doubleMatch;
 389         chi2MatchTrigger=doubleChi2;
 390       }
 391     }
 392
 393     track->SetMatchTrigger(matchTrigger);
 394     track->SetLoTrgNum(loTrgNum);
 395     track->SetChi2MatchTrigger(chi2MatchTrigger);
 396
 397     AliMUONLocalTrigger* locTrg = static_cast<AliMUONLocalTrigger*>(triggerStore.FindLocal(loTrgNum));
 398
 399     if (locTrg)
 400     {
 401       track->SetLocalTrigger(locTrg->LoCircuit(),
 402                              locTrg->LoStripX(),
 403                              locTrg->LoStripY(),
 404                              locTrg->LoDev(),
 405                              locTrg->LoLpt(),
 406                              locTrg->LoHpt());
 407     }
 408   }
 409
 410   trackHitPattern.GetHitPattern(trackStore,triggerStore);
 411 }
 412
 413 //__________________________________________________________________________
 414 void
 415 AliMUONVTrackReconstructor::EventReconstructTrigger(const AliMUONTriggerCircuit& circuit,
 416                                                     const AliMUONVTriggerStore& triggerStore,
 417                                                     AliMUONVTriggerTrackStore& triggerTrackStore)
 418 {
 419   /// To make the trigger tracks from Local Trigger
 420   AliDebug(1, "");
 421
 422   AliMUONGlobalTrigger* globalTrigger = triggerStore.Global();
 423
 424   UChar_t gloTrigPat = 0;
 425
 426   if (globalTrigger)
 427   {
 428     gloTrigPat = globalTrigger->GetGlobalResponse();
 429   }
 430
 431   TIter next(triggerStore.CreateIterator());
 432   AliMUONLocalTrigger* locTrg(0x0);
 433
 434   Float_t z11 = AliMUONConstants::DefaultChamberZ(10);
 435   Float_t z21 = AliMUONConstants::DefaultChamberZ(12);
 436
 437   AliMUONTriggerTrack triggerTrack;
 438
 439   while ( ( locTrg = static_cast<AliMUONLocalTrigger*>(next()) ) )
 440   {
 441     Bool_t xTrig=kFALSE;
 442     Bool_t yTrig=kFALSE;
 443
 444     Int_t localBoardId = locTrg->LoCircuit();
 445     if ( locTrg->LoSdev()==1 && locTrg->LoDev()==0 &&
 446          locTrg->LoStripX()==0) xTrig=kFALSE; // no trigger in X
 447     else xTrig=kTRUE;                         // trigger in X
 448     if (locTrg->LoTrigY()==1 &&
 449         locTrg->LoStripY()==15 ) yTrig = kFALSE; // no trigger in Y
 450     else yTrig = kTRUE;                          // trigger in Y
 451
 452     if (xTrig && yTrig)
 453     { // make Trigger Track if trigger in X and Y
 454
 455       Float_t y11 = circuit.GetY11Pos(localBoardId, locTrg->LoStripX());
 456       // need first to convert deviation to [0-30]
 457       // (see AliMUONLocalTriggerBoard::LocalTrigger)
 458       Int_t deviation = locTrg->LoDev();
 459       Int_t sign = 0;
 460       if ( !locTrg->LoSdev() &&  deviation ) sign=-1;
 461       if ( !locTrg->LoSdev() && !deviation ) sign= 0;
 462       if (  locTrg->LoSdev() == 1 )          sign=+1;
 463       deviation *= sign;
 464       deviation += 15;
 465       Int_t stripX21 = locTrg->LoStripX()+deviation+1;
 466       Float_t y21 = circuit.GetY21Pos(localBoardId, stripX21);
 467       Float_t x11 = circuit.GetX11Pos(localBoardId, locTrg->LoStripY());
 468
 469       AliDebug(1, Form(" MakeTriggerTrack %d %d %d %d %f %f %f \n",locTrg->LoCircuit(),
 470                        locTrg->LoStripX(),locTrg->LoStripX()+locTrg->LoDev()+1,locTrg->LoStripY(),y11, y21, x11));
 471
 472       Float_t thetax = TMath::ATan2( x11 , z11 );
 473       Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) );
 474
 475       triggerTrack.SetX11(x11);
 476       triggerTrack.SetY11(y11);
 477       triggerTrack.SetThetax(thetax);
 478       triggerTrack.SetThetay(thetay);
 479       triggerTrack.SetGTPattern(gloTrigPat);
 480       triggerTrack.SetLoTrgNum(localBoardId);
 481
 482       triggerTrackStore.Add(triggerTrack);
 483     } // board is fired
 484   } // end of loop on Local Trigger
 485 }