X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONReconstructor.cxx;h=117e24f58068ebf84d1bae36af1a92c0b7884509;hb=81bef2cea632e2a882d2a2d9b2a780a8737f7304;hp=9d47fa0b183d838180bd6757979b2a139720fddd;hpb=1197ff51a23f6c4c675e522e0da5003a117b55ac;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONReconstructor.cxx b/MUON/AliMUONReconstructor.cxx index 9d47fa0b183..117e24f5806 100644 --- a/MUON/AliMUONReconstructor.cxx +++ b/MUON/AliMUONReconstructor.cxx @@ -14,112 +14,296 @@ **************************************************************************/ /* $Id$ */ -/////////////////////////////////////////////////////////////////////////////// -// // -// class for MUON reconstruction // -// // -/////////////////////////////////////////////////////////////////////////////// -#include -#include +//----------------------------- +// Class AliMUONReconstructor +//----------------------------- +// Class for the +// MUON track reconstruction -#include "AliRunLoader.h" -#include "AliHeader.h" -#include "AliGenEventHeader.h" -#include "AliESD.h" #include "AliMUONReconstructor.h" -#include "AliMUONData.h" -#include "AliMUONEventReconstructor.h" +#include "AliESD.h" +#include "AliESDMuonTrack.h" +#include "AliLog.h" +#include "AliMUONConstants.h" +#include "AliMUONCalibrationData.h" +#include "AliMUONClusterFinderAZ.h" #include "AliMUONClusterReconstructor.h" -#include "AliMUONClusterFinderVS.h" +#include "AliMUONData.h" +#include "AliMUONDigitCalibrator.h" +#include "AliMUONEventRecoCombi.h" +#include "AliMUONDigitMaker.h" #include "AliMUONTrack.h" #include "AliMUONTrackParam.h" +#include "AliMUONTrackExtrap.h" +#include "AliMUONVTrackReconstructor.h" +#include "AliMUONTrackReconstructor.h" +#include "AliMUONTrackReconstructorK.h" #include "AliMUONTriggerTrack.h" -#include "AliESDMuonTrack.h" -#include "AliMUONRawData.h" +#include "AliMUONTriggerCircuit.h" +#include "AliMUONTriggerCrateStore.h" + +#include "AliMpSegmentation.h" +#include "AliMUONPreClusterFinder.h" +#include "AliMUONClusterFinderCOG.h" +#include "AliMUONClusterFinderSimpleFit.h" +#include "AliMUONClusterFinderMLEM.h" + #include "AliRawReader.h" +#include "AliRun.h" +#include "AliRunLoader.h" +#include "TTask.h" +#include "TStopwatch.h" +/// \cond CLASSIMP ClassImp(AliMUONReconstructor) +/// \endcond + //_____________________________________________________________________________ AliMUONReconstructor::AliMUONReconstructor() - : AliReconstructor() + : AliReconstructor(), + fRunLoader(0x0), + fDigitMaker(new AliMUONDigitMaker()), + fCalibrationData(0x0), + fCrateManager(new AliMUONTriggerCrateStore()), + fTriggerCircuit(new TClonesArray("AliMUONTriggerCircuit", 234)), + fTransformer(new AliMUONGeometryTransformer(kTRUE)) + { +/// Default constructor + + AliDebug(1,""); + // Crate manager + fCrateManager->ReadFromFile(); + + // set to digit maker + fDigitMaker->SetCrateManager(fCrateManager); + + // transformater + fTransformer->ReadGeometryData("volpath.dat", "geometry.root"); + + // trigger circuit + for (Int_t i = 0; i < AliMUONConstants::NTriggerCircuit(); i++) { + AliMUONTriggerCircuit* c = new AliMUONTriggerCircuit(); + c->SetTransformer(fTransformer); + c->Init(i,*fCrateManager); + TClonesArray& circuit = *fTriggerCircuit; + new(circuit[circuit.GetEntriesFast()])AliMUONTriggerCircuit(*c); + delete c; + } + + } + //_____________________________________________________________________________ AliMUONReconstructor::~AliMUONReconstructor() { +/// Destructor + + AliDebug(1,""); + delete fCalibrationData; + delete fDigitMaker; + delete fCrateManager; + delete fTriggerCircuit; + delete fTransformer; +} + +//_____________________________________________________________________________ +TTask* +AliMUONReconstructor::GetCalibrationTask(AliMUONData* data) const +{ +/// Create the calibration task(s). + + const AliRun* run = fRunLoader->GetAliRun(); + + AliInfo("Calibration will occur."); + Int_t runNumber = run->GetRunNumber(); + fCalibrationData = new AliMUONCalibrationData(runNumber); + if ( !fCalibrationData->IsValid() ) + { + AliError("Could not retrieve calibrations !"); + delete fCalibrationData; + fCalibrationData = 0x0; + return 0x0; + } + TTask* calibration = new TTask("MUONCalibrator","MUON Digit calibrator"); + calibration->Add(new AliMUONDigitCalibrator(data,fCalibrationData)); + //FIXME: calibration->Add(something about dead channels should go here). + return calibration; + +} + +//_____________________________________________________________________________ +void +AliMUONReconstructor::Init(AliRunLoader* runLoader) +{ +/// Initialize + + fRunLoader = runLoader; } + +//_____________________________________________________________________________ +AliMUONClusterReconstructor* +AliMUONReconstructor::CreateClusterReconstructor(AliMUONData* data) const +{ +/// Create cluster reconstructor + + AliMUONVClusterFinder* clusterFinder(0x0); + + TString opt(GetOption()); + opt.ToUpper(); + + if ( strstr(opt,"PRECLUSTER") ) + { + clusterFinder = new AliMUONPreClusterFinder; + } + else if ( strstr(opt,"COG") ) + { + clusterFinder = new AliMUONClusterFinderCOG; + } + else if ( strstr(opt,"SIMPLEFIT") ) + { + clusterFinder = new AliMUONClusterFinderSimpleFit; + } + else if ( strstr(opt,"MLEM:DRAW") ) + { + clusterFinder = new AliMUONClusterFinderMLEM(kTRUE); + } + else if ( strstr(opt,"MLEM") ) + { + clusterFinder = new AliMUONClusterFinderMLEM(kFALSE); + } + + if ( clusterFinder) + { + AliInfo(Form("Will use %s for clusterizing",clusterFinder->ClassName())); + } + + AliMUONClusterReconstructor* clusterReco = + new AliMUONClusterReconstructor(data,clusterFinder,fTransformer); + return clusterReco; +} + //_____________________________________________________________________________ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const { -// AliLoader +/// Reconstruct +/// \todo add more + AliLoader* loader = runLoader->GetLoader("MUONLoader"); - Int_t nEvents = runLoader->GetNumberOfEvents(); + Int_t nEvents = runLoader->GetNumberOfEvents(); + Int_t evtNumber = runLoader->GetEventNumber(); -// used local container for each method -// passing fLoader as argument, could be avoided ??? - AliMUONEventReconstructor* recoEvent = new AliMUONEventReconstructor(loader); - AliMUONData* dataEvent = recoEvent->GetMUONData(); + AliMUONData* data = new AliMUONData(loader,"MUON","MUON"); - AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader); - AliMUONData* dataCluster = recoCluster->GetMUONData(); +// passing loader as argument. + AliMUONVTrackReconstructor* recoEvent; + if (strstr(GetOption(),"Original")) recoEvent = new AliMUONTrackReconstructor(data); + else if (strstr(GetOption(),"Combi")) recoEvent = new AliMUONTrackReconstructorK(data,"Combi"); + else recoEvent = new AliMUONTrackReconstructorK(data,"Kalman"); + + recoEvent->SetTriggerCircuit(fTriggerCircuit); + + AliMUONClusterReconstructor* recoCluster = CreateClusterReconstructor(data); + AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel(); + + if (!strstr(GetOption(),"VS")) { + recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ(); + recoCluster->SetRecoModel(recModel); + } recModel->SetGhostChi2Cut(10); + recModel->SetEventNumber(evtNumber); loader->LoadDigits("READ"); loader->LoadRecPoints("RECREATE"); loader->LoadTracks("RECREATE"); + TTask* calibration = GetCalibrationTask(data); + + Int_t chBeg = (strstr(GetOption(),"Combi") ? 6 : 0); + // Loop over events for(Int_t ievent = 0; ievent < nEvents; ievent++) { - printf("Event %d\n",ievent); + + AliDebug(1,Form("Event %d",ievent)); + runLoader->GetEvent(ievent); //----------------------- digit2cluster & Trigger2Trigger ------------------- if (!loader->TreeR()) loader->MakeRecPointsContainer(); // tracking branch - dataCluster->MakeBranch("RC"); - dataCluster->SetTreeAddress("D,RC"); - recoCluster->Digits2Clusters(); - dataCluster->Fill("RC"); + if (!strstr(GetOption(),"Combi")) { + data->MakeBranch("RC"); + data->SetTreeAddress("D,RC"); + } else { + data->SetTreeAddress("D"); + data->SetTreeAddress("RCC"); + } + // Important for avoiding a memory leak when reading digits ( to be investigated more in detail) + // In any case the reading of GLT is needed for the Trigger2Tigger method below + data->SetTreeAddress("GLT"); + + data->GetDigits(); + + if ( calibration ) + { + calibration->ExecuteTask(); + } + + recoCluster->Digits2Clusters(chBeg); + + if (strstr(GetOption(),"Combi")) { + // Combined cluster / track finder + AliMUONEventRecoCombi::Instance()->FillEvent(data, (AliMUONClusterFinderAZ*)recModel); + ((AliMUONClusterFinderAZ*) recModel)->SetReco(2); + } + else data->Fill("RC"); // trigger branch - dataCluster->MakeBranch("TC"); - dataCluster->SetTreeAddress("TC"); + data->MakeBranch("TC"); + data->SetTreeAddress("TC"); recoCluster->Trigger2Trigger(); - dataCluster->Fill("TC"); + data->Fill("TC"); - loader->WriteRecPoints("OVERWRITE"); + //AZ loader->WriteRecPoints("OVERWRITE"); //---------------------------- Track & TriggerTrack --------------------- if (!loader->TreeT()) loader->MakeTracksContainer(); // trigger branch - dataEvent->MakeBranch("RL"); //trigger track - dataEvent->SetTreeAddress("RL"); + data->MakeBranch("RL"); //trigger track + data->SetTreeAddress("RL"); recoEvent->EventReconstructTrigger(); - dataEvent->Fill("RL"); + data->Fill("RL"); // tracking branch - dataEvent->MakeBranch("RT"); //track - dataEvent->SetTreeAddress("RT"); + data->MakeBranch("RT"); //track + data->SetTreeAddress("RT"); recoEvent->EventReconstruct(); - dataEvent->Fill("RT"); + data->Fill("RT"); - loader->WriteTracks("OVERWRITE"); + loader->WriteTracks("OVERWRITE"); + if (strstr(GetOption(),"Combi")) { + // Combined cluster / track + ((AliMUONClusterFinderAZ*) recModel)->SetReco(1); + data->MakeBranch("RC"); + data->SetTreeAddress("RC"); + AliMUONEventRecoCombi::Instance()->FillRecP(data, (AliMUONTrackReconstructorK*)recoEvent); + data->Fill("RC"); + } + loader->WriteRecPoints("OVERWRITE"); + //--------------------------- Resetting branches ----------------------- - dataCluster->ResetDigits(); - dataCluster->ResetRawClusters(); - dataCluster->ResetTrigger(); - - dataEvent->ResetRawClusters(); - dataEvent->ResetTrigger(); - dataEvent->ResetRecTracks(); - dataEvent->ResetRecTriggerTracks(); - + data->ResetDigits(); + data->ResetRawClusters(); + data->ResetTrigger(); + data->ResetRecTracks(); + data->ResetRecTriggerTracks(); + } loader->UnloadDigits(); loader->UnloadRecPoints(); @@ -127,92 +311,177 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const delete recoCluster; delete recoEvent; + delete data; + delete calibration; } //_____________________________________________________________________________ -void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const +void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, + AliRawReader* rawReader) const { -// AliLoader +/// Recontruct +/// \todo add more + + // AliLoader AliLoader* loader = runLoader->GetLoader("MUONLoader"); + Int_t evtNumber = runLoader->GetEventNumber(); + + AliMUONData data(loader,"MUON","MUON"); -// used local container for each method -// passing fLoader as argument, could be avoided ??? - AliMUONEventReconstructor* recoEvent = new AliMUONEventReconstructor(loader); - AliMUONData* dataEvent = recoEvent->GetMUONData(); + // passing loader as argument. + fDigitMaker->SetMUONData(&data); - AliMUONRawData* rawData = new AliMUONRawData(loader); - AliMUONData* dataCluster = rawData->GetMUONData(); + AliMUONClusterReconstructor* recoCluster = CreateClusterReconstructor(&data); + + AliMUONVTrackReconstructor *recoEvent; + if (strstr(GetOption(),"Original")) recoEvent = new AliMUONTrackReconstructor(&data); + else if (strstr(GetOption(),"Combi")) recoEvent = new AliMUONTrackReconstructorK(&data,"Combi"); + else recoEvent = new AliMUONTrackReconstructorK(&data,"Kalman"); + + recoEvent->SetTriggerCircuit(fTriggerCircuit); - AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader, dataCluster); AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel(); + + if (!strstr(GetOption(),"VS")) + { + recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ(); + recoCluster->SetRecoModel(recModel); + } recModel->SetGhostChi2Cut(10); + recModel->SetEventNumber(evtNumber); + TTask* calibration = GetCalibrationTask(&data); + loader->LoadRecPoints("RECREATE"); loader->LoadTracks("RECREATE"); - + + // Digits are not stored on disk and created on flight from rawdata. + // In order to write digits on disk the following line should be uncommented + // loader->LoadDigits("RECREATE"); + // Loop over events Int_t iEvent = 0; - - while (rawReader->NextEvent()) { - printf("Event %d\n",iEvent); + + TStopwatch totalTimer; + TStopwatch rawTimer; + TStopwatch calibTimer; + TStopwatch clusterTimer; + TStopwatch trackingTimer; + + rawTimer.Start(kTRUE); rawTimer.Stop(); + calibTimer.Start(kTRUE); calibTimer.Stop(); + clusterTimer.Start(kTRUE); clusterTimer.Stop(); + trackingTimer.Start(kTRUE); trackingTimer.Stop(); + + totalTimer.Start(kTRUE); + + while (rawReader->NextEvent()) + { + AliDebug(1,Form("Event %d",iEvent)); + runLoader->GetEvent(iEvent++); - rawData->Raw2Digits(rawReader); + //----------------------- raw2digits & raw2trigger------------------- +// if (!loader->TreeD()) +// { +// AliDebug(1,Form("Making Digit Container for event %d",iEvent)); +// loader->MakeDigitsContainer(); +// } + // Digits are not stored on disk and created on flight from rawdata. + // In order to write digits on disk the following lines should be uncommented + // data.MakeBranch("D,GLT"); + // data.SetTreeAddress("D,GLT"); + data.SetDataContainer("D, GLT"); + rawTimer.Start(kFALSE); + fDigitMaker->Raw2Digits(rawReader); + rawTimer.Stop(); + + if ( calibration ) + { + calibTimer.Start(kFALSE); + calibration->ExecuteTask(); + calibTimer.Stop(); + } + // Digits are not stored on disk and created on flight from rawdata. + // In order to write digits on disk the following lines should be uncommented + // data.Fill("D,GLT"); + // loader->WriteDigits("OVERWRITE"); //----------------------- digit2cluster & Trigger2Trigger ------------------- + clusterTimer.Start(kFALSE); + if (!loader->TreeR()) loader->MakeRecPointsContainer(); // tracking branch - dataCluster->MakeBranch("RC"); - dataCluster->SetTreeAddress("D,RC"); + data.MakeBranch("RC"); + data.SetTreeAddress("RC"); recoCluster->Digits2Clusters(); - dataCluster->Fill("RC"); + data.Fill("RC"); // trigger branch - dataCluster->MakeBranch("TC"); - dataCluster->SetTreeAddress("TC"); - recoCluster->Trigger2Trigger(); - dataCluster->Fill("TC"); - + data.MakeBranch("TC"); + data.SetTreeAddress("TC"); + data.Fill("TC"); + loader->WriteRecPoints("OVERWRITE"); + clusterTimer.Stop(); + //---------------------------- Track & TriggerTrack --------------------- + trackingTimer.Start(kFALSE); if (!loader->TreeT()) loader->MakeTracksContainer(); // trigger branch - dataEvent->MakeBranch("RL"); //trigger track - dataEvent->SetTreeAddress("RL"); + data.MakeBranch("RL"); //trigger track + data.SetTreeAddress("RL"); recoEvent->EventReconstructTrigger(); - dataEvent->Fill("RL"); + data.Fill("RL"); // tracking branch - dataEvent->MakeBranch("RT"); //track - dataEvent->SetTreeAddress("RT"); + data.MakeBranch("RT"); //track + data.SetTreeAddress("RT"); recoEvent->EventReconstruct(); - dataEvent->Fill("RT"); + data.Fill("RT"); loader->WriteTracks("OVERWRITE"); - + trackingTimer.Stop(); + //--------------------------- Resetting branches ----------------------- - dataCluster->ResetDigits(); - dataCluster->ResetRawClusters(); - dataCluster->ResetTrigger(); - - dataEvent->ResetRawClusters(); - dataEvent->ResetTrigger(); - dataEvent->ResetRecTracks(); - dataEvent->ResetRecTriggerTracks(); + data.ResetDigits(); + data.ResetRawClusters(); + data.ResetTrigger(); + data.ResetRecTracks(); + data.ResetRecTriggerTracks(); } + + totalTimer.Stop(); + loader->UnloadRecPoints(); loader->UnloadTracks(); + loader->UnloadDigits(); + + delete recoEvent; delete recoCluster; - delete recoEvent; + + AliInfo(Form("Execution time for converting RAW data to digits in MUON : R:%.2fs C:%.2fs", + rawTimer.RealTime(),rawTimer.CpuTime())); + AliInfo(Form("Execution time for calibrating MUON : R:%.2fs C:%.2fs", + calibTimer.RealTime(),calibTimer.CpuTime())); + AliInfo(Form("Execution time for clusterizing MUON : R:%.2fs C:%.2fs", + clusterTimer.RealTime(),clusterTimer.CpuTime())); + AliInfo(Form("Execution time for tracking MUON : R:%.2fs C:%.2fs", + trackingTimer.RealTime(),trackingTimer.CpuTime())); + AliInfo(Form("Total Execution time for Reconstruct(from raw) MUON : R:%.2fs C:%.2fs", + totalTimer.RealTime(),totalTimer.CpuTime())); } //_____________________________________________________________________________ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const { +/// Fill ESD +/// \todo add more + TClonesArray* recTracksArray = 0; TClonesArray* recTrigTracksArray = 0; @@ -221,10 +490,9 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON"); // declaration - Int_t iEvent, nPart; - Int_t nTrackHits, nPrimary; + Int_t iEvent;// nPart; + Int_t nTrackHits;// nPrimary; Double_t fitFmin; - TArrayF vertex(3); Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum; Double_t xRec, yRec, zRec, chi2MatchTrigger; @@ -234,34 +502,16 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const AliMUONTrack* recTrack = 0; AliMUONTrackParam* trackParam = 0; AliMUONTriggerTrack* recTriggerTrack = 0; - TParticle* particle = new TParticle(); - AliGenEventHeader* header = 0; + iEvent = runLoader->GetEventNumber(); runLoader->GetEvent(iEvent); - // vertex calculation (maybe it exists already somewhere else) - vertex[0] = vertex[1] = vertex[2] = 0.; - nPrimary = 0; - if ( (header = runLoader->GetHeader()->GenEventHeader()) ) { - header->PrimaryVertex(vertex); - } else { - runLoader->LoadKinematics("READ"); - runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle); - nPart = (Int_t)runLoader->TreeK()->GetEntries(); - for(Int_t iPart = 0; iPart < nPart; iPart++) { - runLoader->TreeK()->GetEvent(iPart); - if (particle->GetFirstMother() == -1) { - vertex[0] += particle->Vx(); - vertex[1] += particle->Vy(); - vertex[2] += particle->Vz(); - nPrimary++; - } - if (nPrimary) { - vertex[0] /= (double)nPrimary; - vertex[1] /= (double)nPrimary; - vertex[2] /= (double)nPrimary; - } - } + // Get vertex + Double_t vertex[3] = {0}; + const AliESDVertex *esdVert = esd->GetVertex(); + if (esdVert->GetNContributors()) { + esdVert->GetXYZ(vertex); + printf("find vertex\n"); } // setting ESD MUON class AliESDMuonTrack* theESDTrack = new AliESDMuonTrack() ; @@ -295,11 +545,13 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks); trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First(); - trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]); + + if (esdVert->GetNContributors()) + AliMUONTrackExtrap::ExtrapToVertex(trackParam, vertex[0],vertex[1],vertex[2]); bendingSlope = trackParam->GetBendingSlope(); nonBendingSlope = trackParam->GetNonBendingSlope(); - inverseBendingMomentum = trackParam->GetInverseBendingMomentum(); + inverseBendingMomentum = trackParam->GetInverseBendingMomentum(); xRec = trackParam->GetNonBendingCoor(); yRec = trackParam->GetBendingCoor(); zRec = trackParam->GetZ(); @@ -326,19 +578,21 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const esd->AddMuonTrack(theESDTrack); } // end loop tracks - // add global trigger pattern - if (nRecTracks != 0) - esd->SetTrigger(trigPat); - // reset muondata muonData->ResetRecTracks(); muonData->ResetRecTriggerTracks(); //} // end loop on event loader->UnloadTracks(); - if (!header) - runLoader->UnloadKinematics(); + delete theESDTrack; delete muonData; - // delete particle; +}//_____________________________________________________________________________ +void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const +{ +/// Fill ESD +/// \todo add more + + // don't need rawReader ??? + FillESD(runLoader, esd); }