/* $Id$ */
-#include "AliRun.h"
-#include "AliRunDigitizer.h"
-#include "AliRunLoader.h"
-
#include "AliMUONDigitizer.h"
#include "AliMUONConstants.h"
-#include "AliMUONChamber.h"
+#include "AliMUONSegmentation.h"
#include "AliMUONHitMapA1.h"
#include "AliMUON.h"
#include "AliMUONLoader.h"
-#include "AliMUONChamber.h"
#include "AliMUONConstants.h"
-#include "AliMUONDigitizer.h"
#include "AliMUONTransientDigit.h"
-#include "AliMUONHitMapA1.h"
#include "AliMUONTriggerDecision.h"
+#include "AliLog.h"
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONGeometryModule.h"
+#include "AliMUONGeometryStore.h"
+#include "AliRun.h"
+#include "AliRunDigitizer.h"
+#include "AliRunLoader.h"
/////////////////////////////////////////////////////////////////////////////////////
//
fTDList(0),
fTDCounter(0),
fMask(0),
- fSignal(0),
- fDebug(0)
+ fSignal(0)
{
// Default constructor.
// Initializes all pointers to NULL.
fMUON = NULL;
fMUONData = NULL;
fTrigDec = NULL;
-};
+}
//___________________________________________
AliMUONDigitizer::AliMUONDigitizer(AliRunDigitizer* manager) :
fTDList(0),
fTDCounter(0),
fMask(0),
- fSignal(0),
- fDebug(0)
+ fSignal(0)
{
// Constructor which should be used rather than the default constructor.
// Initializes all pointers to NULL.
fMUON = NULL;
fMUONData = NULL;
fTrigDec = NULL;
-};
+}
//___________________________________________
AliMUONDigitizer::AliMUONDigitizer(const AliMUONDigitizer& rhs)
{
// Protected copy constructor
- Fatal("AliMUONDigitizer", "Not implemented.");
+ AliFatal("Not implemented.");
}
//___________________________________________
if (this == &rhs) return *this;
- Fatal("operator=", "Not implemented.");
+ AliFatal("Not implemented.");
return *this;
}
-
+
//------------------------------------------------------------------------
Bool_t AliMUONDigitizer::Init()
{
}
//------------------------------------------------------------------------
-void AliMUONDigitizer::Exec(Option_t* option)
+void AliMUONDigitizer::Exec(Option_t* /*option*/)
{
// The main work loop starts here.
// The digitization process is broken up into two steps:
// 2) Loop over the generated transient digits and write them to the output
// stream. Done in CreateDigits()
- if (GetDebug() > 0) Info("Exec", "Running digitiser.");
- ParseOptions(option);
+ AliDebug(1, "Running digitiser.");
if (fManager->GetNinputs() == 0)
{
- Warning("Exec", "No inputs set, nothing to do.");
+ AliWarning("No inputs set, nothing to do.");
return;
- };
+ }
if (!FetchLoaders(fManager->GetInputFolderName(0), fRunLoader, fGime) ) return;
if (! FetchGlobalPointers(fRunLoader) ) return;
InitArrays();
- if (GetDebug() > 1) Info("Exec", "Event Number is %d.", fManager->GetOutputEventNr());
+ AliDebug(2, Form("Event Number is %d.", fManager->GetOutputEventNr()));
// Loop over files to merge and to digitize
fSignal = kTRUE;
for (Int_t inputFile = 0; inputFile < fManager->GetNinputs(); inputFile++)
{
fMask = fManager->GetMask(inputFile);
- if (GetDebug() > 1)
- Info("Exec", "Digitising folder %d, with fMask = %d: %s", inputFile, fMask,
- (const char*)fManager->GetInputFolderName(inputFile));
+ AliDebug(2, Form("Digitising folder %d, with fMask = %d: %s", inputFile, fMask,
+ (const char*)fManager->GetInputFolderName(inputFile)));
if (inputFile != 0)
// If this is the first file then we already have the loaders loaded.
if (! InitInputData(fGime) ) continue;
GenerateTransientDigits();
CleanupInputData(fGime);
- };
+ }
Bool_t ok = FetchLoaders(fManager->GetOutputFolderName(), fRunLoader, fGime);
if (ok) ok = InitOutputData(fGime);
CleanupArrays();
CleanupTriggerArrays();
-};
+}
//--------------------------------------------------------------------------
void AliMUONDigitizer::AddOrUpdateTransientDigit(AliMUONTransientDigit* mTD)
}
else
AddTransientDigit(mTD);
-};
+}
//------------------------------------------------------------------------
void AliMUONDigitizer::UpdateTransientDigit(AliMUONTransientDigit* mTD)
// Update the transient digit that is already in the fTDList by adding the new
// transient digits charges and track lists to the existing one.
- if (GetDebug() > 3)
- Info("UpdateTransientDigit", "Updating transient digit 0x%X", (void*)mTD);
+ AliDebug(4,Form( "Updating transient digit 0x%X", (void*)mTD));
// Choosing the maping of the cathode plane of the chamber:
- Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
+ Int_t detElemId = mTD->DetElemId();
+
+ Int_t iNchCpl= fNDetElemId[detElemId] + (mTD->Cathode()-1) * AliMUONConstants::NDetElem();
+
AliMUONTransientDigit *pdigit =
static_cast<AliMUONTransientDigit*>( fHitMap[iNchCpl]->GetHit(mTD->PadX(),mTD->PadY()) );
Int_t ntracks = mTD->GetNTracks();
if (ntracks > kMAXTRACKS) // Truncate the number of tracks to kMAXTRACKS if we have to.
{
- if (GetDebug() > 0)
- {
- Warning("UpdateTransientDigit",
- "TransientDigit returned the number of tracks to be %d, which is bigger than kMAXTRACKS.",
- ntracks);
- Warning("UpdateTransientDigit", "Reseting the number of tracks to be %d.", kMAXTRACKS);
- }
+ AliDebug(1,Form(
+ "TransientDigit returned the number of tracks to be %d, which is bigger than kMAXTRACKS.",
+ ntracks));
+ AliDebug(1,Form( "Reseting the number of tracks to be %d.", kMAXTRACKS));
ntracks = kMAXTRACKS;
- };
+ }
for (Int_t i = 0; i < ntracks; i++)
{
pdigit->UpdateTrackList( mTD->GetTrack(i), mTD->GetCharge(i) );
- };
-};
+ }
+}
//------------------------------------------------------------------------
void AliMUONDigitizer::AddTransientDigit(AliMUONTransientDigit* mTD)
// Adds the transient digit to the fTDList and sets the appropriate entry
// in the fHitMap arrays.
- if (GetDebug() > 3)
- Info("AddTransientDigit", "Adding transient digit 0x%X", (void*)mTD);
+ AliDebug(4,Form( "Adding transient digit 0x%X", (void*)mTD));
// Choosing the maping of the cathode plane of the chamber:
- Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
+
+ Int_t detElemId = mTD->DetElemId();
+ Int_t iNchCpl= fNDetElemId[detElemId] + (mTD->Cathode()-1) * AliMUONConstants::NDetElem();
+
fTDList->AddAtAndExpand(mTD, fTDCounter);
- fHitMap[iNchCpl]->SetHit( mTD->PadX(), mTD->PadY(), fTDCounter);
- fTDCounter++;
-};
+ if (iNchCpl>-1 && iNchCpl<2*AliMUONConstants::NDetElem()) {
+ fHitMap[iNchCpl]->SetHit( mTD->PadX(), mTD->PadY(), fTDCounter);
+ fTDCounter++;
+ }
+}
//------------------------------------------------------------------------
Bool_t AliMUONDigitizer::ExistTransientDigit(AliMUONTransientDigit* mTD)
// as mTD. If yes then kTRUE is returned else kFASLE is returned.
// Choosing the maping of the cathode plane of the chamber:
- Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
- return( fHitMap[iNchCpl]->TestHit(mTD->PadX(), mTD->PadY()) );
-};
+ Int_t detElemId = mTD->DetElemId();
+
+ Int_t iNchCpl= fNDetElemId[detElemId] + (mTD->Cathode()-1) *AliMUONConstants::NDetElem() ;
+
+ // Int_t iNchCpl= mTD->Chamber() + (mTD->Cathode()-1) * AliMUONConstants::NCh();
+ if (iNchCpl>-1 && iNchCpl<2*AliMUONConstants::NDetElem())
+ return( fHitMap[iNchCpl]->TestHit(mTD->PadX(), mTD->PadY()) );
+ else return kFALSE;
+}
//-----------------------------------------------------------------------
void AliMUONDigitizer::CreateDigits()
// Loops over the fTDList for each cathode, gets the correct signal for the
// digit and adds the new digit to the output stream.
- if (GetDebug() > 1) Info("CreateDigits", "Creating digits...");
- for (Int_t icat = 0; icat < 2; icat++)
- {
- //
- // Filling Digit List
- Int_t nentries = fTDList->GetEntriesFast();
- for (Int_t nent = 0; nent < nentries; nent++)
- {
- AliMUONTransientDigit* td = (AliMUONTransientDigit*)fTDList->At(nent);
- if (td == NULL) continue;
+ fTDList->Sort(); // sort by idDE
+ AliDebug(2, "Creating digits...");
+ // for (Int_t icat = 0; icat < 2; icat++) {
+
+ Int_t digitindex[14] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+
+ //
+ // Filling Digit List
+ Int_t nentries = fTDList->GetEntriesFast();
+ for (Int_t nent = 0; nent < nentries; nent++) {
+
+ AliMUONTransientDigit* td = (AliMUONTransientDigit*)fTDList->At(nent);
+ if (td == NULL) continue;
- // Must be the same cathode, otherwise we will fill a mixture
- // of digits from both cathodes.
- if (icat != td->Cathode() - 1) continue;
+ // Must be the same cathode, otherwise we will fill a mixture
+ // of digits from both cathodes.
+ //if (icat != td->Cathode() - 1) continue;
- if (GetDebug() > 2)
- Info("CreateDigits", "Creating digit from transient digit 0x%X", (void*)td);
-
- Int_t q = GetSignalFrom(td);
- if (q > 0) AddDigit(td, q);
- };
- FillOutputData();
- };
-};
+ AliDebug(3,Form( "Creating digit from transient digit 0x%X", (void*)td));
+
+ Int_t q = GetSignalFrom(td);
+ if (q > 0) {
+ Int_t chamber = td->Chamber();
+ if (0 <= chamber && chamber <= 13 )
+ AddDigit(td, q, digitindex[chamber]++);
+ else
+ AliError(Form("Invalid chamber %d\n",chamber));
+ }
+ }
+ FillOutputData();
+ // }
+ fTDCounter = 0;
+}
//------------------------------------------------------------------------
-void AliMUONDigitizer::AddDigit(AliMUONTransientDigit* td, Int_t responseCharge)
+void AliMUONDigitizer::AddDigit(
+ AliMUONTransientDigit* td, Int_t responseCharge,
+ Int_t digitindex
+ )
{
// Prepares the digits, track and charge arrays in preparation for a call to
// AddDigit(Int_t, Int_t[kMAXTRACKS], Int_t[kMAXTRACKS], Int_t[6])
Int_t tracks[kMAXTRACKS];
Int_t charges[kMAXTRACKS];
- Int_t digits[6];
+ Int_t digits[7];
digits[0] = td->PadX();
digits[1] = td->PadY();
digits[3] = responseCharge;
digits[4] = td->Physics();
digits[5] = td->Hit();
-
+ digits[6] = td->DetElemId();
+
Int_t nptracks = td->GetNTracks();
- if (nptracks > kMAXTRACKS)
- {
- if (GetDebug() > 0)
- {
- Warning("AddDigit",
- "TransientDigit returned the number of tracks to be %d, which is bigger than kMAXTRACKS.",
- nptracks);
- Warning("AddDigit", "Reseting the number of tracks to be %d.", kMAXTRACKS);
- }
- nptracks = kMAXTRACKS;
- };
+ if (nptracks > kMAXTRACKS) {
+
+ AliDebug(1, Form(
+ "TransientDigit returned the number of tracks to be %d, which is bigger than kMAXTRACKS.",
+ nptracks));
+ AliDebug(1, Form("Reseting the number of tracks to be %d.", kMAXTRACKS));
+ nptracks = kMAXTRACKS;
+ }
- for (Int_t i = 0; i < nptracks; i++)
- {
- tracks[i] = td->GetTrack(i);
- charges[i] = td->GetCharge(i);
- };
+ for (Int_t i = 0; i < nptracks; i++) {
+
+ tracks[i] = td->GetTrack(i);
+ charges[i] = td->GetCharge(i);
+ }
// Sort list of tracks according to charge
SortTracks(tracks,charges,nptracks);
- if (nptracks < kMAXTRACKS )
- {
- for (Int_t i = nptracks; i < kMAXTRACKS; i++)
- {
- tracks[i] = -1;
- charges[i] = 0;
- };
- };
+ if (nptracks < kMAXTRACKS ) {
- if (GetDebug() > 3) Info("AddDigit", "Adding digit with charge %d.", responseCharge);
+ for (Int_t i = nptracks; i < kMAXTRACKS; i++) {
+ tracks[i] = -1;
+ charges[i] = 0;
+ }
+ }
+
+ AliDebug(4,Form( "Adding digit with charge %d.", responseCharge));
OnWriteTransientDigit(td);
AddDigit(td->Chamber(), tracks, charges, digits);
- AddDigitTrigger(td->Chamber(), tracks, charges, digits);
-};
+ AddDigitTrigger(td->Chamber(), tracks, charges, digits, digitindex);
+}
//------------------------------------------------------------------------
void AliMUONDigitizer::OnCreateTransientDigit(AliMUONTransientDigit* /*digit*/, TObject* /*source_object*/)
//
// This is derived by Digitisers that want to trace which digits were made from
// which hits.
-};
+}
//------------------------------------------------------------------------
void AliMUONDigitizer::OnWriteTransientDigit(AliMUONTransientDigit* /*digit*/)
//
// This is derived by Digitisers that want to trace which digits were made from
// which hits.
-};
+}
//------------------------------------------------------------------------
Bool_t AliMUONDigitizer::FetchLoaders(const char* foldername, AliRunLoader*& runloader, AliMUONLoader*& muonloader)
// The muon loader is then loaded from the fetched run loader.
// kTRUE is returned if no error occurred otherwise kFALSE is returned.
- if (GetDebug() > 2)
- Info("FetchLoaders", "Fetching run loader and muon loader from folder: %s", foldername);
+ AliDebug(3, Form("Fetching run loader and muon loader from folder: %s", foldername));
runloader = AliRunLoader::GetRunLoader(foldername);
if (runloader == NULL)
{
- Error("FetchLoaders", "RunLoader not found in folder: %s", foldername);
+ AliError(Form("RunLoader not found in folder: %s", foldername));
return kFALSE;
}
muonloader = (AliMUONLoader*) runloader->GetLoader("MUONLoader");
if (muonloader == NULL)
{
- Error("FetchLoaders", "MUONLoader not found in folder: %s", foldername);
+ AliError(Form("MUONLoader not found in folder: %s", foldername));
return kFALSE;
}
return kTRUE;
-};
+}
//------------------------------------------------------------------------
Bool_t AliMUONDigitizer::FetchGlobalPointers(AliRunLoader* runloader)
// AliMUONData fetched from the MUON module.
// kTRUE is returned if no error occurred otherwise kFALSE is returned.
- if (GetDebug() > 2)
- Info("FetchGlobalPointers", "Fetching gAlice, MUON module and AliMUONData from runloader 0x%X.",
+ AliDebug(3, Form("Fetching gAlice, MUON module and AliMUONData from runloader 0x%X.",
(void*)runloader
- );
+ ));
if (runloader->GetAliRun() == NULL) runloader->LoadgAlice();
gAlice = runloader->GetAliRun();
if (gAlice == NULL)
{
- Error("FetchGlobalPointers", "Could not find the AliRun object in runloader 0x%X.", (void*)runloader);
+ AliError(Form("Could not find the AliRun object in runloader 0x%X.", (void*)runloader));
return kFALSE;
- };
+ }
fMUON = (AliMUON*) gAlice->GetDetector("MUON");
if (fMUON == NULL)
{
- Error("FetchGlobalPointers", "Could not find the MUON module in runloader 0x%X.", (void*)runloader);
+ AliError(Form("Could not find the MUON module in runloader 0x%X.", (void*)runloader));
return kFALSE;
- };
+ }
AliMUONLoader *muonloader = (AliMUONLoader*) runloader->GetLoader("MUONLoader");
if (muonloader == NULL)
{
- Error("FetchGlobalPointers", "MUONLoader not found ");
+ AliError( "MUONLoader not found ");
return kFALSE;
}
if (fMUONData == NULL) fMUONData = new AliMUONData(muonloader,"MUON","MUON");
if (fMUONData == NULL)
{
- Error("FetchGlobalPointers", "Could not find AliMUONData object in runloader 0x%X.", (void*)runloader);
+ AliError(Form("Could not find AliMUONData object in runloader 0x%X.", (void*)runloader));
return kFALSE;
- };
+ }
return kTRUE;
}
Bool_t AliMUONDigitizer::FetchTriggerPointer(AliMUONLoader* loader)
{
if (fMUONData == NULL) {
- Error("FetchTriggerPointer", "MUONData not found");
+ AliError("MUONData not found");
return kFALSE;
}
return kTRUE;
}
-//------------------------------------------------------------------------
-void AliMUONDigitizer::ParseOptions(Option_t* options)
-{
-// Called by the Exec method. ParseOptions should parse the option string given to the Exec method.
-//
-// The following options are defined:
-// "debug" - Sets the debug level to 99, which will show all debug messages.
-// "deb" - Same as "debug", implemented for backward comparability.
-//
-// If an invalid option is specified it is simply ignored.
-
- TString optionString = options;
- if (optionString.Data() == "debug" ||
- optionString.Data() == "deb" // maintained for compatability.
- )
- {
- Info("ParseOptions", "Called with option \"debug\".");
- SetDebug(99);
- };
-};
//------------------------------------------------------------------------
void AliMUONDigitizer::InitArrays()
//
// Note: the fTDList and fHitMap arrays must be NULL before calling this method.
- if (GetDebug() > 1) Info("InitArrays", "Initialising internal arrays.");
- if (GetDebug() > 3) Info("InitArrays", "Creating transient digits list.");
- fTDList = new TObjArray;
+ AliDebug(2, "Initialising internal arrays.");
+ AliDebug(4, "Creating transient digits list.");
+ fTDList = new TObjArray;
- // Array of pointer of the AliMUONHitMapA1:
- // two HitMaps per chamber, or one HitMap per cahtode plane
- fHitMap = new AliMUONHitMapA1* [2*AliMUONConstants::NCh()];
-
- // Loop over chambers for the definition AliMUONHitMap
- for (Int_t i = 0; i < AliMUONConstants::NCh(); i++)
- {
- if (GetDebug() > 3) Info("InitArrays", "Creating hit map for chamber %d, cathode 1.", i+1);
- AliMUONChamber* chamber = &(fMUON->Chamber(i));
- AliSegmentation* c1Segmentation = chamber->SegmentationModel(1); // Cathode plane 1
- fHitMap[i] = new AliMUONHitMapA1(c1Segmentation, fTDList);
- if (GetDebug() > 3) Info("InitArrays", "Creating hit map for chamber %d, cathode 2.", i+1);
- AliSegmentation* c2Segmentation = chamber->SegmentationModel(2); // Cathode plane 2
- fHitMap[i+AliMUONConstants::NCh()] = new AliMUONHitMapA1(c2Segmentation, fTDList);
- };
-};
+ // Array of pointer of the AliMUONHitMapA1:
+ // two HitMaps per chamber, or one HitMap per cahtode plane
+ fHitMap = new AliMUONHitMapA1* [2*AliMUONConstants::NDetElem()];
+ for (Int_t i=0; i<2*AliMUONConstants::NDetElem(); i++) fHitMap[i] = 0x0;
+
+ Int_t k = 0;
+ Int_t idDE;
+
+ for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
+
+
+ AliDebug(4,Form( "Creating hit map for chamber %d, cathode 1.", i+1));
+ AliMUONSegmentation* segmentation = fMUON->GetSegmentation();
+ AliMUONGeometrySegmentation* c1Segmentation
+ = segmentation->GetModuleSegmentation(i, 0); // Cathode plane 1
+ AliDebug(4,Form( "Creating hit map for chamber %d, cathode 2.", i+1));
+ AliMUONGeometrySegmentation* c2Segmentation
+ = segmentation->GetModuleSegmentation(i, 1); // Cathode plane 2
+
+ const AliMUONGeometryTransformer* kGeometryTransformer
+ = fMUON->GetGeometryTransformer();
+
+ AliMUONGeometryStore* detElements
+ = kGeometryTransformer->GetModuleTransformer(i)->GetDetElementStore();
+
+ // Loop over detection elements
+ for (Int_t j=0; j<detElements->GetNofEntries(); j++) {
+
+ idDE = detElements->GetEntry(j)->GetUniqueID();
+ fNDetElemId[idDE] = k;
+
+ Int_t npx1 = c1Segmentation->Npx(idDE)+1;
+ Int_t npy1 = c1Segmentation->Npy(idDE)+1;
+ fHitMap[k] = new AliMUONHitMapA1(npx1, npy1, fTDList);
+
+ Int_t npx2 = c2Segmentation->Npx(idDE)+1;
+ Int_t npy2 = c2Segmentation->Npy(idDE)+1;
+ fHitMap[k+AliMUONConstants::NDetElem()] = new AliMUONHitMapA1(npx2, npy2, fTDList);
+ k++;
+ }
+ }
+}
//------------------------------------------------------------------------
void AliMUONDigitizer::CleanupArrays()
{
// The arrays fTDList and fHitMap are deleted and the pointers set to NULL.
- if (GetDebug() > 1) Info("CleanupArrays", "Deleting internal arrays.");
- for(Int_t i = 0; i < 2*AliMUONConstants::NCh(); i++)
- {
- if (GetDebug() > 3) Info("CleanupArrays", "Deleting hit map for chamber %d, cathode %d.",
- i%AliMUONConstants::NCh()+1, i/AliMUONConstants::NCh()+1);
+ AliDebug(2, "Deleting internal arrays.");
+ for(Int_t i = 0; i < 2*AliMUONConstants::NDetElem(); i++) {
delete fHitMap[i];
- };
+ }
delete [] fHitMap;
fHitMap = NULL;
- if (GetDebug() > 3) Info("CleanupArrays", "Deleting transient digits list.");
+ AliDebug(4, "Deleting transient digits list.");
fTDList->Delete();
delete fTDList;
fTDList = NULL;
-};
+}
//------------------------------------------------------------------------
void AliMUONDigitizer::SortTracks(Int_t *tracks, Int_t *charges, Int_t ntr) const
// Only the 3 most significant tracks are actually sorted
//
- if (ntr <= 1) return;
+ if (ntr <= 1) return;
- //
- // Loop over signals, only 3 times
- //
+ //
+ // Loop over signals, only 3 times
+ //
- Int_t qmax;
- Int_t jmax;
- Int_t idx[3] = {-2,-2,-2};
- Int_t jch[3] = {-2,-2,-2};
- Int_t jtr[3] = {-2,-2,-2};
- Int_t i, j, imax;
+ Int_t qmax;
+ Int_t jmax;
+ Int_t idx[3] = {-2,-2,-2};
+ Int_t jch[3] = {-2,-2,-2};
+ Int_t jtr[3] = {-2,-2,-2};
+ Int_t i, j, imax;
- if (ntr < 3) imax = ntr;
- else imax=3;
+ if (ntr < 3) imax = ntr;
+ else imax=3;
- for(i = 0; i < imax; i++)
- {
- qmax=0;
- jmax=0;
-
- for(j = 0; j < ntr; j++)
- {
- if ( (i == 1 && j == idx[i-1]) ||
- (i == 2 && (j == idx[i-1] || j == idx[i-2]))
- )
- continue;
-
- if(charges[j] > qmax)
- {
- qmax = charges[j];
- jmax = j;
- }
- }
-
- if(qmax > 0)
- {
- idx[i] = jmax;
- jch[i] = charges[jmax];
- jtr[i] = tracks[jmax];
- }
-
- }
-
- for(i = 0; i < 3; i++)
- {
- if (jtr[i] == -2)
- {
- charges[i] = 0;
- tracks[i] = 0;
- }
- else
- {
- charges[i] = jch[i];
- tracks[i] = jtr[i];
- }
- }
-};
+ for(i = 0; i < imax; i++)
+ {
+ qmax=0;
+ jmax=0;
+
+ for(j = 0; j < ntr; j++)
+ {
+ if ( (i == 1 && j == idx[i-1]) ||
+ (i == 2 && (j == idx[i-1] || j == idx[i-2]))
+ )
+ continue;
+
+ if(charges[j] > qmax)
+ {
+ qmax = charges[j];
+ jmax = j;
+ }
+ }
+
+ if(qmax > 0)
+ {
+ idx[i] = jmax;
+ jch[i] = charges[jmax];
+ jtr[i] = tracks[jmax];
+ }
+
+ }
+
+ for(i = 0; i < 3; i++)
+ {
+ if (jtr[i] == -2)
+ {
+ charges[i] = 0;
+ tracks[i] = 0;
+ }
+ else
+ {
+ charges[i] = jch[i];
+ tracks[i] = jtr[i];
+ }
+ }
+}