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