X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONEventReconstructor.cxx;h=d048b13b8cebc5a5735d875482f4ebfad229fd75;hb=d12a7158071fadfaa7bffffba4e08a03b4c12342;hp=9bd2bffba1e7c5cbfddcd6c04d70519de3d4944d;hpb=d0f40f2310e6a9cffbb30fb9e8848715a0ee5e57;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONEventReconstructor.cxx b/MUON/AliMUONEventReconstructor.cxx index 9bd2bffba1e..d048b13b8ce 100644 --- a/MUON/AliMUONEventReconstructor.cxx +++ b/MUON/AliMUONEventReconstructor.cxx @@ -13,129 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.26 2001/05/03 08:11:31 hristov -stdlib.h included to define exit() - -Revision 1.25 2001/04/25 14:50:42 gosset -Corrections to violations of coding conventions - -Revision 1.24 2001/03/30 09:37:58 gosset -Initialisations of pointers... for GEANT background events in the constructor - -Revision 1.23 2001/01/26 21:44:45 morsch -Use access functions to AliMUONDigit, ... member data. - -Revision 1.22 2001/01/26 20:00:53 hristov -Major upgrade of AliRoot code -Revision 1.20 2001/01/08 11:01:02 gosset -Modifications used for addendum to Dimuon TDR (JP Cussonneau): -*. MaxBendingMomentum to make both a segment and a track (default 500) -*. MaxChi2 per degree of freedom to make a track (default 100) -*. MinBendingMomentum used also to make a track - and not only a segment (default 3) -*. wider roads for track search in stations 1 to 3 -*. extrapolation to actual Z instead of Z(chamber) in FollowTracks -*. in track fit: - - limits on parameters X and Y (+/-500) - - covariance matrices in double precision - - normalization of covariance matrices before inversion - - suppression of Minuit printouts -*. correction against memory leak (delete extrapHit) in FollowTracks -*. RMax to 10 degrees with Z(chamber) instead of fixed values; - RMin and Rmax cuts suppressed in NewHitForRecFromGEANT, - because useless with realistic geometry - -Revision 1.19 2000/11/20 11:24:10 gosset -New package for reconstructed tracks (A. Gheata): -* tree output in the file "tree_reco.root" -* display events and make histograms from this file - -Revision 1.18 2000/10/26 12:47:03 gosset -Real distance between chambers of each station taken into account -for the reconstruction parameters "fSegmentMaxDistBending[5]" - -Revision 1.17 2000/10/24 09:26:20 gosset -Comments updated - -Revision 1.16 2000/10/24 09:22:35 gosset -Method AddHitsForRecFromRawClusters: real Z of raw cluster and not Z of chamber - -Revision 1.15 2000/10/12 15:17:30 gosset -Sign of fSimpleBValue corrected: sign ox Bx and not Bz (thanks to Galina) - -Revision 1.14 2000/10/04 18:21:26 morsch -Include stdlib.h - -Revision 1.13 2000/10/02 21:28:09 fca -Removal of useless dependecies via forward declarations - -Revision 1.12 2000/10/02 16:58:29 egangler -Cleaning of the code : --> coding conventions --> void Streamers --> some useless includes removed or replaced by "class" statement - -Revision 1.11 2000/09/22 09:16:33 hristov -Casting needed on DEC - -Revision 1.10 2000/09/19 09:49:50 gosset -AliMUONEventReconstructor package -* track extrapolation independent from reco_muon.F, use of AliMagF... -* possibility to use new magnetic field (automatic from generated root file) - -Revision 1.9 2000/07/20 12:45:27 gosset -New "EventReconstructor..." structure, - hopefully more adapted to tree/streamer. -"AliMUONEventReconstructor::RemoveDoubleTracks" - to keep only one track among similar ones. - -Revision 1.8 2000/07/18 16:04:06 gosset -AliMUONEventReconstructor package: -* a few minor modifications and more comments -* a few corrections - * right sign for Z of raw clusters - * right loop over chambers inside station - * symmetrized covariance matrix for measurements (TrackChi2MCS) - * right sign of charge in extrapolation (ExtrapToZ) - * right zEndAbsorber for Branson correction below 3 degrees -* use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit -* no parameter for AliMUONTrack::Fit() but more fit parameters in Track object - -Revision 1.7 2000/07/03 12:28:06 gosset -Printout at the right place after extrapolation to vertex - -Revision 1.6 2000/06/30 12:01:06 gosset -Correction for hit search in the right chamber (JPC) - -Revision 1.5 2000/06/30 10:15:48 gosset -Changes to EventReconstructor...: -precision fit with multiple Coulomb scattering; -extrapolation to vertex with Branson correction in absorber (JPC) - -Revision 1.4 2000/06/27 14:11:36 gosset -Corrections against violations of coding conventions - -Revision 1.3 2000/06/16 07:27:08 gosset -To remove problem in running RuleChecker, like in MUON-dev - -Revision 1.1.2.5 2000/06/16 07:00:26 gosset -To remove problem in running RuleChecker - -Revision 1.1.2.4 2000/06/12 08:00:07 morsch -Dummy streamer to solve CINT compilation problem (to be investigated !) - -Revision 1.1.2.3 2000/06/09 20:59:57 morsch -Make includes consistent with new file structure. - -Revision 1.1.2.2 2000/06/09 12:58:05 gosset -Removed comment beginnings in Log sections of .cxx files -Suppressed most violations of coding rules - -Revision 1.1.2.1 2000/06/07 14:44:53 gosset -Addition of files for track reconstruction in C++ -*/ +/* $Id$ */ //////////////////////////////////// // @@ -154,59 +32,71 @@ Addition of files for track reconstruction in C++ // //////////////////////////////////// -#include // for cout -#include // for exit() +#include +#include +#include +#include +#include //AZ -#include - -#include "AliMUON.h" -#include "AliMUONChamber.h" #include "AliMUONEventReconstructor.h" +#include "AliMUON.h" +#include "AliMUONHit.h" #include "AliMUONHitForRec.h" +#include "AliMUONTriggerTrack.h" +#include "AliMUONTriggerCircuit.h" #include "AliMUONRawCluster.h" +#include "AliMUONLocalTrigger.h" +#include "AliMUONGlobalTrigger.h" #include "AliMUONRecoEvent.h" #include "AliMUONSegment.h" #include "AliMUONTrack.h" #include "AliMUONTrackHit.h" #include "AliMagF.h" #include "AliRun.h" // for gAlice +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliMUONTrackK.h" //AZ +#include "AliMC.h" +#include "AliLog.h" //************* Defaults parameters for reconstruction -static const Double_t kDefaultMinBendingMomentum = 3.0; -static const Double_t kDefaultMaxBendingMomentum = 500.0; -static const Double_t kDefaultMaxChi2 = 100.0; -static const Double_t kDefaultMaxSigma2Distance = 16.0; -static const Double_t kDefaultBendingResolution = 0.01; -static const Double_t kDefaultNonBendingResolution = 0.144; -static const Double_t kDefaultChamberThicknessInX0 = 0.03; +const Double_t AliMUONEventReconstructor::fgkDefaultMinBendingMomentum = 3.0; +const Double_t AliMUONEventReconstructor::fgkDefaultMaxBendingMomentum = 500.0; +const Double_t AliMUONEventReconstructor::fgkDefaultMaxChi2 = 100.0; +const Double_t AliMUONEventReconstructor::fgkDefaultMaxSigma2Distance = 16.0; +const Double_t AliMUONEventReconstructor::fgkDefaultBendingResolution = 0.01; +const Double_t AliMUONEventReconstructor::fgkDefaultNonBendingResolution = 0.144; +const Double_t AliMUONEventReconstructor::fgkDefaultChamberThicknessInX0 = 0.03; // Simple magnetic field: // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG // Length and Position from reco_muon.F, with opposite sign: // Length = ZMAGEND-ZCOIL // Position = (ZMAGEND+ZCOIL)/2 // to be ajusted differently from real magnetic field ???? -static const Double_t kDefaultSimpleBValue = 7.0; -static const Double_t kDefaultSimpleBLength = 428.0; -static const Double_t kDefaultSimpleBPosition = 1019.0; -static const Int_t kDefaultRecGeantHits = 0; -static const Double_t kDefaultEfficiency = 0.95; +const Double_t AliMUONEventReconstructor::fgkDefaultSimpleBValue = 7.0; +const Double_t AliMUONEventReconstructor::fgkDefaultSimpleBLength = 428.0; +const Double_t AliMUONEventReconstructor::fgkDefaultSimpleBPosition = 1019.0; +const Int_t AliMUONEventReconstructor::fgkDefaultRecGeantHits = 0; +const Double_t AliMUONEventReconstructor::fgkDefaultEfficiency = 0.95; -static const Int_t kDefaultPrintLevel = 0; +const Int_t AliMUONEventReconstructor::fgkDefaultPrintLevel = -1; ClassImp(AliMUONEventReconstructor) // Class implementation in ROOT context //__________________________________________________________________________ -AliMUONEventReconstructor::AliMUONEventReconstructor(void) +AliMUONEventReconstructor::AliMUONEventReconstructor(AliLoader* loader) + : TObject() { // Constructor for class AliMUONEventReconstructor SetReconstructionParametersToDefaults(); + fTrackMethod = 1; //AZ - tracking method (1-default, 2-Kalman) // Memory allocation for the TClonesArray of hits for reconstruction // Is 10000 the right size ???? fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000); fNHitsForRec = 0; // really needed or GetEntriesFast sufficient ???? // Memory allocation for the TClonesArray's of segments in stations // Is 2000 the right size ???? - for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) { + for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) { fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000); fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ???? } @@ -214,16 +104,20 @@ AliMUONEventReconstructor::AliMUONEventReconstructor(void) // Is 10 the right size ???? fRecTracksPtr = new TClonesArray("AliMUONTrack", 10); fNRecTracks = 0; // really needed or GetEntriesFast sufficient ???? +// trigger tracks + fRecTriggerTracksPtr = new TClonesArray("AliMUONTriggerTrack", 10); + fNRecTriggerTracks = 0; // really needed or GetEntriesFast sufficient ???? // Memory allocation for the TClonesArray of hits on reconstructed tracks // Is 100 the right size ???? fRecTrackHitsPtr = new TClonesArray("AliMUONTrack", 100); fNRecTrackHits = 0; // really needed or GetEntriesFast sufficient ???? - // Sign of fSimpleBValue according to sign of Bx value at (50,50,950). + // Sign of fSimpleBValue according to sign of Bx value at (50,50,-950). Float_t b[3], x[3]; - x[0] = 50.; x[1] = 50.; x[2] = 950.; + x[0] = 50.; x[1] = 50.; x[2] = -950.; gAlice->Field()->Field(x, b); fSimpleBValue = TMath::Sign(fSimpleBValue,(Double_t) b[0]); + fSimpleBPosition = TMath::Sign(fSimpleBPosition,(Double_t) x[2]); // See how to get fSimple(BValue, BLength, BPosition) // automatically calculated from the actual magnetic field ???? @@ -246,19 +140,39 @@ AliMUONEventReconstructor::AliMUONEventReconstructor(void) fRecoEvent = 0; fEventTree = 0; fTreeFile = 0; - + + // initialize loader's + fLoader = loader; + + // initialize container + fMUONData = new AliMUONData(fLoader,"MUON","MUON"); + + // Loading AliRun master + AliRunLoader* runloader = fLoader->GetRunLoader(); + if (runloader->GetAliRun() == 0x0) runloader->LoadgAlice(); + gAlice = runloader->GetAliRun(); + return; } - -AliMUONEventReconstructor::AliMUONEventReconstructor (const AliMUONEventReconstructor& Reconstructor) + //__________________________________________________________________________ +AliMUONEventReconstructor::AliMUONEventReconstructor (const AliMUONEventReconstructor& rhs) + : TObject(rhs) { - // Dummy copy constructor +// Protected copy constructor + + AliFatal("Not implemented."); } -AliMUONEventReconstructor & AliMUONEventReconstructor::operator=(const AliMUONEventReconstructor& Reconstructor) +AliMUONEventReconstructor & +AliMUONEventReconstructor::operator=(const AliMUONEventReconstructor& rhs) { - // Dummy assignment operator - return *this; +// Protected assignement operator + + if (this == &rhs) return *this; + + AliFatal("Not implemented."); + + return *this; } //__________________________________________________________________________ @@ -272,7 +186,7 @@ AliMUONEventReconstructor::~AliMUONEventReconstructor(void) // if (fEventTree) delete fEventTree; if (fRecoEvent) delete fRecoEvent; delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ???? - for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) + for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) delete fSegmentsPtr[st]; // Correct destruction of everything ???? return; } @@ -282,10 +196,10 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void) { // Set reconstruction parameters to default values // Would be much more convenient with a structure (or class) ???? - fMinBendingMomentum = kDefaultMinBendingMomentum; - fMaxBendingMomentum = kDefaultMaxBendingMomentum; - fMaxChi2 = kDefaultMaxChi2; - fMaxSigma2Distance = kDefaultMaxSigma2Distance; + fMinBendingMomentum = fgkDefaultMinBendingMomentum; + fMaxBendingMomentum = fgkDefaultMaxBendingMomentum; + fMaxChi2 = fgkDefaultMaxChi2; + fMaxSigma2Distance = fgkDefaultMaxSigma2Distance; AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // ******** Parameters for making HitsForRec @@ -293,7 +207,7 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void) // like in TRACKF_STAT: // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3; // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9 - for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) { + for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) { if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) * 2.0 * TMath::Pi() / 180.0; else fRMin[ch] = 30.0; @@ -308,31 +222,31 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void) // Maximum distance in non bending plane // 5 * 0.22 just to remember the way it was made in TRACKF_STAT // SIGCUT*DYMAX(IZ) - for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) + for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) fSegmentMaxDistNonBending[st] = 5. * 0.22; // Maximum distance in bending plane: // values from TRACKF_STAT, corresponding to (J psi 20cm), // scaled to the real distance between chambers in a station - fSegmentMaxDistBending[0] = 1.5 * - ((&(pMUON->Chamber(1)))->Z() - (&(pMUON->Chamber(0)))->Z()) / 20.0; - fSegmentMaxDistBending[1] = 1.5 * - ((&(pMUON->Chamber(3)))->Z() - (&(pMUON->Chamber(2)))->Z()) / 20.0; - fSegmentMaxDistBending[2] = 3.0 * - ((&(pMUON->Chamber(5)))->Z() - (&(pMUON->Chamber(4)))->Z()) / 20.0; - fSegmentMaxDistBending[3] = 6.0 * - ((&(pMUON->Chamber(7)))->Z() - (&(pMUON->Chamber(6)))->Z()) / 20.0; - fSegmentMaxDistBending[4] = 6.0 * - ((&(pMUON->Chamber(9)))->Z() - (&(pMUON->Chamber(8)))->Z()) / 20.0; + fSegmentMaxDistBending[0] = TMath::Abs( 1.5 * + ((&(pMUON->Chamber(1)))->Z() - (&(pMUON->Chamber(0)))->Z()) / 20.0); + fSegmentMaxDistBending[1] = TMath::Abs( 1.5 * + ((&(pMUON->Chamber(3)))->Z() - (&(pMUON->Chamber(2)))->Z()) / 20.0); + fSegmentMaxDistBending[2] = TMath::Abs( 3.0 * + ((&(pMUON->Chamber(5)))->Z() - (&(pMUON->Chamber(4)))->Z()) / 20.0); + fSegmentMaxDistBending[3] = TMath::Abs( 6.0 * + ((&(pMUON->Chamber(7)))->Z() - (&(pMUON->Chamber(6)))->Z()) / 20.0); + fSegmentMaxDistBending[4] = TMath::Abs( 6.0 * + ((&(pMUON->Chamber(9)))->Z() - (&(pMUON->Chamber(8)))->Z()) / 20.0); - fBendingResolution = kDefaultBendingResolution; - fNonBendingResolution = kDefaultNonBendingResolution; - fChamberThicknessInX0 = kDefaultChamberThicknessInX0; - fSimpleBValue = kDefaultSimpleBValue; - fSimpleBLength = kDefaultSimpleBLength; - fSimpleBPosition = kDefaultSimpleBPosition; - fRecGeantHits = kDefaultRecGeantHits; - fEfficiency = kDefaultEfficiency; - fPrintLevel = kDefaultPrintLevel; + fBendingResolution = fgkDefaultBendingResolution; + fNonBendingResolution = fgkDefaultNonBendingResolution; + fChamberThicknessInX0 = fgkDefaultChamberThicknessInX0; + fSimpleBValue = fgkDefaultSimpleBValue; + fSimpleBLength = fgkDefaultSimpleBLength; + fSimpleBPosition = fgkDefaultSimpleBPosition; + fRecGeantHits = fgkDefaultRecGeantHits; + fEfficiency = fgkDefaultEfficiency; + fPrintLevel = fgkDefaultPrintLevel; return; } @@ -431,9 +345,8 @@ void AliMUONEventReconstructor::NextBkgGeantEvent(void) sprintf(treeName, "TreeK%d", fBkgGeantEventNumber); fBkgGeantTK = (TTree*)gDirectory->Get(treeName); if (!fBkgGeantTK) { - cout << "ERROR: cannot find Kine Tree for background event: " << - fBkgGeantEventNumber << endl; - exit(0); + AliError(Form("cannot find Kine Tree for background event: %d",fBkgGeantEventNumber)); + exit(0); } } if (fBkgGeantTK) @@ -443,8 +356,7 @@ void AliMUONEventReconstructor::NextBkgGeantEvent(void) sprintf(treeName, "TreeH%d", fBkgGeantEventNumber); fBkgGeantTH = (TTree*)gDirectory->Get(treeName); if (!fBkgGeantTH) { - cout << "ERROR: cannot find Hits Tree for background event: " << - fBkgGeantEventNumber << endl; + AliError(Form("cannot find Hits Tree for background event: %d",fBkgGeantEventNumber)); exit(0); } if (fBkgGeantTH && fBkgGeantHits) { @@ -466,6 +378,24 @@ void AliMUONEventReconstructor::EventReconstruct(void) MakeEventToBeReconstructed(); MakeSegments(); MakeTracks(); + if (fMUONData->IsTriggerTrackBranchesInTree()) + ValidateTracksWithTrigger(); + + // Add tracks to MUON data container + for(Int_t i=0; iAt(i); + fMUONData->AddRecTrack(*track); + } + + return; +} + +//__________________________________________________________________________ +void AliMUONEventReconstructor::EventReconstructTrigger(void) +{ + // To reconstruct one event + if (fPrintLevel >= 1) cout << "enter EventReconstructTrigger" << endl; + MakeTriggerTracks(); return; } @@ -477,7 +407,7 @@ void AliMUONEventReconstructor::ResetHitsForRec(void) // and the index of the first HitForRec per chamber if (fHitsForRecPtr) fHitsForRecPtr->Clear(); fNHitsForRec = 0; - for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) + for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0; return; } @@ -487,7 +417,7 @@ void AliMUONEventReconstructor::ResetSegments(void) { // To reset the TClonesArray of segments and the number of Segments // for all stations - for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) { + for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) { if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear(); fNSegments[st] = 0; } @@ -505,6 +435,17 @@ void AliMUONEventReconstructor::ResetTracks(void) return; } + //__________________________________________________________________________ +void AliMUONEventReconstructor::ResetTriggerTracks(void) +{ + // To reset the TClonesArray of reconstructed trigger tracks + if (fRecTriggerTracksPtr) fRecTriggerTracksPtr->Delete(); + // Delete in order that the Track destructors are called, + // hence the space for the TClonesArray of pointers to TrackHit's is freed + fNRecTriggerTracks = 0; + return; +} + //__________________________________________________________________________ void AliMUONEventReconstructor::ResetTrackHits(void) { @@ -520,16 +461,47 @@ void AliMUONEventReconstructor::MakeEventToBeReconstructed(void) // To make the list of hits to be reconstructed, // either from the GEANT hits or from the raw clusters // according to the parameter set for the reconstructor +// TString evfoldname = AliConfig::GetDefaultEventFolderName();//to be interfaced properly + +// AliRunLoader* rl = AliRunLoader::GetRunLoader(evfoldname); +// if (rl == 0x0) +// { +// Error("MakeEventToBeReconstructed", +// "Can not find Run Loader in Event Folder named %s.", +// evfoldname.Data()); +// return; +// } +// AliLoader* gime = rl->GetLoader("MUONLoader"); +// if (gime == 0x0) +// { +// Error("MakeEventToBeReconstructed","Can not get MUON Loader from Run Loader."); +// return; +// } + if (fPrintLevel >= 1) cout << "enter MakeEventToBeReconstructed" << endl; ResetHitsForRec(); if (fRecGeantHits == 1) { // Reconstruction from GEANT hits // Back to the signal file - ((gAlice->TreeK())->GetCurrentFile())->cd(); - // Signal hits - // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? - // Security on MUON ???? - AddHitsForRecFromGEANT(gAlice->TreeH()); + TTree* treeH = fLoader->TreeH(); + if (treeH == 0x0) + { + Int_t retval = fLoader->LoadHits(); + if ( retval) + { + AliError("Error occured while loading hits."); + return; + } + treeH = fLoader->TreeH(); + if (treeH == 0x0) + { + AliError("Can not get TreeH"); + return; + } + } + fMUONData->SetTreeAddress("H"); + AddHitsForRecFromGEANT(treeH); + // Background hits AddHitsForRecFromBkgGEANT(fBkgGeantTH, fBkgGeantHits); // Sort HitsForRec in increasing order with respect to chamber number @@ -541,14 +513,15 @@ void AliMUONEventReconstructor::MakeEventToBeReconstructed(void) // Security on MUON ???? // TreeR assumed to be be "prepared" in calling function // by "MUON->GetTreeR(nev)" ???? - TTree *treeR = gAlice->TreeR(); + TTree *treeR = fLoader->TreeR(); + fMUONData->SetTreeAddress("RC"); AddHitsForRecFromRawClusters(treeR); // No sorting: it is done automatically in the previous function } if (fPrintLevel >= 10) { cout << "end of MakeEventToBeReconstructed" << endl; cout << "NHitsForRec: " << fNHitsForRec << endl; - for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) { + for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) { cout << "chamber(0...): " << ch << " NHitsForRec: " << fNHitsForRecPerChamber[ch] << " index(first HitForRec): " << fIndexOfFirstHitForRecPerChamber[ch] @@ -569,26 +542,46 @@ void AliMUONEventReconstructor::AddHitsForRecFromGEANT(TTree *TH) { // To add to the list of hits for reconstruction // the GEANT signal hits from a hit tree TH. + Int_t hitBits, chamBits; //AZ if (fPrintLevel >= 2) cout << "enter AddHitsForRecFromGEANT with TH: " << TH << endl; if (TH == NULL) return; - AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? + // AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? + //AliMUONData * muondata = pMUON->GetMUONData(); // Security on MUON ???? // See whether it could be the same for signal and background ???? // Loop over tracks in tree Int_t ntracks = (Int_t) TH->GetEntries(); if (fPrintLevel >= 2) cout << "ntracks: " << ntracks << endl; + fMuons = 0; //AZ for (Int_t track = 0; track < ntracks; track++) { - gAlice->ResetHits(); + fMUONData->ResetHits(); TH->GetEvent(track); // Loop over hits Int_t hit = 0; - for (AliMUONHit* mHit = (AliMUONHit*) pMUON->FirstHit(-1); - mHit; - mHit = (AliMUONHit*) pMUON->NextHit(), hit++) { - NewHitForRecFromGEANT(mHit,track, hit, 1); - } // end of hit loop + hitBits = 0; // AZ + chamBits = 0; // AZ + Int_t itrack = track; //AZ + + Int_t ihit, nhits=0; + nhits = (Int_t) fMUONData->Hits()->GetEntriesFast(); + AliMUONHit* mHit=0x0; + + for(ihit=0; ihit(fMUONData->Hits()->At(ihit)); + Int_t ipart = TMath::Abs ((Int_t) mHit->Particle()); //AZ + if (NewHitForRecFromGEANT(mHit,track, hit, 1) && ipart == 13 + //if (NewHitForRecFromGEANT(mHit,itrack-1, hit, 1) && ipart == 13 + && itrack <= 2 && !BIT(mHit->Chamber()-1) ) chamBits |= BIT(mHit->Chamber()-1); //AZ - set bit + } + + if (chamBits&3 && chamBits>>2&3 && chamBits>>4&3 && chamBits>>6&3 && + chamBits>>8&3 && ((chamBits>>6&3)==3 || (chamBits>>8&3)==3)) + fMuons += 1; //AZ + //if (chamBits&3 && chamBits>>2&3 && chamBits>>4&3 && chamBits>>6&3 && + // chamBits>>8&3 && ((chamBits>>6&3)==3 || (chamBits>>8&3)==3) && + // ((chamBits&3)==3 || (chamBits>>2&3)==3)) fMuons += 1; } // end of track loop return; } @@ -636,7 +629,7 @@ AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* H Double_t bendCoor, nonBendCoor, radius; Int_t chamber = Hit->Chamber() - 1; // chamber(0...) // only in tracking chambers (fChamber starts at 1) - if (chamber >= kMaxMuonTrackingChambers) return NULL; + if (chamber >= fgkMaxMuonTrackingChambers) return NULL; // only if hit is efficient (keep track for checking ????) if (gRandom->Rndm() > fEfficiency) return NULL; // only if radius between RMin and RMax @@ -679,7 +672,7 @@ void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber() // Update the information for HitsForRec per chamber too. Int_t ch, nhits, prevch; fHitsForRecPtr->Sort(); - for (ch = 0; ch < kMaxMuonTrackingChambers; ch++) { + for (ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) { fNHitsForRecPerChamber[ch] = 0; fIndexOfFirstHitForRecPerChamber[ch] = 0; } @@ -714,7 +707,7 @@ void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber() // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? // // Security on MUON ???? // // Loop over tracking chambers -// for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) { +// for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) { // // number of HitsForRec to 0 for the chamber // fNHitsForRecPerChamber[ch] = 0; // // index of first HitForRec for the chamber @@ -758,21 +751,49 @@ void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR) // on the radius between RMin and RMax. AliMUONHitForRec *hitForRec; AliMUONRawCluster *clus; - Int_t iclus, nclus; + Int_t iclus, nclus, nTRentries; TClonesArray *rawclusters; if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromRawClusters" << endl; - AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? - // Security on MUON ???? + +// TString evfoldname = AliConfig::GetDefaultEventFolderName();//to be interfaced properly +// AliRunLoader* rl = AliRunLoader::GetRunLoader(evfoldname); +// if (rl == 0x0) +// { +// Error("MakeEventToBeReconstructed", +// "Can not find Run Loader in Event Folder named %s.", +// evfoldname.Data()); +// return; +// } +// AliLoader* gime = rl->GetLoader("MUONLoader"); +// if (gime == 0x0) +// { +// Error("MakeEventToBeReconstructed","Can not get MUON Loader from Run Loader."); +// return; +// } +// // Loading AliRun master +// rl->LoadgAlice(); +// gAlice = rl->GetAliRun(); + + // Loading MUON subsystem + // AliMUON * pMUON = (AliMUON *) gAlice->GetDetector("MUON"); + + nTRentries = Int_t(TR->GetEntries()); + if (nTRentries != 1) { + AliError(Form("nTRentries = %d not equal to 1 ",nTRentries)); + exit(0); + } + fLoader->TreeR()->GetEvent(0); // only one entry + // Loop over tracking chambers - for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) { + for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) { // number of HitsForRec to 0 for the chamber fNHitsForRecPerChamber[ch] = 0; // index of first HitForRec for the chamber if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0; else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec; - rawclusters = pMUON->RawClustAddress(ch); - pMUON->ResetRawClusters(); - TR->GetEvent((Int_t) (TR->GetEntries()) - 1); // to be checked ???? + rawclusters =fMUONData->RawClusters(ch); +// pMUON->ResetRawClusters(); +// TR->GetEvent((Int_t) (TR->GetEntries()) - 1); // to be checked ???? nclus = (Int_t) (rawclusters->GetEntries()); // Loop over (cathode correlated) raw clusters for (iclus = 0; iclus < nclus; iclus++) { @@ -790,7 +811,7 @@ void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR) hitForRec->SetChamberNumber(ch); hitForRec->SetHitNumber(iclus); // Z coordinate of the raw cluster (cm) - hitForRec->SetZ(clus->fZ[0]); + hitForRec->SetZ(clus->GetZ(0)); if (fPrintLevel >= 10) { cout << "chamber (0...): " << ch << " raw cluster (0...): " << iclus << endl; @@ -799,6 +820,7 @@ void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR) hitForRec->Dump();} } // end of cluster loop } // end of chamber loop + SortHitsForRecWithIncreasingChamber(); //AZ return; } @@ -810,11 +832,13 @@ void AliMUONEventReconstructor::MakeSegments(void) if (fPrintLevel >= 1) cout << "enter MakeSegments" << endl; ResetSegments(); // Loop over stations - for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) + Int_t nb = (fTrackMethod == 2) ? 3 : 0; //AZ + //AZ for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) + for (Int_t st = nb; st < fgkMaxMuonTrackingStations; st++) //AZ MakeSegmentsPerStation(st); if (fPrintLevel >= 10) { cout << "end of MakeSegments" << endl; - for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) { + for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) { cout << "station(0...): " << st << " Segments: " << fNSegments[st] << endl; @@ -894,7 +918,7 @@ void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station) (hit1Ptr->GetZ() - hit2Ptr->GetZ()); // absolute value of impact parameter impactParam = - TMath::Abs(hit1Ptr->GetBendingCoor() - hit2Ptr->GetZ() * bendingSlope); + TMath::Abs(hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope); } // check for distances not too large, // and impact parameter not too big if stations downstream of the dipole. @@ -948,15 +972,132 @@ void AliMUONEventReconstructor::MakeTracks(void) // The order may be important for the following Reset's ResetTracks(); ResetTrackHits(); - // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5 - MakeTrackCandidates(); - // Follow tracks in stations(1..) 3, 2 and 1 - FollowTracks(); - // Remove double tracks - RemoveDoubleTracks(); + if (fTrackMethod == 2) { //AZ - Kalman filter + MakeTrackCandidatesK(); + // Follow tracks in stations(1..) 3, 2 and 1 + FollowTracksK(); + // Remove double tracks + RemoveDoubleTracksK(); + // Propagate tracks to the vertex thru absorber + GoToVertex(); + } else { + // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5 + MakeTrackCandidates(); + // Follow tracks in stations(1..) 3, 2 and 1 + FollowTracks(); + // Remove double tracks + RemoveDoubleTracks(); + UpdateTrackParamAtHit(); + } return; } + //__________________________________________________________________________ +void AliMUONEventReconstructor::ValidateTracksWithTrigger(void) +{ + // Try to match track from tracking system with trigger track + AliMUONTrack *track; + TClonesArray *recTriggerTracks; + + fMUONData->ResetRecTriggerTracks(); + fMUONData->SetTreeAddress("RL"); + fMUONData->GetRecTriggerTracks(); + recTriggerTracks = fMUONData->RecTriggerTracks(); + + track = (AliMUONTrack*) fRecTracksPtr->First(); + while (track) { + track->MatchTriggerTrack(recTriggerTracks); + track = (AliMUONTrack*) fRecTracksPtr->After(track); + } + +} + + //__________________________________________________________________________ +Bool_t AliMUONEventReconstructor::MakeTriggerTracks(void) +{ + // To make the trigger tracks from Local Trigger + if (fPrintLevel >= 1) cout << "enter MakeTriggerTracks" << endl; + // ResetTriggerTracks(); + + Int_t nTRentries; + Long_t gloTrigPat; + TClonesArray *localTrigger; + TClonesArray *globalTrigger; + AliMUONLocalTrigger *locTrg; + AliMUONGlobalTrigger *gloTrg; + AliMUONTriggerCircuit *circuit; + AliMUONTriggerTrack *recTriggerTrack = 0; + + TTree* treeR = fLoader->TreeR(); + + // Loading MUON subsystem + AliMUON * pMUON = (AliMUON *) gAlice->GetDetector("MUON"); + + nTRentries = Int_t(treeR->GetEntries()); + + treeR->GetEvent(0); // only one entry + + if (!(fMUONData->IsTriggerBranchesInTree())) { + AliWarning(Form("Trigger information is not avalaible, nTRentries = %d not equal to 1",nTRentries)); + return kFALSE; + } + + fMUONData->SetTreeAddress("TC"); + fMUONData->GetTrigger(); + + // global trigger for trigger pattern + gloTrigPat = 0; + globalTrigger = fMUONData->GlobalTrigger(); + gloTrg = (AliMUONGlobalTrigger*)globalTrigger->UncheckedAt(0); + if (gloTrg) { + if (gloTrg->SinglePlusLpt()) gloTrigPat|= 0x1; + if (gloTrg->SinglePlusHpt()) gloTrigPat|= 0x2; + if (gloTrg->SinglePlusApt()) gloTrigPat|= 0x4; + + if (gloTrg->SingleMinusLpt()) gloTrigPat|= 0x8; + if (gloTrg->SingleMinusHpt()) gloTrigPat|= 0x10; + if (gloTrg->SingleMinusApt()) gloTrigPat|= 0x20; + + if (gloTrg->SingleUndefLpt()) gloTrigPat|= 0x40; + if (gloTrg->SingleUndefHpt()) gloTrigPat|= 0x80; + if (gloTrg->SingleUndefApt()) gloTrigPat|= 0x100; + + if (gloTrg->PairUnlikeLpt()) gloTrigPat|= 0x200; + if (gloTrg->PairUnlikeHpt()) gloTrigPat|= 0x400; + if (gloTrg->PairUnlikeApt()) gloTrigPat|= 0x800; + + if (gloTrg->PairLikeLpt()) gloTrigPat|= 0x1000; + if (gloTrg->PairLikeHpt()) gloTrigPat|= 0x2000; + if (gloTrg->PairLikeApt()) gloTrigPat|= 0x4000; + } + + + + // local trigger for tracking + localTrigger = fMUONData->LocalTrigger(); + Int_t nlocals = (Int_t) (localTrigger->GetEntries()); + Float_t z11 = ( &(pMUON->Chamber(10)) )->Z(); + Float_t z21 = ( &(pMUON->Chamber(12)) )->Z(); + + for (Int_t i=0; iUncheckedAt(i); + circuit = &(pMUON->TriggerCircuit(locTrg->LoCircuit())); + Float_t y11 = circuit->GetY11Pos(locTrg->LoStripX()); + Int_t stripX21 = locTrg->LoStripX()+locTrg->LoDev()+1; + Float_t y21 = circuit->GetY21Pos(stripX21); + Float_t x11 = circuit->GetX11Pos(locTrg->LoStripY()); + Float_t thetax = TMath::ATan2( x11 , z11 ); + Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) ); + + recTriggerTrack = new AliMUONTriggerTrack(x11,y11,thetax,thetay,gloTrigPat,this); + // since static statement does not work, set gloTrigPat for each track + + // fNRecTriggerTracks++; + fMUONData->AddRecTriggerTrack(*recTriggerTrack); + } // end of loop on Local Trigger + return kTRUE; +} + //__________________________________________________________________________ Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment) { @@ -1419,6 +1560,26 @@ void AliMUONEventReconstructor::RemoveDoubleTracks(void) return; } + //__________________________________________________________________________ +void AliMUONEventReconstructor::UpdateTrackParamAtHit() +{ + // Set track parameters after track fitting. Fill fTrackParamAtHit of AliMUONTrack's + AliMUONTrack *track; + AliMUONTrackHit *trackHit; + AliMUONTrackParam *trackParam; + track = (AliMUONTrack*) fRecTracksPtr->First(); + while (track) { + trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->First(); + while (trackHit) { + trackParam = trackHit->GetTrackParam(); + track->AddTrackParamAtHit(trackParam); + trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->After(trackHit); + } // trackHit + track = (AliMUONTrack*) fRecTracksPtr->After(track); + } // track + return; +} + //__________________________________________________________________________ void AliMUONEventReconstructor::EventDump(void) { @@ -1427,7 +1588,6 @@ void AliMUONEventReconstructor::EventDump(void) AliMUONTrack *track; AliMUONTrackParam *trackParam, *trackParam1; - TParticle *p; Double_t bendingSlope, nonBendingSlope, pYZ; Double_t pX, pY, pZ, x, y, z, c; Int_t np, trackIndex, nTrackHits; @@ -1439,6 +1599,7 @@ void AliMUONEventReconstructor::EventDump(void) fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole // Loop over reconstructed tracks for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) { + if (fTrackMethod != 1) continue; //AZ - skip the rest for now if (fPrintLevel >= 1) cout << " track number: " << trackIndex << endl; // function for each track for modularity ???? @@ -1480,17 +1641,42 @@ void AliMUONEventReconstructor::EventDump(void) z, x, y, pX, pY, pZ, c); } // informations about generated particles - np = gAlice->GetNtrack(); + np = gAlice->GetMCApp()->GetNtrack(); printf(" **** number of generated particles: %d \n", np); - for (Int_t iPart = 0; iPart < np; iPart++) { - p = gAlice->Particle(iPart); - printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n", - iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode()); - } +// for (Int_t iPart = 0; iPart < np; iPart++) { +// p = gAlice->Particle(iPart); +// printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n", +// iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode()); +// } return; } + +//__________________________________________________________________________ +void AliMUONEventReconstructor::EventDumpTrigger(void) +{ + // Dump reconstructed trigger event + // and the particle parameters + + AliMUONTriggerTrack *triggertrack ; + Int_t nTriggerTracks = fMUONData->RecTriggerTracks()->GetEntriesFast(); + + if (fPrintLevel >= 1) { + cout << "****** enter EventDumpTrigger ******" << endl; + cout << " Number of Reconstructed tracks :" << nTriggerTracks << endl; + } + // Loop over reconstructed tracks + for (Int_t trackIndex = 0; trackIndex < nTriggerTracks; trackIndex++) { + triggertrack = (AliMUONTriggerTrack*)fMUONData->RecTriggerTracks()->At(trackIndex); + printf(" trigger track number %i x11=%f y11=%f thetax=%f thetay=%f \n", + trackIndex, + triggertrack->GetX11(),triggertrack->GetY11(), + triggertrack->GetThetax(),triggertrack->GetThetay()); + } +} + +//__________________________________________________________________________ void AliMUONEventReconstructor::FillEvent() { // Create a new AliMUONRecoEvent, fill its track list, then add it as a @@ -1506,7 +1692,8 @@ void AliMUONEventReconstructor::FillEvent() TDirectory *current = gDirectory; if (!fTreeFile) fTreeFile = new TFile("tree_reco.root", "RECREATE"); if (!fEventTree) fEventTree = new TTree("TreeRecoEvent", "MUON reconstructed events"); - if (fRecoEvent->MakeDumpTracks(fRecTracksPtr)) { + //AZif (fRecoEvent->MakeDumpTracks(fRecTracksPtr)) { + if (fRecoEvent->MakeDumpTracks(fMuons, fRecTracksPtr, this)) { //AZ if (fPrintLevel > 1) fRecoEvent->EventInfo(); TBranch *branch = fEventTree->GetBranch("Event"); if (!branch) branch = fEventTree->Branch("Event", "AliMUONRecoEvent", &fRecoEvent, 64000); @@ -1518,3 +1705,239 @@ void AliMUONEventReconstructor::FillEvent() // restore directory current->cd(); } + +//__________________________________________________________________________ +void AliMUONEventReconstructor::MakeTrackCandidatesK(void) +{ + // To make initial tracks for Kalman filter from the list of segments + Int_t istat, iseg; + AliMUONSegment *segment; + AliMUONTrackK *trackK; + + if (fPrintLevel >= 1) cout << "enter MakeTrackCandidatesK" << endl; + // Reset the TClonesArray of reconstructed tracks + if (fRecTracksPtr) fRecTracksPtr->Delete(); + // Delete in order that the Track destructors are called, + // hence the space for the TClonesArray of pointers to TrackHit's is freed + fNRecTracks = 0; + + AliMUONTrackK a(this, fHitsForRecPtr); // bad idea ??? + // Loop over stations(1...) 5 and 4 + for (istat=4; istat>=3; istat--) { + // Loop over segments in the station + for (iseg=0; isegGetModule("MUON"); + for (Int_t i1=0; i1GetTHTrack() > 1 || hit->GetGeantSignal() == 0) continue; + /* + cout << " Hit #" << hit->GetChamberNumber() << " "; + cout << hit->GetBendingCoor() << " "; + cout << hit->GetNonBendingCoor() << " "; + cout << hit->GetZ() << " "; + cout << hit->GetGeantSignal() << " "; + cout << hit->GetTHTrack() << endl; + */ + /* + printf(" Hit # %d %10.4f %10.4f %10.4f", + hit->GetChamberNumber(), hit->GetBendingCoor(), + hit->GetNonBendingCoor(), hit->GetZ()); + if (fRecGeantHits) { + // from GEANT hits + printf(" %3d %3d \n", hit->GetGeantSignal(), hit->GetTHTrack()); + } else { + // from raw clusters + rawclusters = pMUON->RawClustAddress(hit->GetChamberNumber()); + clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit-> + GetHitNumber()); + printf("%3d", clus->fTracks[1]-1); + if (clus->fTracks[2] != 0) printf("%3d \n", clus->fTracks[2]-1); + else printf("\n"); + } + */ + } + + icand = -1; + Int_t nSeeds = fNRecTracks; // starting number of seeds + // Loop over track candidates + while (icand < fNRecTracks-1) { + icand ++; + trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]); + + // Discard candidate which will produce the double track + if (icand > 0) { + ok = CheckCandidateK(icand,nSeeds); + if (!ok) { + //trackK->SetRecover(-1); // mark candidate to be removed + //continue; + } + } + + ok = kTRUE; + if (trackK->GetRecover() == 0) hit = (AliMUONHitForRec*) + trackK->GetHitOnTrack()->Last(); // last hit + else hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[1]; // 2'nd hit + ichamBeg = hit->GetChamberNumber(); + ichamEnd = 0; + // Check propagation direction + if (trackK->GetTrackDir() > 0) { + ichamEnd = 9; // forward propagation + ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2); + if (ok) { + ichamBeg = ichamEnd; + ichamEnd = 6; // backward propagation + // Change weight matrix and zero fChi2 for backpropagation + trackK->StartBack(); + ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2); + ichamBeg = ichamEnd; + ichamEnd = 0; + } + } else { + if (trackK->GetBPFlag()) { + // backpropagation + ichamEnd = 6; // backward propagation + // Change weight matrix and zero fChi2 for backpropagation + trackK->StartBack(); + ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2); + ichamBeg = ichamEnd; + ichamEnd = 0; + } + } + + if (ok) { + trackK->SetTrackDir(-1); + trackK->SetBPFlag(kFALSE); + ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2); + } + if (ok) trackK->SetTrackQuality(0); // compute "track quality" + else trackK->SetRecover(-1); // mark candidate to be removed + + // Majority 3 of 4 in first 2 stations + chamBits = 0; + for (Int_t i=0; iGetNTrackHits(); i++) { + hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[i]; + chamBits |= BIT(hit->GetChamberNumber()-1); + } + //if (!((chamBits&3)==3 || (chamBits>>2&3)==3)) trackK->SetRecover(-1); + //mark candidate to be removed + } // while + + for (Int_t i=0; iGetRecover() < 0) fRecTracksPtr->RemoveAt(i); + } + + // Compress TClonesArray + fRecTracksPtr->Compress(); + fNRecTracks = fRecTracksPtr->GetEntriesFast(); + return; +} + +//__________________________________________________________________________ +Bool_t AliMUONEventReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds) const +{ + // Discards track candidate if it will produce the double track (having + // the same seed segment hits as hits of a good track found before) + AliMUONTrackK *track1, *track2; + AliMUONHitForRec *hit1, *hit2, *hit; + + track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]); + hit1 = (AliMUONHitForRec*) (*track1->GetHitOnTrack())[0]; // 1'st hit + hit2 = (AliMUONHitForRec*) (*track1->GetHitOnTrack())[1]; // 2'nd hit + + for (Int_t i=0; iGetRecover() < 0) continue; + if (track2->GetRecover() < 0 && icand >= nSeeds) continue; + + if (track1->GetStartSegment() == track2->GetStartSegment()) { + return kFALSE; + } else { + Int_t nSame = 0; + for (Int_t j=0; jGetNTrackHits(); j++) { + hit = (AliMUONHitForRec*) (*track2->GetHitOnTrack())[j]; + if (hit == hit1 || hit == hit2) { + nSame++; + if (nSame == 2) return kFALSE; + } + } // for (Int_t j=0; + } + } // for (Int_t i=0; + return kTRUE; +} + +//__________________________________________________________________________ +void AliMUONEventReconstructor::RemoveDoubleTracksK(void) +{ + // Removes double tracks (sharing more than half of their hits). Keeps + // the track with higher quality + AliMUONTrackK *track1, *track2, *trackToKill; + + // Sort tracks according to their quality + fRecTracksPtr->Sort(); + + // Loop over first track of the pair + track1 = (AliMUONTrackK*) fRecTracksPtr->First(); + while (track1) { + // Loop over second track of the pair + track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1); + while (track2) { + // Check whether or not to keep track2 + if (!track2->KeepTrack(track1)) { + cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) << + " " << track2->GetTrackQuality() << endl; + trackToKill = track2; + track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2); + trackToKill->Kill(); + fRecTracksPtr->Compress(); + } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2); + } // track2 + track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1); + } // track1 + + fNRecTracks = fRecTracksPtr->GetEntriesFast(); + cout << " Number of Kalman tracks: " << fNRecTracks << endl; +} + +//__________________________________________________________________________ +void AliMUONEventReconstructor::GoToVertex(void) +{ + // Propagates track to the vertex thru absorber + // (using Branson correction for now) + + Double_t zVertex; + zVertex = 0; + for (Int_t i=0; iBranson(); + //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber + ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2 + ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(); // with absorber + } +}