MUON/AliMUONEventReconstructor.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 Revision 1.13  2000/10/02 21:28:09  fca
  19 Removal of useless dependecies via forward declarations
  20
  21 Revision 1.12  2000/10/02 16:58:29  egangler
  22 Cleaning of the code :
  23 -> coding conventions
  24 -> void Streamers
  25 -> some useless includes removed or replaced by "class" statement
  26
  27 Revision 1.11  2000/09/22 09:16:33  hristov
  28 Casting needed on DEC
  29
  30 Revision 1.10  2000/09/19 09:49:50  gosset
  31 AliMUONEventReconstructor package
  32 * track extrapolation independent from reco_muon.F, use of AliMagF...
  33 * possibility to use new magnetic field (automatic from generated root file)
  34
  35 Revision 1.9  2000/07/20 12:45:27  gosset
  36 New "EventReconstructor..." structure,
  37         hopefully more adapted to tree/streamer.
  38 "AliMUONEventReconstructor::RemoveDoubleTracks"
  39         to keep only one track among similar ones.
  40
  41 Revision 1.8  2000/07/18 16:04:06  gosset
  42 AliMUONEventReconstructor package:
  43 * a few minor modifications and more comments
  44 * a few corrections
  45   * right sign for Z of raw clusters
  46   * right loop over chambers inside station
  47   * symmetrized covariance matrix for measurements (TrackChi2MCS)
  48   * right sign of charge in extrapolation (ExtrapToZ)
  49   * right zEndAbsorber for Branson correction below 3 degrees
  50 * use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit
  51 * no parameter for AliMUONTrack::Fit() but more fit parameters in Track object
  52
  53 Revision 1.7  2000/07/03 12:28:06  gosset
  54 Printout at the right place after extrapolation to vertex
  55
  56 Revision 1.6  2000/06/30 12:01:06  gosset
  57 Correction for hit search in the right chamber (JPC)
  58
  59 Revision 1.5  2000/06/30 10:15:48  gosset
  60 Changes to EventReconstructor...:
  61 precision fit with multiple Coulomb scattering;
  62 extrapolation to vertex with Branson correction in absorber (JPC)
  63
  64 Revision 1.4  2000/06/27 14:11:36  gosset
  65 Corrections against violations of coding conventions
  66
  67 Revision 1.3  2000/06/16 07:27:08  gosset
  68 To remove problem in running RuleChecker, like in MUON-dev
  69
  70 Revision 1.1.2.5  2000/06/16 07:00:26  gosset
  71 To remove problem in running RuleChecker
  72
  73 Revision 1.1.2.4  2000/06/12 08:00:07  morsch
  74 Dummy streamer to solve CINT compilation problem (to be investigated !)
  75
  76 Revision 1.1.2.3  2000/06/09 20:59:57  morsch
  77 Make includes consistent with new file structure.
  78
  79 Revision 1.1.2.2  2000/06/09 12:58:05  gosset
  80 Removed comment beginnings in Log sections of .cxx files
  81 Suppressed most violations of coding rules
  82
  83 Revision 1.1.2.1  2000/06/07 14:44:53  gosset
  84 Addition of files for track reconstruction in C++
  85 */
  86
  87 //__________________________________________________________________________
  88 //
  89 // MUON event reconstructor in ALICE
  90 //
  91 // This class contains as data:
  92 // * the parameters for the event reconstruction
  93 // * a pointer to the array of hits to be reconstructed (the event)
  94 // * a pointer to the array of segments made with these hits inside each station
  95 // * a pointer to the array of reconstructed tracks
  96 //
  97 // It contains as methods, among others:
  98 // * MakeEventToBeReconstructed to build the array of hits to be reconstructed
  99 // * MakeSegments to build the segments
 100 // * MakeTracks to build the tracks
 101 //__________________________________________________________________________
 102
 103 #include <iostream.h>
 104 #include <stdlib.h>
 105
 106 #include <TRandom.h>
 107 #include <TFile.h>
 108 #include <TTree.h>
 109
 110 #include "AliMUONEventReconstructor.h"
 111 #include "AliMUON.h"
 112 #include "AliMUONHitForRec.h"
 113 #include "AliMUONSegment.h"
 114 #include "AliMUONHit.h"
 115 #include "AliMUONRawCluster.h"
 116 #include "AliMUONTrack.h"
 117 #include "AliMUONChamber.h"
 118 #include "AliMUONTrackHit.h"
 119 #include "AliMagF.h"
 120 #include "AliRun.h"
 121 #include "TParticle.h"
 122
 123 //************* Defaults parameters for reconstruction
 124 static const Double_t kDefaultMinBendingMomentum = 3.0;
 125 static const Double_t kDefaultMaxSigma2Distance = 16.0;
 126 static const Double_t kDefaultBendingResolution = 0.01;
 127 static const Double_t kDefaultNonBendingResolution = 0.144;
 128 static const Double_t kDefaultChamberThicknessInX0 = 0.03;
 129 // Simple magnetic field:
 130 // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
 131 // Length and Position from reco_muon.F, with opposite sign:
 132 // Length = ZMAGEND-ZCOIL
 133 // Position = (ZMAGEND+ZCOIL)/2
 134 // to be ajusted differently from real magnetic field ????
 135 static const Double_t kDefaultSimpleBValue = 7.0;
 136 static const Double_t kDefaultSimpleBLength = 428.0;
 137 static const Double_t kDefaultSimpleBPosition = 1019.0;
 138 static const Int_t kDefaultRecGeantHits = 0;
 139 static const Double_t kDefaultEfficiency = 0.95;
 140
 141 static const Int_t kDefaultPrintLevel = 0;
 142
 143 ClassImp(AliMUONEventReconstructor) // Class implementation in ROOT context
 144
 145   //__________________________________________________________________________
 146 AliMUONEventReconstructor::AliMUONEventReconstructor(void)
 147 {
 148   // Constructor for class AliMUONEventReconstructor
 149   SetReconstructionParametersToDefaults();
 150   // Memory allocation for the TClonesArray of hits for reconstruction
 151   // Is 10000 the right size ????
 152   fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
 153   fNHitsForRec = 0; // really needed or GetEntriesFast sufficient ????
 154   // Memory allocation for the TClonesArray's of segments in stations
 155   // Is 2000 the right size ????
 156   for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
 157     fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
 158     fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
 159   }
 160   // Memory allocation for the TClonesArray of reconstructed tracks
 161   // Is 10 the right size ????
 162   fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
 163   fNRecTracks = 0; // really needed or GetEntriesFast sufficient ????
 164   // Memory allocation for the TClonesArray of hits on reconstructed tracks
 165   // Is 100 the right size ????
 166   fRecTrackHitsPtr = new TClonesArray("AliMUONTrack", 100);
 167   fNRecTrackHits = 0; // really needed or GetEntriesFast sufficient ????
 168
 169   // Sign of fSimpleBValue according to sign of value at (50,50,950).
 170   Float_t b[3], x[3];
 171   x[0] = 50.; x[1] = 50.; x[2] = 950.;
 172   gAlice->Field()->Field(x, b);
 173   fSimpleBValue = TMath::Sign(fSimpleBValue,(Double_t) b[2]);
 174   // See how to get fSimple(BValue, BLength, BPosition)
 175   // automatically calculated from the actual magnetic field ????
 176
 177   if (fPrintLevel >= 0) {
 178     cout << "AliMUONEventReconstructor constructed with defaults" << endl; Dump();
 179     cout << endl << "Magnetic field from root file:" << endl;
 180     gAlice->Field()->Dump();
 181     cout << endl;
 182   }
 183   return;
 184 }
 185
 186 AliMUONEventReconstructor::AliMUONEventReconstructor (const AliMUONEventReconstructor& Reconstructor)
 187 {
 188   // Dummy copy constructor
 189 }
 190
 191 AliMUONEventReconstructor & AliMUONEventReconstructor::operator=(const AliMUONEventReconstructor& Reconstructor)
 192 {
 193   // Dummy assignment operator
 194     return *this;
 195 }
 196
 197   //__________________________________________________________________________
 198 AliMUONEventReconstructor::~AliMUONEventReconstructor(void)
 199 {
 200   // Destructor for class AliMUONEventReconstructor
 201   delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ????
 202   for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
 203     delete fSegmentsPtr[st]; // Correct destruction of everything ????
 204   return;
 205 }
 206
 207   //__________________________________________________________________________
 208 void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
 209 {
 210   // Set reconstruction parameters to default values
 211   // Would be much more convenient with a structure (or class) ????
 212   fMinBendingMomentum = kDefaultMinBendingMomentum;
 213   fMaxSigma2Distance = kDefaultMaxSigma2Distance;
 214
 215   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
 216   // ******** Parameters for making HitsForRec
 217   // minimum radius,
 218   // like in TRACKF_STAT:
 219   // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3;
 220   // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9
 221   for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
 222     if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
 223                   2.0 * TMath::Pi() / 180.0;
 224     else fRMin[ch] = 30.0;
 225   }
 226   // maximum radius
 227   // like in TRACKF_STAT (10 degrees ????)
 228   fRMax[0] = fRMax[1] = 91.5;
 229   fRMax[2] = fRMax[3] = 122.5;
 230   fRMax[4] = fRMax[5] = 158.3;
 231   fRMax[6] = fRMax[7] = 260.0;
 232   fRMax[8] = fRMax[9] = 260.0;
 233
 234   // ******** Parameters for making segments
 235   // should be parametrized ????
 236   // according to interval between chambers in a station ????
 237   // Maximum distance in non bending plane
 238   // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
 239   // SIGCUT*DYMAX(IZ)
 240   for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
 241     fSegmentMaxDistNonBending[st] = 5. * 0.22;
 242   // Maximum distance in bending plane
 243   // values from TRACKF_STAT corresponding to (J psi 20cm)
 244   fSegmentMaxDistBending[0] = 1.5;
 245   fSegmentMaxDistBending[1] = 1.5;
 246   fSegmentMaxDistBending[2] = 3.0;
 247   fSegmentMaxDistBending[3] = 6.0;
 248   fSegmentMaxDistBending[4] = 6.0;
 249
 250   fBendingResolution = kDefaultBendingResolution;
 251   fNonBendingResolution = kDefaultNonBendingResolution;
 252   fChamberThicknessInX0 = kDefaultChamberThicknessInX0;
 253   fSimpleBValue = kDefaultSimpleBValue;
 254   fSimpleBLength = kDefaultSimpleBLength;
 255   fSimpleBPosition = kDefaultSimpleBPosition;
 256   fRecGeantHits = kDefaultRecGeantHits;
 257   fEfficiency = kDefaultEfficiency;
 258   fPrintLevel = kDefaultPrintLevel;
 259   return;
 260 }
 261
 262 //__________________________________________________________________________
 263 Double_t AliMUONEventReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum)
 264 {
 265   // Returns impact parameter at vertex in bending plane (cm),
 266   // from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
 267   // using simple values for dipole magnetic field.
 268   // The sign is the sign of the charge.
 269   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
 270           BendingMomentum);
 271 }
 272
 273 //__________________________________________________________________________
 274 Double_t AliMUONEventReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam)
 275 {
 276   // Returns signed bending momentum in bending plane (GeV/c),
 277   // from the impact parameter "ImpactParam" at vertex in bending plane (cm),
 278   // using simple values for dipole magnetic field.
 279   // The sign is the sign of the charge.
 280   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
 281           ImpactParam);
 282 }
 283
 284 //__________________________________________________________________________
 285 void AliMUONEventReconstructor::SetBkgGeantFile(Text_t *BkgGeantFileName)
 286 {
 287   // Set background file ... for GEANT hits
 288   // Must be called after having loaded the firts signal event
 289   if (fPrintLevel >= 0) {
 290     cout << "Enter SetBkgGeantFile with BkgGeantFileName ``"
 291          << BkgGeantFileName << "''" << endl;}
 292   if (strlen(BkgGeantFileName)) {
 293     // BkgGeantFileName not empty: try to open the file
 294     if (fPrintLevel >= 2) {cout << "Before File(Bkg)" << endl; gDirectory->Dump();}
 295     fBkgGeantFile = new TFile(BkgGeantFileName);
 296     if (fPrintLevel >= 2) {cout << "After File(Bkg)" << endl; gDirectory->Dump();}
 297     if (fBkgGeantFile-> IsOpen()) {
 298       if (fPrintLevel >= 0) {
 299         cout << "Background for GEANT hits in file: ``" << BkgGeantFileName
 300              << "'' successfully opened" << endl;}
 301     }
 302     else {
 303       cout << "Background for GEANT hits in file: " << BkgGeantFileName << endl;
 304       cout << "NOT FOUND: EXIT" << endl;
 305       exit(0); // right instruction for exit ????
 306     }
 307     // Arrays for "particles" and "hits"
 308     fBkgGeantParticles = new TClonesArray("TParticle", 200);
 309     fBkgGeantHits = new TClonesArray("AliMUONHit", 2000);
 310     // Event number to -1 for initialization
 311     fBkgGeantEventNumber = -1;
 312     // Back to the signal file:
 313     // first signal event must have been loaded previously,
 314     // otherwise, Segmentation violation at the next instruction
 315     // How is it possible to do smething better ????
 316     ((gAlice->TreeK())->GetCurrentFile())->cd();
 317     if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();}
 318   }
 319   return;
 320 }
 321
 322 //__________________________________________________________________________
 323 void AliMUONEventReconstructor::NextBkgGeantEvent(void)
 324 {
 325   // Get next event in background file for GEANT hits
 326   // Goes back to event number 0 when end of file is reached
 327   char treeName[20];
 328   TBranch *branch;
 329   if (fPrintLevel >= 0) {
 330     cout << "Enter NextBkgGeantEvent" << endl;}
 331   // Clean previous event
 332   if(fBkgGeantTK) delete fBkgGeantTK;
 333   fBkgGeantTK = NULL;
 334   if(fBkgGeantParticles) fBkgGeantParticles->Clear();
 335   if(fBkgGeantTH) delete fBkgGeantTH;
 336   fBkgGeantTH = NULL;
 337   if(fBkgGeantHits) fBkgGeantHits->Clear();
 338   // Increment event number
 339   fBkgGeantEventNumber++;
 340   // Get access to Particles and Hits for event from background file
 341   if (fPrintLevel >= 2) {cout << "Before cd(Bkg)" << endl; gDirectory->Dump();}
 342   fBkgGeantFile->cd();
 343   if (fPrintLevel >= 2) {cout << "After cd(Bkg)" << endl; gDirectory->Dump();}
 344   // Particles: TreeK for event and branch "Particles"
 345   sprintf(treeName, "TreeK%d", fBkgGeantEventNumber);
 346   fBkgGeantTK = (TTree*)gDirectory->Get(treeName);
 347   if (!fBkgGeantTK) {
 348     if (fPrintLevel >= 0) {
 349       cout << "Cannot find Kine Tree for background event: " <<
 350         fBkgGeantEventNumber << endl;
 351       cout << "Goes back to event 0" << endl;
 352     }
 353     fBkgGeantEventNumber = 0;
 354     sprintf(treeName, "TreeK%d", fBkgGeantEventNumber);
 355     fBkgGeantTK = (TTree*)gDirectory->Get(treeName);
 356     if (!fBkgGeantTK) {
 357       cout << "ERROR: cannot find Kine Tree for background event: " <<
 358         fBkgGeantEventNumber << endl;
 359       exit(0);
 360     }
 361   }
 362   if (fBkgGeantTK)
 363     fBkgGeantTK->SetBranchAddress("Particles", &fBkgGeantParticles);
 364   fBkgGeantTK->GetEvent(0); // why event 0 ???? necessary ????
 365   // Hits: TreeH for event and branch "MUON"
 366   sprintf(treeName, "TreeH%d", fBkgGeantEventNumber);
 367   fBkgGeantTH = (TTree*)gDirectory->Get(treeName);
 368   if (!fBkgGeantTH) {
 369     cout << "ERROR: cannot find Hits Tree for background event: " <<
 370       fBkgGeantEventNumber << endl;
 371       exit(0);
 372   }
 373   if (fBkgGeantTH && fBkgGeantHits) {
 374     branch = fBkgGeantTH->GetBranch("MUON");
 375     if (branch) branch->SetAddress(&fBkgGeantHits);
 376   }
 377   fBkgGeantTH->GetEntries(); // necessary ????
 378   // Back to the signal file
 379   ((gAlice->TreeK())->GetCurrentFile())->cd();
 380   if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();}
 381   return;
 382 }
 383
 384 //__________________________________________________________________________
 385 void AliMUONEventReconstructor::EventReconstruct(void)
 386 {
 387   // To reconstruct one event
 388   if (fPrintLevel >= 1) cout << "enter EventReconstruct" << endl;
 389   MakeEventToBeReconstructed();
 390   MakeSegments();
 391   MakeTracks();
 392   return;
 393 }
 394
 395   //__________________________________________________________________________
 396 void AliMUONEventReconstructor::ResetHitsForRec(void)
 397 {
 398   // To reset the array and the number of HitsForRec,
 399   // and also the number of HitsForRec
 400   // and the index of the first HitForRec per chamber
 401   if (fHitsForRecPtr) fHitsForRecPtr->Clear();
 402   fNHitsForRec = 0;
 403   for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++)
 404     fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
 405   return;
 406 }
 407
 408   //__________________________________________________________________________
 409 void AliMUONEventReconstructor::ResetSegments(void)
 410 {
 411   // To reset the TClonesArray of segments and the number of Segments
 412   // for all stations
 413   for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
 414     if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
 415     fNSegments[st] = 0;
 416   }
 417   return;
 418 }
 419
 420   //__________________________________________________________________________
 421 void AliMUONEventReconstructor::ResetTracks(void)
 422 {
 423   // To reset the TClonesArray of reconstructed tracks
 424   if (fRecTracksPtr) fRecTracksPtr->Delete();
 425   // Delete in order that the Track destructors are called,
 426   // hence the space for the TClonesArray of pointers to TrackHit's is freed
 427   fNRecTracks = 0;
 428   return;
 429 }
 430
 431   //__________________________________________________________________________
 432 void AliMUONEventReconstructor::ResetTrackHits(void)
 433 {
 434   // To reset the TClonesArray of hits on reconstructed tracks
 435   if (fRecTrackHitsPtr) fRecTrackHitsPtr->Clear();
 436   fNRecTrackHits = 0;
 437   return;
 438 }
 439
 440   //__________________________________________________________________________
 441 void AliMUONEventReconstructor::MakeEventToBeReconstructed(void)
 442 {
 443   // To make the list of hits to be reconstructed,
 444   // either from the GEANT hits or from the raw clusters
 445   // according to the parameter set for the reconstructor
 446   if (fPrintLevel >= 1) cout << "enter MakeEventToBeReconstructed" << endl;
 447   ResetHitsForRec();
 448   if (fRecGeantHits == 1) {
 449     // Reconstruction from GEANT hits
 450     // Back to the signal file
 451     ((gAlice->TreeK())->GetCurrentFile())->cd();
 452     // Signal hits
 453     // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 454     // Security on MUON ????
 455     AddHitsForRecFromGEANT(gAlice->TreeH());
 456     // Background hits
 457     AddHitsForRecFromBkgGEANT(fBkgGeantTH, fBkgGeantHits);
 458     // Sort HitsForRec in increasing order with respect to chamber number
 459     SortHitsForRecWithIncreasingChamber();
 460   }
 461   else {
 462     // Reconstruction from raw clusters
 463     // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 464     // Security on MUON ????
 465     // TreeR assumed to be be "prepared" in calling function
 466     // by "MUON->GetTreeR(nev)" ????
 467     TTree *treeR = gAlice->TreeR();
 468     AddHitsForRecFromRawClusters(treeR);
 469     // No sorting: it is done automatically in the previous function
 470   }
 471   if (fPrintLevel >= 10) {
 472     cout << "end of MakeEventToBeReconstructed" << endl;
 473     cout << "NHitsForRec: " << fNHitsForRec << endl;
 474     for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
 475       cout << "chamber(0...): " << ch
 476            << "  NHitsForRec: " << fNHitsForRecPerChamber[ch]
 477            << "  index(first HitForRec): " << fIndexOfFirstHitForRecPerChamber[ch]
 478            << endl;
 479       for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
 480            hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
 481            hit++) {
 482         cout << "HitForRec index(0...): " << hit << endl;
 483         ((*fHitsForRecPtr)[hit])->Dump();
 484       }
 485     }
 486   }
 487   return;
 488 }
 489
 490   //__________________________________________________________________________
 491 void AliMUONEventReconstructor::AddHitsForRecFromGEANT(TTree *TH)
 492 {
 493   // To add to the list of hits for reconstruction
 494   // the GEANT signal hits from a hit tree TH.
 495   if (fPrintLevel >= 2)
 496     cout << "enter AddHitsForRecFromGEANT with TH: " << TH << endl;
 497   if (TH == NULL) return;
 498   AliMUON *pMUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 499   // Security on MUON ????
 500   // See whether it could be the same for signal and background ????
 501   // Loop over tracks in tree
 502   Int_t ntracks = (Int_t) TH->GetEntries();
 503   if (fPrintLevel >= 2)
 504     cout << "ntracks: " << ntracks << endl;
 505   for (Int_t track = 0; track < ntracks; track++) {
 506     gAlice->ResetHits();
 507     TH->GetEvent(track);
 508     // Loop over hits
 509     Int_t hit = 0;
 510     for (AliMUONHit* mHit = (AliMUONHit*) pMUON->FirstHit(-1);
 511          mHit;
 512          mHit = (AliMUONHit*) pMUON->NextHit(), hit++) {
 513       NewHitForRecFromGEANT(mHit,track, hit, 1);
 514     } // end of hit loop
 515   } // end of track loop
 516   return;
 517 }
 518
 519   //__________________________________________________________________________
 520 void AliMUONEventReconstructor::AddHitsForRecFromBkgGEANT(TTree *TH, TClonesArray *Hits)
 521 {
 522   // To add to the list of hits for reconstruction
 523   // the GEANT background hits from a hit tree TH and a pointer Hits to a hit list.
 524   // How to have only one function "AddHitsForRecFromGEANT" ????
 525   if (fPrintLevel >= 2)
 526     cout << "enter AddHitsForRecFromBkgGEANT with TH: " << TH << endl;
 527   if (TH == NULL) return;
 528   // Loop over tracks in tree
 529   Int_t ntracks = (Int_t) TH->GetEntries();
 530   if (fPrintLevel >= 2)
 531     cout << "ntracks: " << ntracks << endl;
 532   for (Int_t track = 0; track < ntracks; track++) {
 533     if (Hits) Hits->Clear();
 534     TH->GetEvent(track);
 535     // Loop over hits
 536     for (Int_t hit = 0; hit < Hits->GetEntriesFast(); hit++) {
 537       NewHitForRecFromGEANT((AliMUONHit*) (*Hits)[hit], track, hit, 0);
 538     } // end of hit loop
 539   } // end of track loop
 540   return;
 541 }
 542
 543   //__________________________________________________________________________
 544 AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* Hit, Int_t TrackNumber, Int_t HitNumber, Int_t Signal)
 545 {
 546   // To make a new hit for reconstruction from a GEANT hit pointed to by "Hit",
 547   // with hit number "HitNumber" in the track numbered "TrackNumber",
 548   // either from signal ("Signal" = 1) or background ("Signal" = 0) event.
 549   // Selects hits in tracking (not trigger) chambers.
 550   // Takes into account the efficiency (fEfficiency)
 551   // and the smearing from resolution (fBendingResolution and fNonBendingResolution).
 552   // Adds a condition on the radius between RMin and RMax
 553   // to better simulate the real chambers.
 554   // Returns the pointer to the new hit for reconstruction,
 555   // or NULL in case of inefficiency or non tracking chamber or bad radius.
 556   // No condition on at most 20 cm from a muon from a resonance
 557   // like in Fortran TRACKF_STAT.
 558   AliMUONHitForRec* hitForRec;
 559   Double_t bendCoor, nonBendCoor, radius;
 560   Int_t chamber = Hit->fChamber - 1; // chamber(0...)
 561   // only in tracking chambers (fChamber starts at 1)
 562   if (chamber >= kMaxMuonTrackingChambers) return NULL;
 563   // only if hit is efficient (keep track for checking ????)
 564   if (gRandom->Rndm() > fEfficiency) return NULL;
 565   // only if radius between RMin and RMax
 566   bendCoor = Hit->Y();
 567   nonBendCoor = Hit->X();
 568   radius = TMath::Sqrt((bendCoor * bendCoor) + (nonBendCoor * nonBendCoor));
 569   if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL;
 570   // new AliMUONHitForRec from GEANT hit and increment number of AliMUONHitForRec's
 571   hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(Hit);
 572   fNHitsForRec++;
 573   // add smearing from resolution
 574   hitForRec->SetBendingCoor(bendCoor + gRandom->Gaus(0., fBendingResolution));
 575   hitForRec->SetNonBendingCoor(nonBendCoor
 576                                + gRandom->Gaus(0., fNonBendingResolution));
 577   // more information into HitForRec
 578   //  resolution: angular effect to be added here ????
 579   hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
 580   hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
 581   //  GEANT track info
 582   hitForRec->SetHitNumber(HitNumber);
 583   hitForRec->SetTHTrack(TrackNumber);
 584   hitForRec->SetGeantSignal(Signal);
 585   if (fPrintLevel >= 10) {
 586     cout << "track: " << TrackNumber << " hit: " << HitNumber << endl;
 587     Hit->Dump();
 588     cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 589     hitForRec->Dump();}
 590   return hitForRec;
 591 }
 592
 593   //__________________________________________________________________________
 594 void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber()
 595 {
 596   // Sort HitsForRec's in increasing order with respect to chamber number.
 597   // Uses the function "Compare".
 598   // Update the information for HitsForRec per chamber too.
 599   Int_t ch, nhits, prevch;
 600   fHitsForRecPtr->Sort();
 601   for (ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
 602     fNHitsForRecPerChamber[ch] = 0;
 603     fIndexOfFirstHitForRecPerChamber[ch] = 0;
 604   }
 605   prevch = 0; // previous chamber
 606   nhits = 0; // number of hits in current chamber
 607   // Loop over HitsForRec
 608   for (Int_t hit = 0; hit < fNHitsForRec; hit++) {
 609     // chamber number (0...)
 610     ch = ((AliMUONHitForRec*)  ((*fHitsForRecPtr)[hit]))->GetChamberNumber();
 611     // increment number of hits in current chamber
 612     (fNHitsForRecPerChamber[ch])++;
 613     // update index of first HitForRec in current chamber
 614     // if chamber number different from previous one
 615     if (ch != prevch) {
 616       fIndexOfFirstHitForRecPerChamber[ch] = hit;
 617       prevch = ch;
 618     }
 619   }
 620   return;
 621 }
 622
 623 //   //__________________________________________________________________________
 624 // void AliMUONEventReconstructor::AddHitsForRecFromCathodeCorrelations(TTree* TC)
 625 // {
 626 //   // OLD VERSION WHEN ONE ONE WAS USING SO CALLED CATHODE CORRELATIONS
 627 //   // To add to the list of hits for reconstruction
 628 //   // the (cathode correlated) raw clusters
 629 //   // No condition added, like in Fortran TRACKF_STAT,
 630 //   // on the radius between RMin and RMax.
 631 //   AliMUONHitForRec *hitForRec;
 632 //   if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromCathodeCorrelations" << endl;
 633 //   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 634 //   // Security on MUON ????
 635 //   // Loop over tracking chambers
 636 //   for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
 637 //     // number of HitsForRec to 0 for the chamber
 638 //     fNHitsForRecPerChamber[ch] = 0;
 639 //     // index of first HitForRec for the chamber
 640 //     if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
 641 //     else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
 642 //     TClonesArray *reconst_hits  = MUON->ReconstHitsAddress(ch);
 643 //     MUON->ResetReconstHits();
 644 //     TC->GetEvent();
 645 //     Int_t ncor = (Int_t)reconst_hits->GetEntries();
 646 //     // Loop over (cathode correlated) raw clusters
 647 //     for (Int_t cor = 0; cor < ncor; cor++) {
 648 //       AliMUONReconstHit * mCor =
 649 //      (AliMUONReconstHit*) reconst_hits->UncheckedAt(cor);
 650 //       // new AliMUONHitForRec from (cathode correlated) raw cluster
 651 //       // and increment number of AliMUONHitForRec's (total and in chamber)
 652 //       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(mCor);
 653 //       fNHitsForRec++;
 654 //       (fNHitsForRecPerChamber[ch])++;
 655 //       // more information into HitForRec
 656 //       hitForRec->SetChamberNumber(ch);
 657 //       hitForRec->SetHitNumber(cor);
 658 //       // Z coordinate of the chamber (cm) with sign opposite to GEANT sign
 659 //       // could (should) be more exact from chamber geometry ????
 660 //       hitForRec->SetZ(-(&(MUON->Chamber(ch)))->Z());
 661 //       if (fPrintLevel >= 10) {
 662 //      cout << "chamber (0...): " << ch <<
 663 //        " cathcorrel (0...): " << cor << endl;
 664 //      mCor->Dump();
 665 //      cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 666 //      hitForRec->Dump();}
 667 //     } // end of cluster loop
 668 //   } // end of chamber loop
 669 //   return;
 670 // }
 671
 672   //__________________________________________________________________________
 673 void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
 674 {
 675   // To add to the list of hits for reconstruction all the raw clusters
 676   // No condition added, like in Fortran TRACKF_STAT,
 677   // on the radius between RMin and RMax.
 678   AliMUONHitForRec *hitForRec;
 679   AliMUONRawCluster *clus;
 680   Int_t iclus, nclus;
 681   TClonesArray *rawclusters;
 682   if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromRawClusters" << endl;
 683   AliMUON *pMUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 684   // Security on MUON ????
 685   // Loop over tracking chambers
 686   for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
 687     // number of HitsForRec to 0 for the chamber
 688     fNHitsForRecPerChamber[ch] = 0;
 689     // index of first HitForRec for the chamber
 690     if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
 691     else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
 692     rawclusters = pMUON->RawClustAddress(ch);
 693     pMUON->ResetRawClusters();
 694     TR->GetEvent((Int_t) (TR->GetEntries()) - 1); // to be checked ????
 695     nclus = (Int_t) (rawclusters->GetEntries());
 696     // Loop over (cathode correlated) raw clusters
 697     for (iclus = 0; iclus < nclus; iclus++) {
 698       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
 699       // new AliMUONHitForRec from raw cluster
 700       // and increment number of AliMUONHitForRec's (total and in chamber)
 701       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
 702       fNHitsForRec++;
 703       (fNHitsForRecPerChamber[ch])++;
 704       // more information into HitForRec
 705       //  resolution: info should be already in raw cluster and taken from it ????
 706       hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
 707       hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
 708       //  original raw cluster
 709       hitForRec->SetChamberNumber(ch);
 710       hitForRec->SetHitNumber(iclus);
 711       // Z coordinate of the chamber (cm)
 712       // could (should) be more exact from chamber geometry ????
 713       hitForRec->SetZ((&(pMUON->Chamber(ch)))->Z());
 714       if (fPrintLevel >= 10) {
 715         cout << "chamber (0...): " << ch <<
 716           " raw cluster (0...): " << iclus << endl;
 717         clus->Dump();
 718         cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 719         hitForRec->Dump();}
 720     } // end of cluster loop
 721   } // end of chamber loop
 722   return;
 723 }
 724
 725   //__________________________________________________________________________
 726 void AliMUONEventReconstructor::MakeSegments(void)
 727 {
 728   // To make the list of segments in all stations,
 729   // from the list of hits to be reconstructed
 730   if (fPrintLevel >= 1) cout << "enter MakeSegments" << endl;
 731   ResetSegments();
 732   // Loop over stations
 733   for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
 734     MakeSegmentsPerStation(st);
 735   if (fPrintLevel >= 10) {
 736     cout << "end of MakeSegments" << endl;
 737     for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
 738       cout << "station(0...): " << st
 739            << "  Segments: " << fNSegments[st]
 740            << endl;
 741       for (Int_t seg = 0;
 742            seg < fNSegments[st];
 743            seg++) {
 744         cout << "Segment index(0...): " << seg << endl;
 745         ((*fSegmentsPtr[st])[seg])->Dump();
 746       }
 747     }
 748   }
 749   return;
 750 }
 751
 752   //__________________________________________________________________________
 753 void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station)
 754 {
 755   // To make the list of segments in station number "Station" (0...)
 756   // from the list of hits to be reconstructed.
 757   // Updates "fNSegments"[Station].
 758   // Segments in stations 4 and 5 are sorted
 759   // according to increasing absolute value of "impact parameter"
 760   AliMUONHitForRec *hit1Ptr, *hit2Ptr;
 761   AliMUONSegment *segment;
 762   Bool_t last2st;
 763   Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
 764       impactParam = 0., maxImpactParam = 0.; // =0 to avoid compilation warnings.
 765   AliMUON *pMUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 766   if (fPrintLevel >= 1)
 767     cout << "enter MakeSegmentsPerStation (0...) " << Station << endl;
 768   // first and second chambers (0...) in the station
 769   Int_t ch1 = 2 * Station;
 770   Int_t ch2 = ch1 + 1;
 771   // variable true for stations downstream of the dipole:
 772   // Station(0..4) equal to 3 or 4
 773   if ((Station == 3) || (Station == 4)) {
 774     last2st = kTRUE;
 775     // maximum impact parameter (cm) according to fMinBendingMomentum...
 776     maxImpactParam =
 777       TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
 778   }
 779   else last2st = kFALSE;
 780   // extrapolation factor from Z of first chamber to Z of second chamber
 781   // dZ to be changed to take into account fine structure of chambers ????
 782   Double_t extrapFact =
 783     (&(pMUON->Chamber(ch2)))->Z() / (&(pMUON->Chamber(ch1)))->Z();
 784   // index for current segment
 785   Int_t segmentIndex = 0;
 786   // Loop over HitsForRec in the first chamber of the station
 787   for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
 788        hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
 789        hit1++) {
 790     // pointer to the HitForRec
 791     hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
 792     // extrapolation,
 793     // on the straight line joining the HitForRec to the vertex (0,0,0),
 794     // to the Z of the second chamber of the station
 795     extBendCoor = extrapFact * hit1Ptr->GetBendingCoor();
 796     extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor();
 797     // Loop over HitsForRec in the second chamber of the station
 798     for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
 799          hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
 800          hit2++) {
 801       // pointer to the HitForRec
 802       hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]);
 803       // absolute values of distances, in bending and non bending planes,
 804       // between the HitForRec in the second chamber
 805       // and the previous extrapolation
 806       distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor);
 807       distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor);
 808       if (last2st) {
 809         // bending slope
 810         bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) /
 811           (hit1Ptr->GetZ() - hit2Ptr->GetZ());
 812         // absolute value of impact parameter
 813         impactParam =
 814           TMath::Abs(hit1Ptr->GetBendingCoor() - hit2Ptr->GetZ() * bendingSlope);
 815       }
 816       // check for distances not too large,
 817       // and impact parameter not too big if stations downstream of the dipole.
 818       // Conditions "distBend" and "impactParam" correlated for these stations ????
 819       if ((distBend < fSegmentMaxDistBending[Station]) &&
 820           (distNonBend < fSegmentMaxDistNonBending[Station]) &&
 821           (!last2st || (impactParam < maxImpactParam))) {
 822         // make new segment
 823         segment = new ((*fSegmentsPtr[Station])[segmentIndex])
 824           AliMUONSegment(hit1Ptr, hit2Ptr);
 825         // update "link" to this segment from the hit in the first chamber
 826         if (hit1Ptr->GetNSegments() == 0)
 827           hit1Ptr->SetIndexOfFirstSegment(segmentIndex);
 828         hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1);
 829         if (fPrintLevel >= 10) {
 830           cout << "segmentIndex(0...): " << segmentIndex
 831                << "  distBend: " << distBend
 832                << "  distNonBend: " << distNonBend
 833                << endl;
 834           segment->Dump();
 835           cout << "HitForRec in first chamber" << endl;
 836           hit1Ptr->Dump();
 837           cout << "HitForRec in second chamber" << endl;
 838           hit2Ptr->Dump();
 839         };
 840         // increment index for current segment
 841         segmentIndex++;
 842       }
 843     } //for (Int_t hit2
 844   } // for (Int_t hit1...
 845   fNSegments[Station] = segmentIndex;
 846   // Sorting according to "impact parameter" if station(1..5) 4 or 5,
 847   // i.e. Station(0..4) 3 or 4, using the function "Compare".
 848   // After this sorting, it is impossible to use
 849   // the "fNSegments" and "fIndexOfFirstSegment"
 850   // of the HitForRec in the first chamber to explore all segments formed with it.
 851   // Is this sorting really needed ????
 852   if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort();
 853   if (fPrintLevel >= 1) cout << "Station: " << Station << "  NSegments: "
 854                              << fNSegments[Station] << endl;
 855   return;
 856 }
 857
 858   //__________________________________________________________________________
 859 void AliMUONEventReconstructor::MakeTracks(void)
 860 {
 861   // To make the tracks,
 862   // from the list of segments and points in all stations
 863   if (fPrintLevel >= 1) cout << "enter MakeTracks" << endl;
 864   // The order may be important for the following Reset's
 865   ResetTracks();
 866   ResetTrackHits();
 867   // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
 868   MakeTrackCandidates();
 869   // Follow tracks in stations(1..) 3, 2 and 1
 870   FollowTracks();
 871   // Remove double tracks
 872   RemoveDoubleTracks();
 873   return;
 874 }
 875
 876   //__________________________________________________________________________
 877 Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
 878 {
 879   // To make track candidates with two segments in stations(1..) 4 and 5,
 880   // the first segment being pointed to by "BegSegment".
 881   // Returns the number of such track candidates.
 882   Int_t endStation, iEndSegment, nbCan2Seg;
 883   AliMUONSegment *endSegment, *extrapSegment;
 884   AliMUONTrack *recTrack;
 885   Double_t mcsFactor;
 886   if (fPrintLevel >= 1) cout << "enter MakeTrackCandidatesWithTwoSegments" << endl;
 887   // Station for the end segment
 888   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
 889   // multiple scattering factor corresponding to one chamber
 890   mcsFactor = 0.0136 /
 891     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
 892   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 893   // linear extrapolation to end station
 894   extrapSegment =
 895     BegSegment->CreateSegmentFromLinearExtrapToStation(endStation, mcsFactor);
 896   // number of candidates with 2 segments to 0
 897   nbCan2Seg = 0;
 898   // Loop over segments in the end station
 899   for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
 900     // pointer to segment
 901     endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
 902     // test compatibility between current segment and "extrapSegment"
 903     // 4 because 4 quantities in chi2
 904     if ((endSegment->
 905          NormalizedChi2WithSegment(extrapSegment,
 906                                    fMaxSigma2Distance)) <= 4.0) {
 907       // both segments compatible:
 908       // make track candidate from "begSegment" and "endSegment"
 909       if (fPrintLevel >= 2)
 910         cout << "TrackCandidate with Segment " << iEndSegment <<
 911           " in Station(0..) " << endStation << endl;
 912       // flag for both segments in one track:
 913       // to be done in track constructor ????
 914       BegSegment->SetInTrack(kTRUE);
 915       endSegment->SetInTrack(kTRUE);
 916       recTrack = new ((*fRecTracksPtr)[fNRecTracks])
 917         AliMUONTrack(BegSegment, endSegment, this);
 918       fNRecTracks++;
 919       if (fPrintLevel >= 10) recTrack->RecursiveDump();
 920       // increment number of track candidates with 2 segments
 921       nbCan2Seg++;
 922     }
 923   } // for (iEndSegment = 0;...
 924   delete extrapSegment; // should not delete HitForRec's it points to !!!!
 925   return nbCan2Seg;
 926 }
 927
 928   //__________________________________________________________________________
 929 Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
 930 {
 931   // To make track candidates with one segment and one point
 932   // in stations(1..) 4 and 5,
 933   // the segment being pointed to by "BegSegment".
 934   Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
 935   AliMUONHitForRec *extrapHitForRec, *hit;
 936   AliMUONTrack *recTrack;
 937   Double_t mcsFactor;
 938   if (fPrintLevel >= 1)
 939     cout << "enter MakeTrackCandidatesWithOneSegmentAndOnePoint" << endl;
 940   // station for the end point
 941   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
 942   // multiple scattering factor corresponding to one chamber
 943   mcsFactor = 0.0136 /
 944     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
 945   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 946   // first and second chambers(0..) in the end station
 947   ch1 = 2 * endStation;
 948   ch2 = ch1 + 1;
 949   // number of candidates to 0
 950   nbCan1Seg1Hit = 0;
 951   // Loop over chambers of the end station
 952   for (ch = ch2; ch >= ch1; ch--) {
 953     // linear extrapolation to chamber
 954     extrapHitForRec =
 955       BegSegment->CreateHitForRecFromLinearExtrapToChamber(ch, mcsFactor);
 956     // limits for the hit index in the loop
 957     iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
 958     iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
 959     // Loop over HitForRec's in the chamber
 960     for (iHit = iHitMin; iHit < iHitMax; iHit++) {
 961       // pointer to HitForRec
 962       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
 963       // test compatibility between current HitForRec and "extrapHitForRec"
 964       // 2 because 2 quantities in chi2
 965       if ((hit->
 966            NormalizedChi2WithHitForRec(extrapHitForRec,
 967                                        fMaxSigma2Distance)) <= 2.0) {
 968         // both HitForRec's compatible:
 969         // make track candidate from begSegment and current HitForRec
 970         if (fPrintLevel >= 2)
 971           cout << "TrackCandidate with HitForRec " << iHit <<
 972             " in Chamber(0..) " << ch << endl;
 973         // flag for beginning segments in one track:
 974         // to be done in track constructor ????
 975         BegSegment->SetInTrack(kTRUE);
 976         recTrack = new ((*fRecTracksPtr)[fNRecTracks])
 977           AliMUONTrack(BegSegment, hit, this);
 978         // the right place to eliminate "double counting" ???? how ????
 979         fNRecTracks++;
 980         if (fPrintLevel >= 10) recTrack->RecursiveDump();
 981         // increment number of track candidates
 982         nbCan1Seg1Hit++;
 983       }
 984     } // for (iHit = iHitMin;...
 985     delete extrapHitForRec;
 986   } // for (ch = ch2;...
 987   return nbCan1Seg1Hit;
 988 }
 989
 990   //__________________________________________________________________________
 991 void AliMUONEventReconstructor::MakeTrackCandidates(void)
 992 {
 993   // To make track candidates
 994   // with at least 3 aligned points in stations(1..) 4 and 5
 995   // (two Segment's or one Segment and one HitForRec)
 996   Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
 997   AliMUONSegment *begSegment;
 998   if (fPrintLevel >= 1) cout << "enter MakeTrackCandidates" << endl;
 999   // Loop over stations(1..) 5 and 4 for the beginning segment
1000   for (begStation = 4; begStation > 2; begStation--) {
1001     // Loop over segments in the beginning station
1002     for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
1003       // pointer to segment
1004       begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
1005       if (fPrintLevel >= 2)
1006         cout << "look for TrackCandidate's with Segment " << iBegSegment <<
1007           " in Station(0..) " << begStation << endl;
1008       // Look for track candidates with two segments,
1009       // "begSegment" and all compatible segments in other station.
1010       // Only for beginning station(1..) 5
1011       // because candidates with 2 segments have to looked for only once.
1012       if (begStation == 4)
1013         nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
1014       // Look for track candidates with one segment and one point,
1015       // "begSegment" and all compatible HitForRec's in other station.
1016       // Only if "begSegment" does not belong already to a track candidate.
1017       // Is that a too strong condition ????
1018       if (!(begSegment->GetInTrack()))
1019         nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
1020     } // for (iBegSegment = 0;...
1021   } // for (begStation = 4;...
1022   return;
1023 }
1024
1025   //__________________________________________________________________________
1026 void AliMUONEventReconstructor::FollowTracks(void)
1027 {
1028   // Follow tracks in stations(1..) 3, 2 and 1
1029   // too long: should be made more modular !!!!
1030   AliMUONHitForRec *bestHit, *extrapHit, *hit;
1031   AliMUONSegment *bestSegment, *extrapSegment, *segment;
1032   AliMUONTrack *track, *nextTrack;
1033   AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
1034   // -1 to avoid compilation warnings
1035   Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex;
1036   Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
1037   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
1038   // local maxSigma2Distance, for easy increase in testing
1039   maxSigma2Distance = fMaxSigma2Distance;
1040   if (fPrintLevel >= 2)
1041     cout << "enter FollowTracks" << endl;
1042   // Loop over track candidates
1043   track = (AliMUONTrack*) fRecTracksPtr->First();
1044   trackIndex = -1;
1045   while (track) {
1046     // Follow function for each track candidate ????
1047     trackIndex++;
1048     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
1049     if (fPrintLevel >= 2)
1050       cout << "FollowTracks: track candidate(0..): " << trackIndex << endl;
1051     // Fit track candidate
1052     track->SetFitMCS(0); // without multiple Coulomb scattering
1053     track->SetFitNParam(3); // with 3 parameters (X = Y = 0)
1054     track->SetFitStart(0); // from parameters at vertex
1055     track->Fit();
1056     if (fPrintLevel >= 10) {
1057       cout << "FollowTracks: track candidate(0..): " << trackIndex
1058            << " after fit in stations(0..) 3 and 4" << endl;
1059       track->RecursiveDump();
1060     }
1061     // Loop over stations(1..) 3, 2 and 1
1062     // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
1063     // otherwise: majority coincidence 2 !!!!
1064     for (station = 2; station >= 0; station--) {
1065       // Track parameters at first track hit (smallest Z)
1066       trackParam1 = ((AliMUONTrackHit*)
1067                      (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1068       // extrapolation to station
1069       trackParam1->ExtrapToStation(station, trackParam);
1070       extrapSegment = new AliMUONSegment(); //  empty segment
1071       // multiple scattering factor corresponding to one chamber
1072       // and momentum in bending plane (not total)
1073       mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
1074       mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
1075       // Z difference from previous station
1076       dZ1 = (&(pMUON->Chamber(2 * station)))->Z() -
1077         (&(pMUON->Chamber(2 * station + 2)))->Z();
1078       // Z difference between the two previous stations
1079       dZ2 = (&(pMUON->Chamber(2 * station + 2)))->Z() -
1080         (&(pMUON->Chamber(2 * station + 4)))->Z();
1081       // Z difference between the two chambers in the previous station
1082       dZ3 = (&(pMUON->Chamber(2 * station)))->Z() -
1083         (&(pMUON->Chamber(2 * station + 1)))->Z();
1084       extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
1085       extrapSegment->
1086         SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
1087       extrapSegment->UpdateFromStationTrackParam
1088         (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
1089          trackParam1->GetInverseBendingMomentum());
1090       bestChi2 = 5.0;
1091       bestSegment = NULL;
1092       if (fPrintLevel >= 10) {
1093         cout << "FollowTracks: track candidate(0..): " << trackIndex
1094              << " Look for segment in station(0..): " << station << endl;
1095       }
1096       // Loop over segments in station
1097       for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
1098         // Look for best compatible Segment in station
1099         // should consider all possibilities ????
1100         // multiple scattering ????
1101         // separation in 2 functions: Segment and HitForRec ????
1102         segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
1103         chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
1104         if (chi2 < bestChi2) {
1105           // update best Chi2 and Segment if better found
1106           bestSegment = segment;
1107           bestChi2 = chi2;
1108         }
1109       }
1110       if (bestSegment) {
1111         // best segment found: add it to track candidate
1112         track->AddSegment(bestSegment);
1113         // set track parameters at these two TrakHit's
1114         track->SetTrackParamAtHit(track->GetNTrackHits() - 2, &(trackParam[0]));
1115         track->SetTrackParamAtHit(track->GetNTrackHits() - 1, &(trackParam[1]));
1116         if (fPrintLevel >= 10) {
1117           cout << "FollowTracks: track candidate(0..): " << trackIndex
1118                << " Added segment in station(0..): " << station << endl;
1119           track->RecursiveDump();
1120         }
1121       }
1122       else {
1123         // No best segment found:
1124         // Look for best compatible HitForRec in station:
1125         // should consider all possibilities ????
1126         // multiple scattering ???? do about like for extrapSegment !!!!
1127         extrapHit = new AliMUONHitForRec(); //  empty hit
1128         bestChi2 = 3.0;
1129         bestHit = NULL;
1130         if (fPrintLevel >= 10) {
1131           cout << "FollowTracks: track candidate(0..): " << trackIndex
1132                << " Segment not found, look for hit in station(0..): " << station
1133                << endl;
1134         }
1135         // Loop over chambers of the station
1136         for (chInStation = 0; chInStation < 2; chInStation++) {
1137           // coordinates of extrapolated hit
1138           extrapHit->
1139             SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
1140           extrapHit->
1141             SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
1142           // resolutions from "extrapSegment"
1143           extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
1144           extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
1145           // Loop over hits in the chamber
1146           ch = 2 * station + chInStation;
1147           for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
1148                iHit < fIndexOfFirstHitForRecPerChamber[ch] +
1149                  fNHitsForRecPerChamber[ch];
1150                iHit++) {
1151             hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
1152             // condition for hit not already in segment ????
1153             chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
1154             if (chi2 < bestChi2) {
1155               // update best Chi2 and HitForRec if better found
1156               bestHit = hit;
1157               bestChi2 = chi2;
1158               chBestHit = chInStation;
1159             }
1160           }
1161         }
1162         if (bestHit) {
1163           // best hit found: add it to track candidate
1164           track->AddHitForRec(bestHit);
1165           // set track parameters at this TrackHit
1166           track->SetTrackParamAtHit(track->GetNTrackHits() - 1,
1167                                     &(trackParam[chBestHit]));
1168           if (fPrintLevel >= 10) {
1169             cout << "FollowTracks: track candidate(0..): " << trackIndex
1170                  << " Added hit in station(0..): " << station << endl;
1171             track->RecursiveDump();
1172           }
1173         }
1174         else {
1175           // Remove current track candidate
1176           // and corresponding TrackHit's, ...
1177           track->Remove();
1178           delete extrapSegment;
1179           break; // stop the search for this candidate:
1180           // exit from the loop over station
1181         }
1182       }
1183       delete extrapSegment;
1184       // Sort track hits according to increasing Z
1185       track->GetTrackHitsPtr()->Sort();
1186       // Update track parameters at first track hit (smallest Z)
1187       trackParam1 = ((AliMUONTrackHit*)
1188                      (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1189       // Track fit
1190       // with multiple Coulomb scattering if all stations
1191       if (station == 0) track->SetFitMCS(1);
1192       // without multiple Coulomb scattering if not all stations
1193       else track->SetFitMCS(0);
1194       track->SetFitNParam(5);  // with 5 parameters (momentum and position)
1195       track->SetFitStart(1);  // from parameters at first hit
1196       track->Fit();
1197       if (fPrintLevel >= 10) {
1198         cout << "FollowTracks: track candidate(0..): " << trackIndex
1199              << " after fit from station(0..): " << station << " to 4" << endl;
1200         track->RecursiveDump();
1201       }
1202       // Track extrapolation to the vertex through the absorber (Branson)
1203       // after going through the first station
1204       if (station == 0) {
1205         trackParamVertex = *trackParam1;
1206         (&trackParamVertex)->ExtrapToVertex();
1207         track->SetTrackParamAtVertex(&trackParamVertex);
1208         if (fPrintLevel >= 1) {
1209           cout << "FollowTracks: track candidate(0..): " << trackIndex
1210                << " after extrapolation to vertex" << endl;
1211           track->RecursiveDump();
1212         }
1213       }
1214     } // for (station = 2;...
1215     // go really to next track
1216     track = nextTrack;
1217   } // while (track)
1218   // Compression of track array (necessary after Remove ????)
1219   fRecTracksPtr->Compress();
1220   return;
1221 }
1222
1223   //__________________________________________________________________________
1224 void AliMUONEventReconstructor::RemoveDoubleTracks(void)
1225 {
1226   // To remove double tracks.
1227   // Tracks are considered identical
1228   // if they have at least half of their hits in common.
1229   // Among two identical tracks, one keeps the track with the larger number of hits
1230   // or, if these numbers are equal, the track with the minimum Chi2.
1231   AliMUONTrack *track1, *track2, *trackToRemove;
1232   Bool_t identicalTracks;
1233   Int_t hitsInCommon, nHits1, nHits2;
1234   identicalTracks = kTRUE;
1235   while (identicalTracks) {
1236     identicalTracks = kFALSE;
1237     // Loop over first track of the pair
1238     track1 = (AliMUONTrack*) fRecTracksPtr->First();
1239     while (track1 && (!identicalTracks)) {
1240       nHits1 = track1->GetNTrackHits();
1241       // Loop over second track of the pair
1242       track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1243       while (track2 && (!identicalTracks)) {
1244         nHits2 = track2->GetNTrackHits();
1245         // number of hits in common between two tracks
1246         hitsInCommon = track1->HitsInCommon(track2);
1247         // check for identical tracks
1248         if ((4 * hitsInCommon) >= (nHits1 + nHits2)) {
1249           identicalTracks = kTRUE;
1250           // decide which track to remove
1251           if (nHits1 > nHits2) trackToRemove = track2;
1252           else if (nHits1 < nHits2) trackToRemove = track1;
1253           else if ((track1->GetFitFMin()) < (track2->GetFitFMin()))
1254             trackToRemove = track2;
1255           else trackToRemove = track1;
1256           // remove it
1257           trackToRemove->Remove();
1258         }
1259         track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
1260       } // track2
1261       track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1262     } // track1
1263   }
1264   return;
1265 }
1266
1267   //__________________________________________________________________________
1268 void AliMUONEventReconstructor::EventDump(void)
1269 {
1270   // Dump reconstructed event (track parameters at vertex and at first hit),
1271   // and the particle parameters
1272
1273   AliMUONTrack *track;
1274   AliMUONTrackParam *trackParam, *trackParam1;
1275   TClonesArray *particles; // pointer to the particle list
1276   TParticle *p;
1277   Double_t bendingSlope, nonBendingSlope, pYZ;
1278   Double_t pX, pY, pZ, x, y, z, c;
1279   Int_t np, trackIndex, nTrackHits;
1280
1281   if (fPrintLevel >= 1) cout << "****** enter EventDump ******" << endl;
1282   if (fPrintLevel >= 1) {
1283     cout << " Number of Reconstructed tracks :" <<  fNRecTracks << endl;
1284   }
1285   fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1286   // Loop over reconstructed tracks
1287   for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1288     if (fPrintLevel >= 1)
1289       cout << " track number: " << trackIndex << endl;
1290     // function for each track for modularity ????
1291     track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1292     nTrackHits = track->GetNTrackHits();
1293     if (fPrintLevel >= 1)
1294       cout << " number of track hits: " << nTrackHits << endl;
1295     // track parameters at Vertex
1296     trackParam = track->GetTrackParamAtVertex();
1297     x = trackParam->GetNonBendingCoor();
1298     y = trackParam->GetBendingCoor();
1299     z = trackParam->GetZ();
1300     bendingSlope = trackParam->GetBendingSlope();
1301     nonBendingSlope = trackParam->GetNonBendingSlope();
1302     pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1303     pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1304     pX = pZ * nonBendingSlope;
1305     pY = pZ * bendingSlope;
1306     c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1307     if (fPrintLevel >= 1)
1308       printf(" track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1309              z, x, y, pX, pY, pZ, c);
1310
1311     // track parameters at first hit
1312     trackParam1 = ((AliMUONTrackHit*)
1313                    (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1314     x = trackParam1->GetNonBendingCoor();
1315     y = trackParam1->GetBendingCoor();
1316     z = trackParam1->GetZ();
1317     bendingSlope = trackParam1->GetBendingSlope();
1318     nonBendingSlope = trackParam1->GetNonBendingSlope();
1319     pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1320     pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1321     pX = pZ * nonBendingSlope;
1322     pY = pZ * bendingSlope;
1323     c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1324     if (fPrintLevel >= 1)
1325       printf(" track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1326              z, x, y, pX, pY, pZ, c);
1327   }
1328   // informations about generated particles
1329   particles = gAlice->Particles();
1330   np = particles->GetEntriesFast();
1331   printf(" **** number of generated particles: %d  \n", np);
1332
1333   for (Int_t iPart = 0; iPart < np; iPart++) {
1334     p = (TParticle*) particles->UncheckedAt(iPart);
1335     printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1336            iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
1337   }
1338   return;
1339 }
1340