1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.6 2000/06/30 12:01:06 gosset
19 Correction for hit search in the right chamber (JPC)
21 Revision 1.5 2000/06/30 10:15:48 gosset
22 Changes to EventReconstructor...:
23 precision fit with multiple Coulomb scattering;
24 extrapolation to vertex with Branson correction in absorber (JPC)
26 Revision 1.4 2000/06/27 14:11:36 gosset
27 Corrections against violations of coding conventions
29 Revision 1.3 2000/06/16 07:27:08 gosset
30 To remove problem in running RuleChecker, like in MUON-dev
32 Revision 1.1.2.5 2000/06/16 07:00:26 gosset
33 To remove problem in running RuleChecker
35 Revision 1.1.2.4 2000/06/12 08:00:07 morsch
36 Dummy streamer to solve CINT compilation problem (to be investigated !)
38 Revision 1.1.2.3 2000/06/09 20:59:57 morsch
39 Make includes consistent with new file structure.
41 Revision 1.1.2.2 2000/06/09 12:58:05 gosset
42 Removed comment beginnings in Log sections of .cxx files
43 Suppressed most violations of coding rules
45 Revision 1.1.2.1 2000/06/07 14:44:53 gosset
46 Addition of files for track reconstruction in C++
49 //__________________________________________________________________________
51 // MUON event reconstructor in ALICE
53 // This class contains as data:
54 // * the parameters for the event reconstruction
55 // * a pointer to the array of hits to be reconstructed (the event)
56 // * a pointer to the array of segments made with these hits inside each station
57 // * a pointer to the array of reconstructed tracks
59 // It contains as methods, among others:
60 // * MakeEventToBeReconstructed to build the array of hits to be reconstructed
61 // * MakeSegments to build the segments
62 // * MakeTracks to build the tracks
63 //__________________________________________________________________________
70 #include "AliCallf77.h"
71 #include "AliMUONEventReconstructor.h"
73 #include "AliMUONHitForRec.h"
74 #include "AliMUONSegment.h"
75 #include "AliMUONHit.h"
76 #include "AliMUONRawCluster.h"
77 #include "AliMUONTrack.h"
78 #include "AliMUONChamber.h"
79 #include "AliMUONTrackHit.h"
81 #include "TParticle.h"
84 # define initfield initfield_
85 # define reco_gufld reco_gufld_
87 # define initfield INITFIELD
88 # define reco_gufld RECO_GUFLD
93 void type_of_call initfield();
94 void type_of_call reco_gufld(Double_t *Coor, Double_t *Field);
97 //************* Defaults parameters for reconstruction
98 static const Double_t kDefaultMinBendingMomentum = 3.0;
99 static const Double_t kDefaultMaxSigma2Distance = 16.0;
100 static const Double_t kDefaultBendingResolution = 0.01;
101 static const Double_t kDefaultNonBendingResolution = 0.144;
102 static const Double_t kDefaultChamberThicknessInX0 = 0.03;
103 // Simple magnetic field:
104 // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
105 // Length and Position from reco_muon.F, with opposite sign:
106 // Length = ZMAGEND-ZCOIL
107 // Position = (ZMAGEND+ZCOIL)/2
108 // to be ajusted differently from real magnetic field ????
109 static const Double_t kDefaultSimpleBValue = 7.0;
110 static const Double_t kDefaultSimpleBLength = 428.0;
111 static const Double_t kDefaultSimpleBPosition = 1019.0;
112 static const Int_t kDefaultRecGeantHits = 0;
113 static const Double_t kDefaultEfficiency = 0.95;
115 static const Int_t kDefaultPrintLevel = 0;
117 ClassImp(AliMUONEventReconstructor) // Class implementation in ROOT context
119 //__________________________________________________________________________
120 AliMUONEventReconstructor::AliMUONEventReconstructor(void)
122 // Constructor for class AliMUONEventReconstructor
123 SetReconstructionParametersToDefaults();
124 // Memory allocation for the TClonesArray of hits for reconstruction
125 // Is 10000 the right size ????
126 fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
127 fNHitsForRec = 0; // really needed or GetEntriesFast sufficient ????
128 // Memory allocation for the TClonesArray's of segments in stations
129 // Is 2000 the right size ????
130 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
131 fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
132 fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
134 // Memory allocation for the TClonesArray of reconstructed tracks
135 // Is 10 the right size ????
136 fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
137 fNRecTracks = 0; // really needed or GetEntriesFast sufficient ????
139 // Initialize magnetic field
140 // using Fortran subroutine INITFIELD in "reco_muon.F".
141 // Should rather use AliMagF ???? and remove prototyping ...
143 // Impression de quelques valeurs
144 Double_t coor[3], field[3];
148 reco_gufld(coor, field);
149 cout << "coor: " << coor[0] << ", " << coor[1] << ", " << coor[2] << endl;
150 cout << "field: " << field[0] << ", " << field[1] << ", " << field[2] << endl;
152 reco_gufld(coor, field);
153 cout << "coor: " << coor[0] << ", " << coor[1] << ", " << coor[2] << endl;
154 cout << "field: " << field[0] << ", " << field[1] << ", " << field[2] << endl;
157 if (fPrintLevel >= 0) {
158 cout << "AliMUONEventReconstructor constructed with defaults" << endl; Dump();}
162 AliMUONEventReconstructor::AliMUONEventReconstructor (const AliMUONEventReconstructor& Reconstructor)
164 // Dummy copy constructor
167 AliMUONEventReconstructor & AliMUONEventReconstructor::operator=(const AliMUONEventReconstructor& Reconstructor)
169 // Dummy assignment operator
173 //__________________________________________________________________________
174 AliMUONEventReconstructor::~AliMUONEventReconstructor(void)
176 // Destructor for class AliMUONEventReconstructor
177 delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ????
178 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
179 delete fSegmentsPtr[st]; // Correct destruction of everything ????
183 //__________________________________________________________________________
184 void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
186 // Set reconstruction parameters to default values
187 // Would be much more convenient with a structure (or class) ????
188 fMinBendingMomentum = kDefaultMinBendingMomentum;
189 fMaxSigma2Distance = kDefaultMaxSigma2Distance;
191 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
192 // ******** Parameters for making HitsForRec
194 // like in TRACKF_STAT:
195 // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3;
196 // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9
197 for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
198 if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
199 2.0 * TMath::Pi() / 180.0;
200 else fRMin[ch] = 30.0;
203 // like in TRACKF_STAT (10 degrees ????)
204 fRMax[0] = fRMax[1] = 91.5;
205 fRMax[2] = fRMax[3] = 122.5;
206 fRMax[4] = fRMax[5] = 158.3;
207 fRMax[6] = fRMax[7] = 260.0;
208 fRMax[8] = fRMax[9] = 260.0;
210 // ******** Parameters for making segments
211 // should be parametrized ????
212 // according to interval between chambers in a station ????
213 // Maximum distance in non bending plane
214 // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
216 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
217 fSegmentMaxDistNonBending[st] = 5. * 0.22;
218 // Maximum distance in bending plane
219 // values from TRACKF_STAT corresponding to (J psi 20cm)
220 fSegmentMaxDistBending[0] = 1.5;
221 fSegmentMaxDistBending[1] = 1.5;
222 fSegmentMaxDistBending[2] = 3.0;
223 fSegmentMaxDistBending[3] = 6.0;
224 fSegmentMaxDistBending[4] = 6.0;
226 fBendingResolution = kDefaultBendingResolution;
227 fNonBendingResolution = kDefaultNonBendingResolution;
228 fChamberThicknessInX0 = kDefaultChamberThicknessInX0;
229 fSimpleBValue = kDefaultSimpleBValue;
230 fSimpleBLength = kDefaultSimpleBLength;
231 fSimpleBPosition = kDefaultSimpleBPosition;
232 fRecGeantHits = kDefaultRecGeantHits;
233 fEfficiency = kDefaultEfficiency;
234 fPrintLevel = kDefaultPrintLevel;
238 //__________________________________________________________________________
239 Double_t AliMUONEventReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum)
241 // Returns impact parameter at vertex in bending plane (cm),
242 // from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
243 // using simple values for dipole magnetic field.
244 // The sign is the sign of the charge.
245 return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
249 //__________________________________________________________________________
250 Double_t AliMUONEventReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam)
252 // Returns signed bending momentum in bending plane (GeV/c),
253 // from the impact parameter "ImpactParam" at vertex in bending plane (cm),
254 // using simple values for dipole magnetic field.
255 // The sign is the sign of the charge.
256 return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
260 //__________________________________________________________________________
261 void AliMUONEventReconstructor::SetBkgGeantFile(Text_t *BkgGeantFileName)
263 // Set background file ... for GEANT hits
264 // Must be called after having loaded the firts signal event
265 if (fPrintLevel >= 0) {
266 cout << "Enter SetBkgGeantFile with BkgGeantFileName ``"
267 << BkgGeantFileName << "''" << endl;}
268 if (strlen(BkgGeantFileName)) {
269 // BkgGeantFileName not empty: try to open the file
270 if (fPrintLevel >= 2) {cout << "Before File(Bkg)" << endl; gDirectory->Dump();}
271 fBkgGeantFile = new TFile(BkgGeantFileName);
272 if (fPrintLevel >= 2) {cout << "After File(Bkg)" << endl; gDirectory->Dump();}
273 if (fBkgGeantFile-> IsOpen()) {
274 if (fPrintLevel >= 0) {
275 cout << "Background for GEANT hits in file: ``" << BkgGeantFileName
276 << "'' successfully opened" << endl;}
279 cout << "Background for GEANT hits in file: " << BkgGeantFileName << endl;
280 cout << "NOT FOUND: EXIT" << endl;
281 exit(0); // right instruction for exit ????
283 // Arrays for "particles" and "hits"
284 fBkgGeantParticles = new TClonesArray("TParticle", 200);
285 fBkgGeantHits = new TClonesArray("AliMUONHit", 2000);
286 // Event number to -1 for initialization
287 fBkgGeantEventNumber = -1;
288 // Back to the signal file:
289 // first signal event must have been loaded previously,
290 // otherwise, Segmentation violation at the next instruction
291 // How is it possible to do smething better ????
292 ((gAlice->TreeK())->GetCurrentFile())->cd();
293 if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();}
298 //__________________________________________________________________________
299 void AliMUONEventReconstructor::NextBkgGeantEvent(void)
301 // Get next event in background file for GEANT hits
302 // Goes back to event number 0 when end of file is reached
305 if (fPrintLevel >= 0) {
306 cout << "Enter NextBkgGeantEvent" << endl;}
307 // Clean previous event
308 if(fBkgGeantTK) delete fBkgGeantTK;
310 if(fBkgGeantParticles) fBkgGeantParticles->Clear();
311 if(fBkgGeantTH) delete fBkgGeantTH;
313 if(fBkgGeantHits) fBkgGeantHits->Clear();
314 // Increment event number
315 fBkgGeantEventNumber++;
316 // Get access to Particles and Hits for event from background file
317 if (fPrintLevel >= 2) {cout << "Before cd(Bkg)" << endl; gDirectory->Dump();}
319 if (fPrintLevel >= 2) {cout << "After cd(Bkg)" << endl; gDirectory->Dump();}
320 // Particles: TreeK for event and branch "Particles"
321 sprintf(treeName, "TreeK%d", fBkgGeantEventNumber);
322 fBkgGeantTK = (TTree*)gDirectory->Get(treeName);
324 if (fPrintLevel >= 0) {
325 cout << "Cannot find Kine Tree for background event: " <<
326 fBkgGeantEventNumber << endl;
327 cout << "Goes back to event 0" << endl;
329 fBkgGeantEventNumber = 0;
330 sprintf(treeName, "TreeK%d", fBkgGeantEventNumber);
331 fBkgGeantTK = (TTree*)gDirectory->Get(treeName);
333 cout << "ERROR: cannot find Kine Tree for background event: " <<
334 fBkgGeantEventNumber << endl;
339 fBkgGeantTK->SetBranchAddress("Particles", &fBkgGeantParticles);
340 fBkgGeantTK->GetEvent(0); // why event 0 ???? necessary ????
341 // Hits: TreeH for event and branch "MUON"
342 sprintf(treeName, "TreeH%d", fBkgGeantEventNumber);
343 fBkgGeantTH = (TTree*)gDirectory->Get(treeName);
345 cout << "ERROR: cannot find Hits Tree for background event: " <<
346 fBkgGeantEventNumber << endl;
349 if (fBkgGeantTH && fBkgGeantHits) {
350 branch = fBkgGeantTH->GetBranch("MUON");
351 if (branch) branch->SetAddress(&fBkgGeantHits);
353 fBkgGeantTH->GetEntries(); // necessary ????
354 // Back to the signal file
355 ((gAlice->TreeK())->GetCurrentFile())->cd();
356 if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();}
360 //__________________________________________________________________________
361 void AliMUONEventReconstructor::EventReconstruct(void)
363 // To reconstruct one event
364 if (fPrintLevel >= 1) cout << "enter EventReconstruct" << endl;
365 MakeEventToBeReconstructed();
371 //__________________________________________________________________________
372 void AliMUONEventReconstructor::ResetHitsForRec(void)
374 // To reset the array and the number of HitsForRec,
375 // and also the number of HitsForRec
376 // and the index of the first HitForRec per chamber
377 if (fHitsForRecPtr) fHitsForRecPtr->Clear();
379 for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++)
380 fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
384 //__________________________________________________________________________
385 void AliMUONEventReconstructor::ResetSegments(void)
387 // To reset the TClonesArray of segments and the number of Segments
389 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
390 if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
396 //__________________________________________________________________________
397 void AliMUONEventReconstructor::ResetTracks(void)
399 // To reset the TClonesArray of reconstructed tracks
400 if (fRecTracksPtr) fRecTracksPtr->Clear();
405 //__________________________________________________________________________
406 void AliMUONEventReconstructor::MakeEventToBeReconstructed(void)
408 // To make the list of hits to be reconstructed,
409 // either from the GEANT hits or from the raw clusters
410 // according to the parameter set for the reconstructor
411 if (fPrintLevel >= 1) cout << "enter MakeEventToBeReconstructed" << endl;
413 if (fRecGeantHits == 1) {
414 // Reconstruction from GEANT hits
415 // Back to the signal file
416 ((gAlice->TreeK())->GetCurrentFile())->cd();
418 // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
419 // Security on MUON ????
420 AddHitsForRecFromGEANT(gAlice->TreeH());
422 AddHitsForRecFromBkgGEANT(fBkgGeantTH, fBkgGeantHits);
423 // Sort HitsForRec in increasing order with respect to chamber number
424 SortHitsForRecWithIncreasingChamber();
427 // Reconstruction from raw clusters
428 // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
429 // Security on MUON ????
430 // TreeR assumed to be be "prepared" in calling function
431 // by "MUON->GetTreeR(nev)" ????
432 TTree *treeR = gAlice->TreeR();
433 AddHitsForRecFromRawClusters(treeR);
434 // No sorting: it is done automatically in the previous function
436 if (fPrintLevel >= 10) {
437 cout << "end of MakeEventToBeReconstructed" << endl;
438 cout << "NHitsForRec: " << fNHitsForRec << endl;
439 for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
440 cout << "chamber(0...): " << ch
441 << " NHitsForRec: " << fNHitsForRecPerChamber[ch]
442 << " index(first HitForRec): " << fIndexOfFirstHitForRecPerChamber[ch]
444 for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
445 hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
447 cout << "HitForRec index(0...): " << hit << endl;
448 ((*fHitsForRecPtr)[hit])->Dump();
455 //__________________________________________________________________________
456 void AliMUONEventReconstructor::AddHitsForRecFromGEANT(TTree *TH)
458 // To add to the list of hits for reconstruction
459 // the GEANT signal hits from a hit tree TH.
460 if (fPrintLevel >= 2)
461 cout << "enter AddHitsForRecFromGEANT with TH: " << TH << endl;
462 if (TH == NULL) return;
463 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
464 // Security on MUON ????
465 // See whether it could be the same for signal and background ????
466 // Loop over tracks in tree
467 Int_t ntracks = (Int_t) TH->GetEntries();
468 if (fPrintLevel >= 2)
469 cout << "ntracks: " << ntracks << endl;
470 for (Int_t track = 0; track < ntracks; track++) {
475 for (AliMUONHit* mHit = (AliMUONHit*) pMUON->FirstHit(-1);
477 mHit = (AliMUONHit*) pMUON->NextHit(), hit++) {
478 NewHitForRecFromGEANT(mHit,track, hit, 1);
480 } // end of track loop
484 //__________________________________________________________________________
485 void AliMUONEventReconstructor::AddHitsForRecFromBkgGEANT(TTree *TH, TClonesArray *Hits)
487 // To add to the list of hits for reconstruction
488 // the GEANT background hits from a hit tree TH and a pointer Hits to a hit list.
489 // How to have only one function "AddHitsForRecFromGEANT" ????
490 if (fPrintLevel >= 2)
491 cout << "enter AddHitsForRecFromBkgGEANT with TH: " << TH << endl;
492 if (TH == NULL) return;
493 // Loop over tracks in tree
494 Int_t ntracks = (Int_t) TH->GetEntries();
495 if (fPrintLevel >= 2)
496 cout << "ntracks: " << ntracks << endl;
497 for (Int_t track = 0; track < ntracks; track++) {
498 if (Hits) Hits->Clear();
501 for (Int_t hit = 0; hit < Hits->GetEntriesFast(); hit++) {
502 NewHitForRecFromGEANT((AliMUONHit*) (*Hits)[hit], track, hit, 0);
504 } // end of track loop
508 //__________________________________________________________________________
509 AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* Hit, Int_t TrackNumber, Int_t HitNumber, Int_t Signal)
511 // To make a new hit for reconstruction from a GEANT hit pointed to by "Hit",
512 // with hit number "HitNumber" in the track numbered "TrackNumber",
513 // either from signal ("Signal" = 1) or background ("Signal" = 0) event.
514 // Selects hits in tracking (not trigger) chambers.
515 // Takes into account the efficiency (fEfficiency)
516 // and the smearing from resolution (fBendingResolution and fNonBendingResolution).
517 // Adds a condition on the radius between RMin and RMax
518 // to better simulate the real chambers.
519 // Returns the pointer to the new hit for reconstruction,
520 // or NULL in case of inefficiency or non tracking chamber or bad radius.
521 // No condition on at most 20 cm from a muon from a resonance
522 // like in Fortran TRACKF_STAT.
523 AliMUONHitForRec* hitForRec;
524 Double_t bendCoor, nonBendCoor, radius;
525 Int_t chamber = Hit->fChamber - 1; // chamber(0...)
526 // only in tracking chambers (fChamber starts at 1)
527 if (chamber >= kMaxMuonTrackingChambers) return NULL;
528 // only if hit is efficient (keep track for checking ????)
529 if (gRandom->Rndm() > fEfficiency) return NULL;
530 // only if radius between RMin and RMax
532 nonBendCoor = Hit->fX;
533 radius = TMath::Sqrt((bendCoor * bendCoor) + (nonBendCoor * nonBendCoor));
534 if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL;
535 // new AliMUONHitForRec from GEANT hit and increment number of AliMUONHitForRec's
536 hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(Hit);
538 // add smearing from resolution
539 hitForRec->SetBendingCoor(bendCoor + gRandom->Gaus(0., fBendingResolution));
540 hitForRec->SetNonBendingCoor(nonBendCoor
541 + gRandom->Gaus(0., fNonBendingResolution));
542 // more information into HitForRec
543 // resolution: angular effect to be added here ????
544 hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
545 hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
547 hitForRec->SetHitNumber(HitNumber);
548 hitForRec->SetTHTrack(TrackNumber);
549 hitForRec->SetGeantSignal(Signal);
550 if (fPrintLevel >= 10) {
551 cout << "track: " << TrackNumber << " hit: " << HitNumber << endl;
553 cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
558 //__________________________________________________________________________
559 void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber()
561 // Sort HitsForRec's in increasing order with respect to chamber number.
562 // Uses the function "Compare".
563 // Update the information for HitsForRec per chamber too.
564 Int_t ch, nhits, prevch;
565 fHitsForRecPtr->Sort();
566 for (ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
567 fNHitsForRecPerChamber[ch] = 0;
568 fIndexOfFirstHitForRecPerChamber[ch] = 0;
570 prevch = 0; // previous chamber
571 nhits = 0; // number of hits in current chamber
572 // Loop over HitsForRec
573 for (Int_t hit = 0; hit < fNHitsForRec; hit++) {
574 // chamber number (0...)
575 ch = ((AliMUONHitForRec*) ((*fHitsForRecPtr)[hit]))->GetChamberNumber();
576 // increment number of hits in current chamber
577 (fNHitsForRecPerChamber[ch])++;
578 // update index of first HitForRec in current chamber
579 // if chamber number different from previous one
581 fIndexOfFirstHitForRecPerChamber[ch] = hit;
588 // //__________________________________________________________________________
589 // void AliMUONEventReconstructor::AddHitsForRecFromCathodeCorrelations(TTree* TC)
591 // // OLD VERSION WHEN ONE ONE WAS USING SO CALLED CATHODE CORRELATIONS
592 // // To add to the list of hits for reconstruction
593 // // the (cathode correlated) raw clusters
594 // // No condition added, like in Fortran TRACKF_STAT,
595 // // on the radius between RMin and RMax.
596 // AliMUONHitForRec *hitForRec;
597 // if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromCathodeCorrelations" << endl;
598 // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
599 // // Security on MUON ????
600 // // Loop over tracking chambers
601 // for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
602 // // number of HitsForRec to 0 for the chamber
603 // fNHitsForRecPerChamber[ch] = 0;
604 // // index of first HitForRec for the chamber
605 // if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
606 // else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
607 // TClonesArray *reconst_hits = MUON->ReconstHitsAddress(ch);
608 // MUON->ResetReconstHits();
610 // Int_t ncor = (Int_t)reconst_hits->GetEntries();
611 // // Loop over (cathode correlated) raw clusters
612 // for (Int_t cor = 0; cor < ncor; cor++) {
613 // AliMUONReconstHit * mCor =
614 // (AliMUONReconstHit*) reconst_hits->UncheckedAt(cor);
615 // // new AliMUONHitForRec from (cathode correlated) raw cluster
616 // // and increment number of AliMUONHitForRec's (total and in chamber)
617 // hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(mCor);
619 // (fNHitsForRecPerChamber[ch])++;
620 // // more information into HitForRec
621 // hitForRec->SetChamberNumber(ch);
622 // hitForRec->SetHitNumber(cor);
623 // // Z coordinate of the chamber (cm) with sign opposite to GEANT sign
624 // // could (should) be more exact from chamber geometry ????
625 // hitForRec->SetZ(-(&(MUON->Chamber(ch)))->Z());
626 // if (fPrintLevel >= 10) {
627 // cout << "chamber (0...): " << ch <<
628 // " cathcorrel (0...): " << cor << endl;
630 // cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
631 // hitForRec->Dump();}
632 // } // end of cluster loop
633 // } // end of chamber loop
637 //__________________________________________________________________________
638 void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
640 // To add to the list of hits for reconstruction all the raw clusters
641 // No condition added, like in Fortran TRACKF_STAT,
642 // on the radius between RMin and RMax.
643 AliMUONHitForRec *hitForRec;
644 AliMUONRawCluster *clus;
646 TClonesArray *rawclusters;
647 if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromRawClusters" << endl;
648 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
649 // Security on MUON ????
650 // Loop over tracking chambers
651 for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
652 // number of HitsForRec to 0 for the chamber
653 fNHitsForRecPerChamber[ch] = 0;
654 // index of first HitForRec for the chamber
655 if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
656 else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
657 rawclusters = pMUON->RawClustAddress(ch);
658 pMUON->ResetRawClusters();
659 TR->GetEvent((Int_t) (TR->GetEntries()) - 1); // to be checked ????
660 nclus = (Int_t) (rawclusters->GetEntries());
661 // Loop over (cathode correlated) raw clusters
662 for (iclus = 0; iclus < nclus; iclus++) {
663 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
664 // new AliMUONHitForRec from raw cluster
665 // and increment number of AliMUONHitForRec's (total and in chamber)
666 hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
668 (fNHitsForRecPerChamber[ch])++;
669 // more information into HitForRec
670 // resolution: info should be already in raw cluster and taken from it ????
671 hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
672 hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
673 // original raw cluster
674 hitForRec->SetChamberNumber(ch);
675 hitForRec->SetHitNumber(iclus);
676 // Z coordinate of the chamber (cm) with sign opposite to GEANT sign
677 // could (should) be more exact from chamber geometry ????
678 hitForRec->SetZ(-(&(pMUON->Chamber(ch)))->Z());
679 if (fPrintLevel >= 10) {
680 cout << "chamber (0...): " << ch <<
681 " raw cluster (0...): " << iclus << endl;
683 cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
685 } // end of cluster loop
686 } // end of chamber loop
690 //__________________________________________________________________________
691 void AliMUONEventReconstructor::MakeSegments(void)
693 // To make the list of segments in all stations,
694 // from the list of hits to be reconstructed
695 if (fPrintLevel >= 1) cout << "enter MakeSegments" << endl;
697 // Loop over stations
698 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
699 MakeSegmentsPerStation(st);
700 if (fPrintLevel >= 10) {
701 cout << "end of MakeSegments" << endl;
702 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
703 cout << "station(0...): " << st
704 << " Segments: " << fNSegments[st]
707 seg < fNSegments[st];
709 cout << "Segment index(0...): " << seg << endl;
710 ((*fSegmentsPtr[st])[seg])->Dump();
717 //__________________________________________________________________________
718 void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station)
720 // To make the list of segments in station number "Station" (0...)
721 // from the list of hits to be reconstructed.
722 // Updates "fNSegments"[Station].
723 // Segments in stations 4 and 5 are sorted
724 // according to increasing absolute value of "impact parameter"
725 AliMUONHitForRec *hit1Ptr, *hit2Ptr;
726 AliMUONSegment *segment;
728 Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
729 impactParam, maxImpactParam;
730 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
731 if (fPrintLevel >= 1)
732 cout << "enter MakeSegmentsPerStation (0...) " << Station << endl;
733 // first and second chambers (0...) in the station
734 Int_t ch1 = 2 * Station;
736 // variable true for stations downstream of the dipole:
737 // Station(0..4) equal to 3 or 4
738 if ((Station == 3) || (Station == 4)) {
740 // maximum impact parameter (cm) according to fMinBendingMomentum...
742 TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
744 else last2st = kFALSE;
745 // extrapolation factor from Z of first chamber to Z of second chamber
746 // dZ to be changed to take into account fine structure of chambers ????
747 Double_t extrapFact =
748 (&(pMUON->Chamber(ch2)))->Z() / (&(pMUON->Chamber(ch1)))->Z();
749 // index for current segment
750 Int_t segmentIndex = 0;
751 // Loop over HitsForRec in the first chamber of the station
752 for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
753 hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
755 // pointer to the HitForRec
756 hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
758 // on the straight line joining the HitForRec to the vertex (0,0,0),
759 // to the Z of the second chamber of the station
760 extBendCoor = extrapFact * hit1Ptr->GetBendingCoor();
761 extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor();
762 // Loop over HitsForRec in the second chamber of the station
763 for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
764 hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
766 // pointer to the HitForRec
767 hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]);
768 // absolute values of distances, in bending and non bending planes,
769 // between the HitForRec in the second chamber
770 // and the previous extrapolation
771 distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor);
772 distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor);
775 bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) /
776 (hit1Ptr->GetZ() - hit2Ptr->GetZ());
777 // absolute value of impact parameter
779 TMath::Abs(hit1Ptr->GetBendingCoor() - hit2Ptr->GetZ() * bendingSlope);
781 // check for distances not too large,
782 // and impact parameter not too big if stations downstream of the dipole.
783 // Conditions "distBend" and "impactParam" correlated for these stations ????
784 if ((distBend < fSegmentMaxDistBending[Station]) &&
785 (distNonBend < fSegmentMaxDistNonBending[Station]) &&
786 (!last2st || (impactParam < maxImpactParam))) {
788 segment = new ((*fSegmentsPtr[Station])[segmentIndex])
789 AliMUONSegment(hit1Ptr, hit2Ptr);
790 // update "link" to this segment from the hit in the first chamber
791 if (hit1Ptr->GetNSegments() == 0)
792 hit1Ptr->SetIndexOfFirstSegment(segmentIndex);
793 hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1);
794 if (fPrintLevel >= 10) {
795 cout << "segmentIndex(0...): " << segmentIndex
796 << " distBend: " << distBend
797 << " distNonBend: " << distNonBend
800 cout << "HitForRec in first chamber" << endl;
802 cout << "HitForRec in second chamber" << endl;
805 // increment index for current segment
809 } // for (Int_t hit1...
810 fNSegments[Station] = segmentIndex;
811 // Sorting according to "impact parameter" if station(1..5) 4 or 5,
812 // i.e. Station(0..4) 3 or 4, using the function "Compare".
813 // After this sorting, it is impossible to use
814 // the "fNSegments" and "fIndexOfFirstSegment"
815 // of the HitForRec in the first chamber to explore all segments formed with it.
816 // Is this sorting really needed ????
817 if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort();
818 if (fPrintLevel >= 1) cout << "Station: " << Station << " NSegments: "
819 << fNSegments[Station] << endl;
823 //__________________________________________________________________________
824 void AliMUONEventReconstructor::MakeTracks(void)
826 // To make the tracks,
827 // from the list of segments and points in all stations
828 if (fPrintLevel >= 1) cout << "enter MakeTracks" << endl;
830 // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
831 MakeTrackCandidates();
832 // Follow tracks in stations(1..) 3, 2 and 1
837 //__________________________________________________________________________
838 Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
840 // To make track candidates with two segments in stations(1..) 4 and 5,
841 // the first segment being pointed to by "BegSegment".
842 // Returns the number of such track candidates.
843 Int_t endStation, iEndSegment, nbCan2Seg;
844 AliMUONSegment *endSegment, *extrapSegment;
845 AliMUONTrack *recTrack;
847 if (fPrintLevel >= 1) cout << "enter MakeTrackCandidatesWithTwoSegments" << endl;
848 // Station for the end segment
849 endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
850 // multiple scattering factor corresponding to one chamber
852 GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
853 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
854 // linear extrapolation to end station
856 BegSegment->CreateSegmentFromLinearExtrapToStation(endStation, mcsFactor);
857 // number of candidates with 2 segments to 0
859 // Loop over segments in the end station
860 for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
861 // pointer to segment
862 endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
863 // test compatibility between current segment and "extrapSegment"
864 // 4 because 4 quantities in chi2
866 NormalizedChi2WithSegment(extrapSegment,
867 fMaxSigma2Distance)) <= 4.0) {
868 // both segments compatible:
869 // make track candidate from "begSegment" and "endSegment"
870 if (fPrintLevel >= 2)
871 cout << "TrackCandidate with Segment " << iEndSegment <<
872 " in Station(0..) " << endStation << endl;
873 // flag for both segments in one track:
874 // to be done in track constructor ????
875 BegSegment->SetInTrack(kTRUE);
876 endSegment->SetInTrack(kTRUE);
877 recTrack = new ((*fRecTracksPtr)[fNRecTracks])
878 AliMUONTrack(BegSegment, endSegment, this);
880 if (fPrintLevel >= 10) recTrack->RecursiveDump();
881 // increment number of track candidates with 2 segments
884 } // for (iEndSegment = 0;...
885 delete extrapSegment; // should not delete HitForRec's it points to !!!!
889 //__________________________________________________________________________
890 Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
892 // To make track candidates with one segment and one point
893 // in stations(1..) 4 and 5,
894 // the segment being pointed to by "BegSegment".
895 Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
896 AliMUONHitForRec *extrapHitForRec, *hit;
897 AliMUONTrack *recTrack;
899 if (fPrintLevel >= 1)
900 cout << "enter MakeTrackCandidatesWithOneSegmentAndOnePoint" << endl;
901 // station for the end point
902 endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
903 // multiple scattering factor corresponding to one chamber
905 GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
906 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
907 // first and second chambers(0..) in the end station
908 ch1 = 2 * endStation;
910 // number of candidates to 0
912 // Loop over chambers of the end station
913 for (ch = ch2; ch >= ch1; ch--) {
914 // linear extrapolation to chamber
916 BegSegment->CreateHitForRecFromLinearExtrapToChamber(ch, mcsFactor);
917 // limits for the hit index in the loop
918 iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
919 iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
920 // Loop over HitForRec's in the chamber
921 for (iHit = iHitMin; iHit < iHitMax; iHit++) {
922 // pointer to HitForRec
923 hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
924 // test compatibility between current HitForRec and "extrapHitForRec"
925 // 2 because 2 quantities in chi2
927 NormalizedChi2WithHitForRec(extrapHitForRec,
928 fMaxSigma2Distance)) <= 2.0) {
929 // both HitForRec's compatible:
930 // make track candidate from begSegment and current HitForRec
931 if (fPrintLevel >= 2)
932 cout << "TrackCandidate with HitForRec " << iHit <<
933 " in Chamber(0..) " << ch << endl;
934 // flag for beginning segments in one track:
935 // to be done in track constructor ????
936 BegSegment->SetInTrack(kTRUE);
937 recTrack = new ((*fRecTracksPtr)[fNRecTracks])
938 AliMUONTrack(BegSegment, hit, this);
939 // the right place to eliminate "double counting" ???? how ????
941 if (fPrintLevel >= 10) recTrack->RecursiveDump();
942 // increment number of track candidates
945 } // for (iHit = iHitMin;...
946 delete extrapHitForRec;
947 } // for (ch = ch2;...
948 return nbCan1Seg1Hit;
951 //__________________________________________________________________________
952 void AliMUONEventReconstructor::MakeTrackCandidates(void)
954 // To make track candidates
955 // with at least 3 aligned points in stations(1..) 4 and 5
956 // (two Segment's or one Segment and one HitForRec)
957 Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
958 AliMUONSegment *begSegment;
959 if (fPrintLevel >= 1) cout << "enter MakeTrackCandidates" << endl;
960 // Loop over stations(1..) 5 and 4 for the beginning segment
961 for (begStation = 4; begStation > 2; begStation--) {
962 // Loop over segments in the beginning station
963 for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
964 // pointer to segment
965 begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
966 if (fPrintLevel >= 2)
967 cout << "look for TrackCandidate's with Segment " << iBegSegment <<
968 " in Station(0..) " << begStation << endl;
969 // Look for track candidates with two segments,
970 // "begSegment" and all compatible segments in other station.
971 // Only for beginning station(1..) 5
972 // because candidates with 2 segments have to looked for only once.
974 nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
975 // Look for track candidates with one segments and one point,
976 // "begSegment" and all compatible HitForRec's in other station.
977 // Only if "begSegment" does not belong already to a track candidate.
978 // Is that a too strong condition ????
979 if (!(begSegment->GetInTrack()))
980 nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
981 } // for (iBegSegment = 0;...
982 } // for (begStation = 4;...
986 //__________________________________________________________________________
987 void AliMUONEventReconstructor::FollowTracks(void)
989 // Follow tracks in stations(1..) 3, 2 and 1
990 // too long: should be made more modular !!!!
991 AliMUONHitForRec *bestHit, *extrapHit, *hit;
992 AliMUONSegment *bestSegment, *extrapSegment, *segment;
993 AliMUONTrack *track, *nextTrack;
994 AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
995 Int_t ch, chInStation, chBestHit, iHit, iSegment, station, trackIndex;
996 Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
997 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
998 // local maxSigma2Distance, for easy increase in testing
999 maxSigma2Distance = fMaxSigma2Distance;
1000 if (fPrintLevel >= 2)
1001 cout << "enter FollowTracks" << endl;
1002 // Loop over track candidates
1003 track = (AliMUONTrack*) fRecTracksPtr->First();
1006 // Follow function for each track candidate ????
1008 nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
1009 if (fPrintLevel >= 2)
1010 cout << "FollowTracks: track candidate(0..): " << trackIndex << endl;
1011 // Fit track candidate from vertex at X = Y = 0
1012 track->SetFitMCS(0); // fit without Multiple Scattering
1013 track->Fit(track->GetTrackParamAtVertex(), 3);
1014 if (fPrintLevel >= 10) {
1015 cout << "FollowTracks: track candidate(0..): " << trackIndex
1016 << " after fit in stations(1..) 4 and 5" << endl;
1017 track->RecursiveDump();
1019 // Loop over stations(1..) 3, 2 and 1
1020 // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
1021 // otherwise: majority coincidence 2 !!!!
1022 for (station = 2; station >= 0; station--) {
1023 // Track parameters at first track hit (smallest Z)
1024 trackParam1 = ((AliMUONTrackHit*)
1025 (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1026 // extrapolation to station
1027 trackParam1->ExtrapToStation(station, trackParam);
1028 extrapSegment = new AliMUONSegment(); // empty segment
1029 // multiple scattering factor corresponding to one chamber
1030 // and momentum in bending plane (not total)
1031 mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
1032 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
1033 // Z difference from previous station
1034 dZ1 = (&(pMUON->Chamber(2 * station)))->Z() -
1035 (&(pMUON->Chamber(2 * station + 2)))->Z();
1036 // Z difference between the two previous stations
1037 dZ2 = (&(pMUON->Chamber(2 * station + 2)))->Z() -
1038 (&(pMUON->Chamber(2 * station + 4)))->Z();
1039 // Z difference between the two chambers in the previous station
1040 dZ3 = (&(pMUON->Chamber(2 * station)))->Z() -
1041 (&(pMUON->Chamber(2 * station + 1)))->Z();
1042 extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
1044 SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
1045 extrapSegment->UpdateFromStationTrackParam
1046 (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
1047 trackParam1->GetInverseBendingMomentum());
1050 if (fPrintLevel >= 10) {
1051 cout << "FollowTracks: track candidate(0..): " << trackIndex
1052 << " Look for segment in station(0..): " << station << endl;
1054 // Loop over segments in station
1055 for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
1056 // Look for best compatible Segment in station
1057 // should consider all possibilities ????
1058 // multiple scattering ????
1059 // separation in 2 functions: Segment and HitForRec ????
1060 segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
1061 chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
1062 if (chi2 < bestChi2) {
1063 // update best Chi2 and Segment if better found
1064 bestSegment = segment;
1069 // best segment found: add it to track candidate
1070 track->AddSegment(bestSegment);
1071 // set track parameters at these two TrakHit's
1072 track->SetTrackParamAtHit(track->GetNTrackHits() - 2, &(trackParam[0]));
1073 track->SetTrackParamAtHit(track->GetNTrackHits() - 1, &(trackParam[1]));
1074 if (fPrintLevel >= 10) {
1075 cout << "FollowTracks: track candidate(0..): " << trackIndex
1076 << " Added segment in station(0..): " << station << endl;
1077 track->RecursiveDump();
1081 // No best segment found:
1082 // Look for best compatible HitForRec in station:
1083 // should consider all possibilities ????
1084 // multiple scattering ???? do about like for extrapSegment !!!!
1085 extrapHit = new AliMUONHitForRec(); // empty hit
1088 if (fPrintLevel >= 10) {
1089 cout << "FollowTracks: track candidate(0..): " << trackIndex
1090 << " Segment not found, look for hit in station(0..): " << station
1093 // Loop over chambers of the station
1094 for (chInStation = 0; chInStation < 2; chInStation++) {
1095 // coordinates of extrapolated hit
1097 SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
1099 SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
1100 // resolutions from "extrapSegment"
1101 extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
1102 extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
1103 // Loop over hits in the chamber
1104 ch = 2 * station + ch;
1105 for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
1106 iHit < fIndexOfFirstHitForRecPerChamber[ch] +
1107 fNHitsForRecPerChamber[ch];
1109 hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
1110 // condition for hit not already in segment ????
1111 chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
1112 if (chi2 < bestChi2) {
1113 // update best Chi2 and HitForRec if better found
1116 chBestHit = chInStation;
1121 // best hit found: add it to track candidate
1122 track->AddHitForRec(bestHit);
1123 // set track parameters at these two TrackHit's
1124 track->SetTrackParamAtHit(track->GetNTrackHits() - 1,
1125 &(trackParam[chBestHit]));
1126 if (fPrintLevel >= 10) {
1127 cout << "FollowTracks: track candidate(0..): " << trackIndex
1128 << " Added hit in station(0..): " << station << endl;
1129 track->RecursiveDump();
1133 // Remove current track candidate and update fNRecTracks
1134 // To be checked: recursive delete of TrackHit's !!!!
1135 // For cleaner implementation: call track->Remove()
1136 // to be coded, with all cleanings,
1137 // including links between HitForRec's and TrackHit's !!!!
1138 fRecTracksPtr->Remove(track);
1140 delete extrapSegment;
1141 break; // stop the search for this candidate:
1142 // exit from the loop over station
1145 delete extrapSegment;
1146 // Sort track hits according to increasing Z
1147 track->GetTrackHitsPtr()->Sort();
1148 // Update track parameters at first track hit (smallest Z)
1149 trackParam1 = ((AliMUONTrackHit*)
1150 (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1151 // Track fit from first track hit varying X and Y,
1152 // with multiple Coulomb scattering if all stations
1153 if (station == 0) track->SetFitMCS(1);
1154 else track->SetFitMCS(0);
1155 track->Fit(trackParam1, 5);
1156 if (fPrintLevel >= 10) {
1157 cout << "FollowTracks: track candidate(0..): " << trackIndex
1158 << " after fit from station(0..): " << station << " to 4" << endl;
1159 track->RecursiveDump();
1161 // Track extrapolation to the vertex through the absorber (Branson)
1163 trackParamVertex = *trackParam1;
1164 (&trackParamVertex)->ExtrapToVertex();
1165 track->SetTrackParamAtVertex(&trackParamVertex);
1166 if (fPrintLevel >= 1) {
1167 cout << "FollowTracks: track candidate(0..): " << trackIndex
1168 << " after extrapolation to vertex" << endl;
1169 track->RecursiveDump();
1172 } // for (station = 2;...
1173 // go really to next track
1176 // Compression of track array (necessary after Remove ????)
1177 fRecTracksPtr->Compress();
1181 //__________________________________________________________________________
1182 void AliMUONEventReconstructor::EventDump(void)
1184 // Dump reconstructed event (track parameters at vertex and at first hit),
1185 // and the particle parameters
1187 AliMUONTrack *track;
1188 AliMUONTrackParam *trackParam, *trackParam1;
1189 TClonesArray *particles; // pointer to the particle list
1191 Double_t bendingSlope, nonBendingSlope, pYZ;
1192 Double_t pX, pY, pZ, x, y, z, c;
1193 Int_t np, trackIndex, nTrackHits;
1195 if (fPrintLevel >= 1) cout << "****** enter EventDump ******" << endl;
1196 if (fPrintLevel >= 1) {
1197 cout << " Number of Reconstructed tracks :" << fNRecTracks << endl;
1199 fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1200 // Loop over reconstructed tracks
1201 for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1202 if (fPrintLevel >= 1)
1203 cout << " track number: " << trackIndex << endl;
1204 // function for each track for modularity ????
1205 track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1206 nTrackHits = track->GetNTrackHits();
1207 if (fPrintLevel >= 1)
1208 cout << " number of track hits: " << nTrackHits << endl;
1209 // track parameters at Vertex
1210 trackParam = track->GetTrackParamAtVertex();
1211 x = trackParam->GetNonBendingCoor();
1212 y = trackParam->GetBendingCoor();
1213 z = trackParam->GetZ();
1214 bendingSlope = trackParam->GetBendingSlope();
1215 nonBendingSlope = trackParam->GetNonBendingSlope();
1216 pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1217 pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1218 pX = pZ * nonBendingSlope;
1219 pY = pZ * bendingSlope;
1220 c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1221 if (fPrintLevel >= 1)
1222 printf(" track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1223 z, x, y, pX, pY, pZ, c);
1225 // track parameters at first hit
1226 trackParam1 = ((AliMUONTrackHit*)
1227 (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1228 x = trackParam1->GetNonBendingCoor();
1229 y = trackParam1->GetBendingCoor();
1230 z = trackParam1->GetZ();
1231 bendingSlope = trackParam1->GetBendingSlope();
1232 nonBendingSlope = trackParam1->GetNonBendingSlope();
1233 pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1234 pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1235 pX = pZ * nonBendingSlope;
1236 pY = pZ * bendingSlope;
1237 c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1238 if (fPrintLevel >= 1)
1239 printf(" track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1240 z, x, y, pX, pY, pZ, c);
1242 // informations about generated particles
1243 particles = gAlice->Particles();
1244 np = particles->GetEntriesFast();
1245 printf(" **** number of generated particles: %d \n", np);
1247 for (Int_t iPart = 0; iPart < np; iPart++) {
1248 p = (TParticle*) particles->UncheckedAt(iPart);
1249 printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1250 iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
1255 void AliMUONEventReconstructor::Streamer(TBuffer &R__b)