X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONReconstructor.cxx;h=117e24f58068ebf84d1bae36af1a92c0b7884509;hb=81bef2cea632e2a882d2a2d9b2a780a8737f7304;hp=c32ab8bafbb43724d4324aefff151b15831edb0d;hpb=30178c30974cdd6a3b59f09e4d479925642e175b;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONReconstructor.cxx b/MUON/AliMUONReconstructor.cxx index c32ab8bafbb..117e24f5806 100644 --- a/MUON/AliMUONReconstructor.cxx +++ b/MUON/AliMUONReconstructor.cxx @@ -12,123 +12,298 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - /* $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 "AliMUONTriggerDecision.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 "AliMUONReconstructor.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; } + //_____________________________________________________________________________ -void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const +TTask* +AliMUONReconstructor::GetCalibrationTask(AliMUONData* data) const { -// AliLoader - AliLoader* loader = runLoader->GetLoader("MUONLoader"); - Int_t nEvents = runLoader->GetNumberOfEvents(); +/// 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; -// used local container for each method -// passing fLoader as argument, could be avoided ??? - AliMUONEventReconstructor* recoEvent = new AliMUONEventReconstructor(loader); - AliMUONData* dataEvent = recoEvent->GetMUONData(); +} - AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader); - AliMUONData* dataCluster = recoCluster->GetMUONData(); +//_____________________________________________________________________________ +void +AliMUONReconstructor::Init(AliRunLoader* runLoader) +{ +/// Initialize - AliMUONTriggerDecision* trigDec = new AliMUONTriggerDecision(loader); - AliMUONData* dataTrig = trigDec->GetMUONData(); + fRunLoader = runLoader; +} +//_____________________________________________________________________________ +AliMUONClusterReconstructor* +AliMUONReconstructor::CreateClusterReconstructor(AliMUONData* data) const +{ +/// Create cluster reconstructor - for (Int_t i = 0; i < 10; i++) { - AliMUONClusterFinderVS *recModel = new AliMUONClusterFinderVS(); - recModel->SetGhostChi2Cut(10); - recoCluster->SetReconstructionModel(i,recModel); + 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 +{ +/// Reconstruct +/// \todo add more + + AliLoader* loader = runLoader->GetLoader("MUONLoader"); + Int_t nEvents = runLoader->GetNumberOfEvents(); + Int_t evtNumber = runLoader->GetEventNumber(); + + AliMUONData* data = new AliMUONData(loader,"MUON","MUON"); + +// 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 & Digits2Trigger ------------------- + //----------------------- 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 - dataTrig->MakeBranch("GLT"); - dataTrig->SetTreeAddress("D,GLT"); - trigDec->Digits2Trigger(); - dataTrig->Fill("GLT"); + data->MakeBranch("TC"); + data->SetTreeAddress("TC"); + recoCluster->Trigger2Trigger(); + 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"); - //--------------------------- Resetting branches ----------------------- - dataCluster->ResetDigits(); - dataCluster->ResetRawClusters(); + 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"); - dataTrig->ResetDigits(); - dataTrig->ResetTrigger(); + //--------------------------- Resetting branches ----------------------- + data->ResetDigits(); + data->ResetRawClusters(); + data->ResetTrigger(); + data->ResetRecTracks(); + data->ResetRecTriggerTracks(); - dataEvent->ResetRawClusters(); - dataEvent->ResetTrigger(); - dataEvent->ResetRecTracks(); - dataEvent->ResetRecTriggerTracks(); - } loader->UnloadDigits(); loader->UnloadRecPoints(); @@ -136,11 +311,177 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const delete recoCluster; delete recoEvent; - delete trigDec; + delete data; + delete calibration; } + +//_____________________________________________________________________________ +void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, + AliRawReader* rawReader) const +{ +/// Recontruct +/// \todo add more + + // AliLoader + AliLoader* loader = runLoader->GetLoader("MUONLoader"); + Int_t evtNumber = runLoader->GetEventNumber(); + + AliMUONData data(loader,"MUON","MUON"); + + // passing loader as argument. + fDigitMaker->SetMUONData(&data); + + 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); + + 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; + + 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++); + + //----------------------- 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 + data.MakeBranch("RC"); + data.SetTreeAddress("RC"); + recoCluster->Digits2Clusters(); + data.Fill("RC"); + + // trigger branch + 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 + data.MakeBranch("RL"); //trigger track + data.SetTreeAddress("RL"); + recoEvent->EventReconstructTrigger(); + data.Fill("RL"); + + // tracking branch + data.MakeBranch("RT"); //track + data.SetTreeAddress("RT"); + recoEvent->EventReconstruct(); + data.Fill("RT"); + + loader->WriteTracks("OVERWRITE"); + trackingTimer.Stop(); + + //--------------------------- Resetting branches ----------------------- + data.ResetDigits(); + data.ResetRawClusters(); + data.ResetTrigger(); + data.ResetRecTracks(); + data.ResetRecTriggerTracks(); + + } + + totalTimer.Stop(); + + loader->UnloadRecPoints(); + loader->UnloadTracks(); + loader->UnloadDigits(); + + delete recoEvent; + + delete recoCluster; + + 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; @@ -149,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; @@ -162,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() ; @@ -222,11 +544,14 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const // reading info from tracks recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks); - trackParam = recTrack->GetTrackParamAtVertex(); + trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First(); + + 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(); @@ -240,9 +565,9 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum); theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope)); theESDTrack->SetThetaY(TMath::ATan(bendingSlope)); - theESDTrack->SetZ(vertex[2]); - theESDTrack->SetBendingCoor(vertex[1]); // calculate vertex at ESD or Tracking level ? - theESDTrack->SetNonBendingCoor(vertex[0]); + theESDTrack->SetZ(zRec); + theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ? + theESDTrack->SetNonBendingCoor(xRec); theESDTrack->SetChi2(fitFmin); theESDTrack->SetNHit(nTrackHits); theESDTrack->SetMatchTrigger(matchTrigger); @@ -253,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); }