X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FMUONClusterInfo.C;h=b62e2eeb62c973480ac12811bbb133323cfbe3b8;hb=1e686a95afbd1f7c6a6c2a75d8fa731a218c7ea3;hp=cdbe9c86580c9d254131a832326dfb3f9fe10fdf;hpb=fcefbac40e93a473ccf625d8a59bb8d835907b40;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/MUONClusterInfo.C b/MUON/MUONClusterInfo.C index cdbe9c86580..b62e2eeb62c 100644 --- a/MUON/MUONClusterInfo.C +++ b/MUON/MUONClusterInfo.C @@ -1,17 +1,17 @@ /************************************************************************** -* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * -* * -* Author: The ALICE Off-line Project. * -* Contributors are mentioned in the code where appropriate. * -* * -* Permission to use, copy, modify and distribute this software and its * -* documentation strictly for non-commercial purposes is hereby granted * -* without fee, provided that the above copyright notice appears in all * -* copies and that both the copyright notice and this permission notice * -* appear in the supporting documentation. The authors make no claims * -* about the suitability of this software for any purpose. It is * -* provided "as is" without express or implied warranty. * -**************************************************************************/ + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ /* $Id$ */ @@ -24,73 +24,119 @@ #if !defined(__CINT__) || defined(__MAKECINT__) #include #include -#include +#include #include #include #include -#include #include #include +#include // STEER includes -#include "AliMagFMaps.h" -#include "AliTracker.h" #include "AliESDEvent.h" -#include "AliESDMuonTrack.h" #include "AliRecoParam.h" #include "AliCDBManager.h" -#include "AliGeomManager.h" +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliHeader.h" // MUON includes -#include "AliMpCDB.h" +#include "AliMpConstants.h" #include "AliMpSegmentation.h" #include "AliMpVSegmentation.h" #include "AliMpPad.h" +#include "AliMUONCDB.h" #include "AliMUONCalibrationData.h" #include "AliMUONVCalibParam.h" #include "AliMUONPadInfo.h" #include "AliMUONClusterInfo.h" #include "AliMUONRecoParam.h" #include "AliMUONESDInterface.h" -#include "AliMUONRefitter.h" #include "AliMUONVDigit.h" #include "AliMUONVDigitStore.h" #include "AliMUONVCluster.h" #include "AliMUONVClusterStore.h" #include "AliMUONTrack.h" -#include "AliMUONVTrackStore.h" #include "AliMUONTrackParam.h" -#include "AliMUONTrackExtrap.h" #endif const Int_t printLevel = 1; -void Prepare(); TTree* GetESDTree(TFile *esdFile); +UInt_t buildClusterMap(AliMUONTrack &track); //----------------------------------------------------------------------- -void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", const char* outFileName = "clusterInfo.root") +void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", const char* inFileName = "galice.root", + const TString ocdbPath = "local://$ALICE_ROOT/OCDB", const char* outFileName = "clusterInfo.root") { - /// refit ESD tracks from ESD clusters to get track parameters at each attached cluster; - /// fill AliMUONESDClusterInfo object with ESD cluster + corresponding track parameters; - /// write results in a new root file. + /// 1) if (esdFileName != "") + /// loop over ESD event and fill AliMUONClusterInfo object with cluster + corresponding track parameters; + /// 2) if (inFileName != "") + /// loop over RecPoints and fill AliMUONClusterInfo object with cluster not attached to a track; + /// 3) write results in a new root file. + /// + /// ******************************************* WARNING ******************************************* /// + /// track parameters at each cluster are recomputed by the interface using Kalman filter + Smoother /// + /// (It can be changed by resetting the tracker in the interface with a new recoParam object) /// + /// and the magnetic field set in the function prepare() /// + /// ******************************************* WARNING ******************************************* /// + + Bool_t useESD = (strcmp(esdFileName,"")); + Bool_t useRecPoints = (strcmp(inFileName,"")); + if (!useESD && !useRecPoints) { + Error("MUONClusterInfo","you must provide ESD and/or galice root file(s)"); + return; + } AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo(); AliMUONPadInfo padInfo; AliMUONCalibrationData* calibData = 0x0; - - // prepare the refitting - gRandom->SetSeed(0); - Prepare(); AliMUONESDInterface esdInterface; - AliMUONRefitter refitter; - refitter.Connect(&esdInterface); + AliMUONVClusterStore* clusterStore = 0x0; + AliMUONVDigitStore* digitStore = 0x0; // open the ESD file and tree and connect the ESD event - TFile* esdFile = TFile::Open(esdFileName); - TTree* esdTree = GetESDTree(esdFile); - AliESDEvent* esd = new AliESDEvent(); - esd->ReadFromTree(esdTree); + TFile* esdFile = 0x0; + TTree* esdTree = 0x0; + AliESDEvent* esd = 0x0; + Int_t runNumber = -1; + if (useESD) { + esdFile = TFile::Open(esdFileName); + esdTree = GetESDTree(esdFile); + esd = new AliESDEvent(); + esd->ReadFromTree(esdTree); + if (esdTree->GetEvent(0) <= 0) { + Error("MUONClusterInfo", "no ESD object found for event 0"); + return; + } + runNumber = esd->GetRunNumber(); + } + + // get the cluster from RecPoints + AliRunLoader * rl = 0x0; + AliLoader* MUONLoader = 0x0; + if (useRecPoints) { + rl = AliRunLoader::Open(inFileName,"MUONLoader"); + MUONLoader = rl->GetDetectorLoader("MUON"); + MUONLoader->LoadRecPoints("READ"); + MUONLoader->LoadDigits("READ"); + rl->LoadHeader(); + if (runNumber < 0) runNumber = rl->GetHeader()->GetRun(); + } + + // load necessary data from OCDB + AliCDBManager::Instance()->SetDefaultStorage(ocdbPath); + AliCDBManager::Instance()->SetRun(runNumber); + if (!AliMUONCDB::LoadField()) return; + if (!AliMUONCDB::LoadMapping(kTRUE)) return; + AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam(); + if (!recoParam) return; + + // reset tracker for track restoring initial track parameters at cluster + AliMUONESDInterface::ResetTracker(recoParam); + + // prepare access to calibration data + calibData = new AliMUONCalibrationData(runNumber); // prepare the output tree gROOT->cd(); @@ -103,50 +149,157 @@ void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root // timer start... TStopwatch timer; - // Loop over ESD events - if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries()); - else nevents = (Int_t)esdTree->GetEntries(); + // Loop over events + if (useESD) { + if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries()); + else nevents = (Int_t)esdTree->GetEntries(); + } else { + if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)rl->GetNumberOfEvents()); + else nevents = (Int_t)rl->GetNumberOfEvents(); + } for (Int_t iEvent = 0; iEvent < nevents; iEvent++) { //----------------------------------------------// // -------------- process event --------------- // //----------------------------------------------// // get the ESD of current event - esdTree->GetEvent(iEvent); - if (!esd) { - Error("CheckESD", "no ESD object found for event %d", iEvent); - return; + if (useESD) { + esdTree->GetEvent(iEvent); + if (!esd) { + Error("MUONClusterInfo", "no ESD object found for event %d", iEvent); + return; + } + // load the current esd event + esdInterface.LoadEvent(*esd); + // get digit store + if (!useRecPoints) digitStore = esdInterface.GetDigits(); } - Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks(); - if (nTracks < 1) continue; - - // prepare access to calibration data - if (!calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber()); - - // load the current event - esdInterface.LoadEvent(*esd); - - // get digit store - AliMUONVDigitStore* digitStore = esdInterface.GetDigits(); - // refit the tracks from clusters - AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromClusters(); + if (useRecPoints) { + if (!(rl->GetEvent(iEvent) == 0)) { + Error("MUONClusterInfo", "unable to load event %d", iEvent); + return; + } + // get the clusters of current event + TTree* treeR = MUONLoader->TreeR(); + clusterStore = AliMUONVClusterStore::Create(*treeR); + if ( clusterStore != 0x0 ) { + clusterStore->Clear(); + clusterStore->Connect(*treeR); + treeR->GetEvent(0); + } + // get the digits of current event + TTree* treeD = MUONLoader->TreeD(); + digitStore = AliMUONVDigitStore::Create(*treeD); + if ( digitStore != 0x0 ) { + digitStore->Clear(); + digitStore->Connect(*treeD); + treeD->GetEvent(0); + } + } //----------------------------------------------// // ------------- fill cluster info ------------ // //----------------------------------------------// - // loop over the refitted tracks - TIter nextTrack(newTrackStore->CreateIterator()); - AliMUONTrack* newTrack; - while ((newTrack = static_cast(nextTrack()))) { + // --- loop over the refitted tracks --- + if (useESD) { + + TIter nextTrack(esdInterface.CreateTrackIterator()); + AliMUONTrack* track; + while ((track = static_cast(nextTrack()))) { + + UInt_t muonClusterMap = buildClusterMap(*track); + + // loop over clusters + AliMUONTrackParam* trackParam = static_cast(track->GetTrackParamAtCluster()->First()); + while (trackParam) { + clusterInfo->Clear("C"); + + // fill cluster info + AliMUONVCluster* cluster = trackParam->GetClusterPtr(); + clusterInfo->SetRunId(esd->GetRunNumber()); + clusterInfo->SetEventId(iEvent); + clusterInfo->SetZ(cluster->GetZ()); + clusterInfo->SetClusterId(cluster->GetUniqueID()); + clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY()); + clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY()); + clusterInfo->SetClusterChi2(cluster->GetChi2()); + clusterInfo->SetClusterCharge(cluster->GetCharge()); + + // fill track info + clusterInfo->SetTrackId(track->GetUniqueID()); + clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor()); + clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope())); + clusterInfo->SetTrackP(trackParam->P()); + const TMatrixD paramCov = trackParam->GetCovariances(); + clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2))); + clusterInfo->SetTrackChi2(track->GetNormalizedChi2()); + clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge()); + clusterInfo->SetTrackNHits(track->GetNClusters()); + clusterInfo->SetTrackChamberHitMap(muonClusterMap); + + // fill pad info if available + for (Int_t i=0; iGetNDigits(); i++) { + AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i)); + if (!digit) continue; + + // pad location + const AliMpVSegmentation* seg = AliMpSegmentation::Instance()-> + GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); + AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY()); + + // calibration parameters + AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId()); + AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId()); + Int_t manuChannel = digit->ManuChannel(); + Int_t planeType = 0; + if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) { + planeType = 1; + } + + // fill pad info + padInfo.SetPadId(digit->GetUniqueID()); + padInfo.SetPadPlaneType(planeType); + padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY()); + padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY()); + padInfo.SetPadCharge((Double_t)digit->Charge()); + padInfo.SetPadADC(digit->ADC()); + padInfo.SetSaturated(digit->IsSaturated()); + padInfo.SetCalibrated(digit->IsCalibrated()); + padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1)); + padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1), + gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3)); + + clusterInfo->AddPad(padInfo); + } + + // remove clusters attached to a track + if (useRecPoints) { + AliMUONVCluster* cl = clusterStore->FindObject(cluster->GetUniqueID()); + if (cl) clusterStore->Remove(*cl); + } + + // fill cluster info tree + clusterInfoTree->Fill(); + + trackParam = static_cast(track->GetTrackParamAtCluster()->After(trackParam)); + } + + } - // loop over clusters - AliMUONTrackParam* trackParam = static_cast(newTrack->GetTrackParamAtCluster()->First()); - while (trackParam) { + } + + // --- loop over clusters not attached to a track --- + if (useRecPoints) { + + TIter nextCluster(clusterStore->CreateIterator()); + AliMUONVCluster *cluster; + while ( ( cluster = static_cast(nextCluster()) ) ) { + clusterInfo->Clear("C"); // fill cluster info - AliMUONVCluster* cluster = trackParam->GetClusterPtr(); + clusterInfo->SetRunId(rl->GetRunNumber()); clusterInfo->SetEventId(iEvent); clusterInfo->SetZ(cluster->GetZ()); clusterInfo->SetClusterId(cluster->GetUniqueID()); @@ -155,16 +308,17 @@ void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root clusterInfo->SetClusterChi2(cluster->GetChi2()); clusterInfo->SetClusterCharge(cluster->GetCharge()); - // fill track info - clusterInfo->SetTrackId(newTrack->GetUniqueID()); - clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor()); - clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetBendingSlope()), TMath::ATan(trackParam->GetNonBendingSlope())); - clusterInfo->SetTrackP(trackParam->P()); - const TMatrixD paramCov = trackParam->GetCovariances(); - clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2))); - clusterInfo->SetTrackChi2(newTrack->GetNormalizedChi2()); - clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge()); - + // fill dummy track info + clusterInfo->SetTrackId(0); + clusterInfo->SetTrackXY(0.,0.); + clusterInfo->SetTrackThetaXY(0.,0.); + clusterInfo->SetTrackP(0.); + clusterInfo->SetTrackXYErr(0.,0.); + clusterInfo->SetTrackChi2(0.); + clusterInfo->SetTrackCharge(0); + clusterInfo->SetTrackNHits(0); + clusterInfo->SetTrackChamberHitMap(0); + // fill pad info if available for (Int_t i=0; iGetNDigits(); i++) { AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i)); @@ -172,18 +326,23 @@ void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root // pad location const AliMpVSegmentation* seg = AliMpSegmentation::Instance()-> - GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); - AliMpPad pad = seg->PadByIndices(AliMpIntPair(digit->PadX(), digit->PadY())); + GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); + AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY()); // calibration parameters AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId()); AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId()); Int_t manuChannel = digit->ManuChannel(); + Int_t planeType = 0; + if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) { + planeType = 1; + } // fill pad info padInfo.SetPadId(digit->GetUniqueID()); - padInfo.SetPadXY(pad.Position().X(), pad.Position().Y()); - padInfo.SetPadDimXY(pad.Dimensions().X(), pad.Dimensions().Y()); + padInfo.SetPadPlaneType(planeType); + padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY()); + padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY()); padInfo.SetPadCharge((Double_t)digit->Charge()); padInfo.SetPadADC(digit->ADC()); padInfo.SetSaturated(digit->IsSaturated()); @@ -198,13 +357,13 @@ void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root // fill cluster info tree clusterInfoTree->Fill(); - trackParam = static_cast(newTrack->GetTrackParamAtCluster()->After(trackParam)); } + delete digitStore; + delete clusterStore; + } - // free memory - delete newTrackStore; } // ...timer stop @@ -223,48 +382,17 @@ void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root delete calibData; printf("Deleting clusterInfo\n"); delete clusterInfo; - printf("Closing esdFile\n"); - esdFile->Close(); - printf("Deleting esd\n"); - delete esd; - // delete padInfo; - cout<UnloadDigits(); + MUONLoader->UnloadRecPoints(); + rl->UnloadHeader(); + delete rl; } - - // set mag field - printf("Loading field map...\n"); - AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG); - AliTracker::SetFieldMap(field, kFALSE); - AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap()); - - // Load mapping - AliCDBManager* man = AliCDBManager::Instance(); - man->SetDefaultStorage("local://$ALICE_ROOT"); - man->SetRun(0); - if ( ! AliMpCDB::LoadDDLStore() ) { - Error("MUONRefit","Could not access mapping from OCDB !"); - exit(-1); + if (useESD) { + esdFile->Close(); + delete esd; } - - // set reconstruction parameters - AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam(); - muonRecoParam->Print("FULL"); - AliRecoParam::Instance()->RegisterRecoParam(muonRecoParam); - + cout<(track.GetTrackParamAtCluster()->First()); + while (trackParam) { + + muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId()); + + trackParam = static_cast(track.GetTrackParamAtCluster()->After(trackParam)); + } + + return muonClusterMap; + +} +