From 00612ffc94f31c6bf1806a2607acd98a005a24ff Mon Sep 17 00:00:00 2001 From: martinez Date: Mon, 22 Nov 2010 06:40:41 +0000 Subject: [PATCH] =?utf8?q?-=20Remove=20the=20holes=20in=20the=20lists=20of?= =?utf8?q?=20output=20histograms=20to=20make=20the=20merging=20working=20w?= =?utf8?q?ith=20the=20new=20root=20version.=20Add=20a=20new=20compiled=20m?= =?utf8?q?acro=20(MuonResolution.C)=20to=20run=20the=20task.=20Add=20the?= =?utf8?q?=20possibility=20to=20fit=20the=20residuals=20to=20extract=20the?= =?utf8?q?=20resolution=20(now=20the=20default=20option=20instead=20of=20t?= =?utf8?q?aking=20the=20RMS).=20Add=20the=20possibility=20to=20select=20ev?= =?utf8?q?ents=20of=20a=20given=20trigger=20class=20(MUON,=20MB,=20?= =?utf8?q?=E2=80=A6).=20Add=20the=20possibility=20to=20select=20tracks=20i?= =?utf8?q?n=20the=20geometrical=20acceptance=20(passing=20Rabs=20and=20eta?= =?utf8?q?=20cuts).=20Add=20mchview=20objects=20to=20the=20output=20file?= =?utf8?q?=20to=20display=20DE=20resolution=20and=20systematic=20shifts=20?= =?utf8?q?(Philippe=20P.)?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit --- PWG3/muondep/AddTaskMuonResolution.C | 8 +- .../muondep/AliAnalysisTaskMuonResolution.cxx | 267 +++++++--- PWG3/muondep/AliAnalysisTaskMuonResolution.h | 156 ++++-- PWG3/muondep/MuonResolution.C | 504 ++++++++++++++++++ PWG3/muondep/RunMuonResolution.C | 406 ++------------ 5 files changed, 875 insertions(+), 466 deletions(-) create mode 100644 PWG3/muondep/MuonResolution.C diff --git a/PWG3/muondep/AddTaskMuonResolution.C b/PWG3/muondep/AddTaskMuonResolution.C index 9fb473d3f62..69f0dc9d242 100644 --- a/PWG3/muondep/AddTaskMuonResolution.C +++ b/PWG3/muondep/AddTaskMuonResolution.C @@ -1,5 +1,7 @@ -AliAnalysisTaskMuonResolution *AddTaskMuonResolution(Bool_t selectPhysics = kFALSE, Bool_t matchTrig = kFALSE, Double_t minMomentum = 0., - Bool_t correctForSystematics = kTRUE, Int_t extrapMode = 1) +AliAnalysisTaskMuonResolution *AddTaskMuonResolution(Bool_t selectPhysics = kFALSE, Bool_t selectTrigger = kFALSE, + Bool_t matchTrig = kTRUE, Bool_t applyAccCut = kTRUE, + Double_t minMomentum = 0., Bool_t correctForSystematics = kTRUE, + Int_t extrapMode = 1) { /// Add AliAnalysisTaskMuonResolution to the train (Philippe Pillot) @@ -25,7 +27,9 @@ AliAnalysisTaskMuonResolution *AddTaskMuonResolution(Bool_t selectPhysics = kFAL return NULL; } task->SelectPhysics(selectPhysics); + task->SelectTrigger(selectTrigger); task->MatchTrigger(matchTrig); + task->ApplyAccCut(applyAccCut); task->SetMinMomentum(minMomentum); task->CorrectForSystematics(correctForSystematics); task->SetExtrapMode(extrapMode); diff --git a/PWG3/muondep/AliAnalysisTaskMuonResolution.cxx b/PWG3/muondep/AliAnalysisTaskMuonResolution.cxx index 93a4681334f..1bc40dd563b 100644 --- a/PWG3/muondep/AliAnalysisTaskMuonResolution.cxx +++ b/PWG3/muondep/AliAnalysisTaskMuonResolution.cxx @@ -25,6 +25,9 @@ #include #include #include +#include +#include +#include // STEER includes #include "AliESDEvent.h" @@ -34,7 +37,6 @@ #include "AliGeomManager.h" // ANALYSIS includes -#include "AliAnalysisTaskSE.h" #include "AliAnalysisDataSlot.h" #include "AliAnalysisManager.h" #include "AliInputEventHandler.h" @@ -71,10 +73,18 @@ AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution() : fChamberRes(NULL), fTrackRes(NULL), fCanvases(NULL), + fDefaultStorage(""), fNEvents(0), + fShowProgressBar(kFALSE), + fPrintClResPerCh(kFALSE), + fPrintClResPerDE(kFALSE), + fGaus(NULL), fMinMomentum(0.), fSelectPhysics(kFALSE), fMatchTrig(kFALSE), + fApplyAccCut(kFALSE), + fSelectTrigger(kFALSE), + fTriggerMask(0), fExtrapMode(1), fCorrectForSystematics(kTRUE), fOCDBLoaded(kFALSE), @@ -83,7 +93,8 @@ AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution() : fOldAlignStorage(""), fNewAlignStorage(""), fOldGeoTransformer(NULL), - fNewGeoTransformer(NULL) + fNewGeoTransformer(NULL), + fSelectTriggerClass(NULL) { /// Default constructor } @@ -97,10 +108,18 @@ AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) : fChamberRes(NULL), fTrackRes(NULL), fCanvases(NULL), + fDefaultStorage("raw://"), fNEvents(0), + fShowProgressBar(kFALSE), + fPrintClResPerCh(kFALSE), + fPrintClResPerDE(kFALSE), + fGaus(NULL), fMinMomentum(0.), fSelectPhysics(kFALSE), fMatchTrig(kFALSE), + fApplyAccCut(kFALSE), + fSelectTrigger(kFALSE), + fTriggerMask(0), fExtrapMode(1), fCorrectForSystematics(kTRUE), fOCDBLoaded(kFALSE), @@ -109,12 +128,15 @@ AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) : fOldAlignStorage(""), fNewAlignStorage(""), fOldGeoTransformer(NULL), - fNewGeoTransformer(NULL) + fNewGeoTransformer(NULL), + fSelectTriggerClass(NULL) { /// Constructor for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) SetStartingResolution(i, -1., -1.); + FitResiduals(); + // Output slot #1 writes into a TObjArray container DefineOutput(1,TObjArray::Class()); // Output slot #2 writes into a TObjArray container @@ -131,13 +153,17 @@ AliAnalysisTaskMuonResolution::AliAnalysisTaskMuonResolution(const char *name) : AliAnalysisTaskMuonResolution::~AliAnalysisTaskMuonResolution() { /// Destructor - SafeDelete(fResiduals); - SafeDelete(fResidualsVsP); + if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { + SafeDelete(fResiduals); + SafeDelete(fResidualsVsP); + SafeDelete(fTrackRes); + } SafeDelete(fChamberRes); - SafeDelete(fTrackRes); SafeDelete(fCanvases); + SafeDelete(fGaus); SafeDelete(fOldGeoTransformer); SafeDelete(fNewGeoTransformer); + SafeDelete(fSelectTriggerClass); } //___________________________________________________________________________ @@ -148,6 +174,15 @@ void AliAnalysisTaskMuonResolution::UserCreateOutputObjects() // do it once the OCDB has been loaded (i.e. from NotifyRun()) if (!fOCDBLoaded) return; + // set the list of trigger classes that can be selected to fill histograms (in case the physics selection is not used) + fSelectTriggerClass = new TList(); + fSelectTriggerClass->SetOwner(); + fSelectTriggerClass->AddLast(new TObjString(" CINT1B-ABCE-NOPF-ALL ")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMB); + fSelectTriggerClass->AddLast(new TObjString(" CMUS1B-ABCE-NOPF-MUON ")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMUON); + fSelectTriggerClass->AddLast(new TObjString(" CINT1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMB); + fSelectTriggerClass->AddLast(new TObjString(" CMUS1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kMUON); + fSelectTriggerClass->AddLast(new TObjString(" CSH1-B-")); fSelectTriggerClass->Last()->SetUniqueID(AliVEvent::kHighMult); + fResiduals = new TObjArray(1000); fResiduals->SetOwner(); fResidualsVsP = new TObjArray(1000); @@ -164,7 +199,7 @@ void AliAnalysisTaskMuonResolution::UserCreateOutputObjects() if (recoParam->GetDefaultBendingReso(i) > maxSigma[1]) maxSigma[1] = recoParam->GetDefaultBendingReso(i); } const char* axes[2] = {"X", "Y"}; - const Int_t nBins = 2000; + const Int_t nBins = 5000; const Int_t nSigma = 10; const Int_t pNBins = 20; const Double_t pEdges[2] = {0., 50.}; @@ -266,7 +301,16 @@ void AliAnalysisTaskMuonResolution::UserExec(Option_t *) AliESDEvent* esd = dynamic_cast(InputEvent()); if (!esd) return; - if ((++fNEvents)%100 == 0) cout<<"\rEvent processing... "<IsEventSelected() : 0; + if (fSelectPhysics && triggerWord == 0) return; + + // skip events that do not pass the trigger selection if required + TString FiredTriggerClasses = esd->GetFiredTriggerClasses(); + if (!fSelectPhysics) triggerWord = BuildTriggerWord(FiredTriggerClasses); + if (fSelectTrigger && (triggerWord & fTriggerMask) == 0) return; // get tracker to refit AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker(); @@ -281,12 +325,14 @@ void AliAnalysisTaskMuonResolution::UserExec(Option_t *) // skip ghost tracks if (!esdTrack->ContainTrackerData()) continue; - // skip tracks that do not pass the physics selection if required - if (fSelectPhysics && fInputHandler && !fInputHandler->IsEventSelected()) continue; - // skip tracks not matched with trigger if required if (fMatchTrig && !esdTrack->ContainTriggerData()) continue; + // skip tracks that do not pass the acceptance cuts if required + Double_t thetaAbs = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg(); + Double_t eta = esdTrack->Eta(); + if (fApplyAccCut && (thetaAbs < 2. || thetaAbs > 9. || eta < -4. || eta > -2.5)) continue; + // skip low momentum tracks if (esdTrack->PUncorrected() < fMinMomentum) continue; @@ -494,6 +540,7 @@ void AliAnalysisTaskMuonResolution::NotifyRun() if (fOCDBLoaded) return; AliCDBManager* cdbm = AliCDBManager::Instance(); + cdbm->SetDefaultStorage(fDefaultStorage.Data()); cdbm->SetRun(fCurrentRunNumber); if (!AliMUONCDB::LoadField()) return; @@ -503,8 +550,6 @@ void AliAnalysisTaskMuonResolution::NotifyRun() AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam(); if (!recoParam) return; - fOCDBLoaded = kTRUE; - AliMUONESDInterface::ResetTracker(recoParam); for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { @@ -527,7 +572,7 @@ void AliAnalysisTaskMuonResolution::NotifyRun() if (fReAlign) { - // recover default storage name + // recover default storage full name (raw:// cannot be used to set specific storage) TString defaultStorage(cdbm->GetDefaultStorage()->GetType()); if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data()); else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data()); @@ -570,13 +615,17 @@ void AliAnalysisTaskMuonResolution::NotifyRun() } - // print starting chamber resolution - printf("\nstarting chamber resolution:\n"); - printf(" - non-bending:"); - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]); - printf("\n - bending:"); - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]); - printf("\n\n"); + // print starting chamber resolution if required + if (fPrintClResPerCh) { + printf("\nstarting chamber resolution:\n"); + printf(" - non-bending:"); + for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",fClusterResNB[i]); + printf("\n - bending:"); + for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",fClusterResB[i]); + printf("\n\n"); + } + + fOCDBLoaded = kTRUE; UserCreateOutputObjects(); @@ -734,7 +783,7 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *) fLocalChi2->AddAtAndExpand(g, kLocalChi2PerDEMean+ia); // compute residual mean and dispersion and averaged local chi2 per chamber and half chamber - Double_t meanIn, meanInErr, meanOut, meanOutErr, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr; + Double_t meanIn, meanInErr, meanOut, meanOutErr, sigma, sigmaIn, sigmaInErr, sigmaOut, sigmaOutErr; Double_t sigmaTrack, sigmaTrackErr, sigmaMCS, sigmaMCSErr, clusterRes, clusterResErr, sigmaCluster, sigmaClusterErr; for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { @@ -750,10 +799,12 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *) delete tmp; if (fCorrectForSystematics) { - sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn); - sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.; - sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut); - sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.; + sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn); + sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.; + sigmaIn = sigma; + sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut); + sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.; + sigmaOut = sigma; } ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia))->SetPoint(i, i+1, sigmaOut); ((TGraphErrors*)fChamberRes->UncheckedAt(kResidualPerChDispersion_ClusterOut+ia))->SetPointError(i, 0., sigmaOutErr); @@ -768,11 +819,11 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *) // method 2 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e"); - GetMean(tmp, sigmaTrack, sigmaTrackErr, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), i, i+1, kFALSE); + GetMean(tmp, sigmaTrack, sigmaTrackErr, (TGraphErrors*)fChamberRes->UncheckedAt(kTrackResPerChMean+ia), i, i+1, kFALSE, kFALSE); delete tmp; tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerCh+ia))->ProjectionY("tmp",i+1,i+1,"e"); - GetMean(tmp, sigmaMCS, sigmaMCSErr, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), i, i+1, kFALSE); + GetMean(tmp, sigmaMCS, sigmaMCSErr, (TGraphErrors*)fChamberRes->UncheckedAt(kMCSPerChMean+ia), i, i+1, kFALSE, kFALSE); delete tmp; sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack; @@ -820,10 +871,12 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *) delete tmp; if (fCorrectForSystematics) { - sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn); - sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.; - sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut); - sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.; + sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn); + sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.; + sigmaIn = sigma; + sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut); + sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.; + sigmaOut = sigma; } clusterRes = TMath::Sqrt(sigmaIn*sigmaOut); @@ -834,11 +887,11 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *) // method 2 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e"); - GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE); + GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE, kFALSE); delete tmp; tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerHalfCh+ia))->ProjectionY("tmp",k+1,k+1,"e"); - GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE); + GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE, kFALSE); delete tmp; sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack; @@ -877,10 +930,12 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *) delete tmp; if (fCorrectForSystematics) { - sigmaIn = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn); - sigmaInErr = (sigmaIn>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigmaIn : 0.; - sigmaOut = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut); - sigmaOutErr = (sigmaOut>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigmaOut : 0.; + sigma = TMath::Sqrt(sigmaIn*sigmaIn + meanIn*meanIn); + sigmaInErr = (sigma>0) ? TMath::Sqrt(sigmaIn*sigmaIn*sigmaInErr*sigmaInErr + meanIn*meanIn*meanInErr*meanInErr) / sigma : 0.; + sigmaIn = sigma; + sigma = TMath::Sqrt(sigmaOut*sigmaOut + meanOut*meanOut); + sigmaOutErr = (sigma>0) ? TMath::Sqrt(sigmaOut*sigmaOut*sigmaOutErr*sigmaOutErr + meanOut*meanOut*meanOutErr*meanOutErr) / sigma : 0.; + sigmaOut = sigma; } clusterRes = TMath::Sqrt(sigmaIn*sigmaOut); @@ -891,11 +946,11 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *) // method 2 tmp = ((TH2F*)fResiduals->UncheckedAt(kTrackResPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e"); - GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE); + GetMean(tmp, sigmaTrack, sigmaTrackErr, 0x0, 0, 0, kFALSE, kFALSE); delete tmp; tmp = ((TH2F*)fResiduals->UncheckedAt(kMCSPerDE+ia))->ProjectionY("tmp",i+1,i+1,"e"); - GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE); + GetMean(tmp, sigmaMCS, sigmaMCSErr, 0x0, 0, 0, kFALSE, kFALSE); delete tmp; sigmaCluster = sigmaOut*sigmaOut - sigmaTrack*sigmaTrack; @@ -1112,12 +1167,30 @@ void AliAnalysisTaskMuonResolution::Terminate(Option_t *) fCanvases->AddAtAndExpand(cResPerChVsP, kResPerChVsP); // print results - printf("\nchamber resolution:\n"); - printf(" - non-bending:"); - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]); - printf("\n - bending:"); - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]); - printf("\n\n"); + if (fPrintClResPerCh) { + printf("\nchamber resolution:\n"); + printf(" - non-bending:"); + for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",newClusterRes[0][i]); + printf("\n - bending:"); + for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",newClusterRes[1][i]); + printf("\n\n"); + } + + if (fPrintClResPerDE) { + Double_t iDE, clRes; + printf("\nDE resolution:\n"); + printf(" - non-bending:"); + for (Int_t i = 0; i < fNDE; i++) { + ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma))->GetPoint(i, iDE, clRes); + printf((i==0)?" %5.3f":", %5.3f", clRes); + } + printf("\n - bending:"); + for (Int_t i = 0; i < fNDE; i++) { + ((TGraphErrors*)fChamberRes->UncheckedAt(kCombinedResidualPerDESigma+1))->GetPoint(i, iDE, clRes); + printf((i==0)?" %6.4f":", %6.4f", clRes); + } + printf("\n\n"); + } // Post final data. PostData(3, fLocalChi2); @@ -1202,35 +1275,85 @@ void AliAnalysisTaskMuonResolution::ZoomRight(TH1* h, Double_t fractionCut) } //________________________________________________________________________ -void AliAnalysisTaskMuonResolution::GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom) +void AliAnalysisTaskMuonResolution::GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom, Bool_t enableFit) { - /// Fill graph with the mean value of the histogram and the corresponding error (zooming if required) - Int_t firstBin = h->GetXaxis()->GetFirst(); - Int_t lastBin = h->GetXaxis()->GetLast(); - if (zoom) Zoom(h); - mean = (h->GetEntries() > fgkMinEntries) ? h->GetMean() : 0.; - meanErr = (h->GetEntries() > fgkMinEntries) ? h->GetMeanError() : 0.; + /// Fill graph with the mean value and the corresponding error (zooming if required) + + if (h->GetEntries() < fgkMinEntries) { // not enough entries + + mean = 0.; + meanErr = 0.; + + } else if (enableFit && fGaus) { // take the mean of a gaussian fit + + fGaus->SetParameters(h->GetEntries(), 0., 0.1); + + h->Fit("fGaus", "WWNQ"); + + mean = fGaus->GetParameter(1); + meanErr = fGaus->GetParError(1); + + } else { // take the mean of the distribution + + Int_t firstBin = h->GetXaxis()->GetFirst(); + Int_t lastBin = h->GetXaxis()->GetLast(); + + if (zoom) Zoom(h); + + mean = h->GetMean(); + meanErr = h->GetMeanError(); + + if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin); + + } + + // fill graph if required if (g) { g->SetPoint(i, x, mean); g->SetPointError(i, 0., meanErr); } - if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin); + } //________________________________________________________________________ void AliAnalysisTaskMuonResolution::GetRMS(TH1* h, Double_t& rms, Double_t& rmsErr, TGraphErrors* g, Int_t i, Double_t x, Bool_t zoom) { - /// Return the RMS of the histogram and the corresponding error (zooming if required) and fill graph if !=0x0 - Int_t firstBin = h->GetXaxis()->GetFirst(); - Int_t lastBin = h->GetXaxis()->GetLast(); - if (zoom) Zoom(h); - rms = (h->GetEntries() > fgkMinEntries) ? h->GetRMS() : 0.; - rmsErr = (h->GetEntries() > fgkMinEntries) ? h->GetRMSError() : 0.; + /// Return the dispersion value and the corresponding error (zooming if required) and fill graph if !=0x0 + + if (h->GetEntries() < fgkMinEntries) { // not enough entries + + rms = 0.; + rmsErr = 0.; + + } else if (fGaus) { // take the sigma of a gaussian fit + + fGaus->SetParameters(h->GetEntries(), 0., 0.1); + + h->Fit("fGaus", "WWNQ"); + + rms = fGaus->GetParameter(2); + rmsErr = fGaus->GetParError(2); + + } else { // take the RMS of the distribution + + Int_t firstBin = h->GetXaxis()->GetFirst(); + Int_t lastBin = h->GetXaxis()->GetLast(); + + if (zoom) Zoom(h); + + rms = h->GetRMS(); + rmsErr = h->GetRMSError(); + + if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin); + + } + + // fill graph if required if (g) { g->SetPoint(i, x, rms); g->SetPointError(i, 0., rmsErr); } - if (zoom) h->GetXaxis()->SetRange(firstBin,lastBin); + } //________________________________________________________________________ @@ -1248,7 +1371,8 @@ void AliAnalysisTaskMuonResolution::FillSigmaClusterVsP(TH2* hIn, TH2* hOut, TGr Double_t p = 0.5 * (hIn->GetBinLowEdge(j) + hIn->GetBinLowEdge(j+1)); Double_t pErr = p - hIn->GetBinLowEdge(j); clusterRes = TMath::Sqrt(sigmaIn*sigmaOut); - clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.; + //clusterResErr = (clusterRes > 0.) ? 0.5 * TMath::Sqrt(sigmaInErr*sigmaInErr*sigmaOut*sigmaOut + sigmaIn*sigmaIn*sigmaOutErr*sigmaOutErr) / clusterRes : 0.; + clusterResErr = TMath::Sqrt(sigmaInErr*sigmaOutErr); g->SetPoint(j, p, clusterRes); g->SetPointError(j, pErr, clusterResErr); } @@ -1286,3 +1410,22 @@ void AliAnalysisTaskMuonResolution::Cov2CovP(const AliMUONTrackParam ¶m, TMa covP.Mult(jacob,tmp); } +//__________________________________________________________________________ +UInt_t AliAnalysisTaskMuonResolution::BuildTriggerWord(TString& FiredTriggerClasses) +{ + /// build the trigger word from the fired trigger classes and the list of selectable trigger + + UInt_t word = 0; + + TObjString* trigClasseName = 0x0; + TIter nextTrigger(fSelectTriggerClass); + while ((trigClasseName = static_cast(nextTrigger()))) { + + TRegexp GenericTriggerClasseName(trigClasseName->String()); + if (FiredTriggerClasses.Contains(GenericTriggerClasseName)) word |= trigClasseName->GetUniqueID(); + + } + + return word; +} + diff --git a/PWG3/muondep/AliAnalysisTaskMuonResolution.h b/PWG3/muondep/AliAnalysisTaskMuonResolution.h index 71749c967a1..6d456e502a0 100644 --- a/PWG3/muondep/AliAnalysisTaskMuonResolution.h +++ b/PWG3/muondep/AliAnalysisTaskMuonResolution.h @@ -10,12 +10,15 @@ #include #include +#include #include "AliMUONConstants.h" +#include "AliAnalysisTaskSE.h" class TH1; class TH2; class TGraphErrors; class TObjArray; +class TList; class AliMUONTrack; class AliMUONTrackParam; class AliMUONGeometryTransformer; @@ -27,6 +30,9 @@ public: AliAnalysisTaskMuonResolution(const char *name); virtual ~AliAnalysisTaskMuonResolution(); + /// Set location of the default OCDB storage (if not set use "raw://") + void SetDefaultStorage(const char* ocdbPath) { fDefaultStorage = ocdbPath; } + void SetStartingResolution(Int_t chId, Double_t valNB, Double_t valB); void SetStartingResolution(Double_t valNB[10], Double_t valB[10]); void GetStartingResolution(Double_t valNB[10], Double_t valB[10]); @@ -40,6 +46,14 @@ public: /// set the flag to use only tracks matched with trigger or not void MatchTrigger(Bool_t flag = kTRUE) { fMatchTrig = flag; } + /// set the flag to use only tracks passing the acceptance cuts (Rabs, eta) + void ApplyAccCut(Bool_t flag = kTRUE) { fApplyAccCut = flag; } + + /// Select events belonging to at least one of the trigger classes selected by the mask to fill histograms: + /// - if the physics selection is used, apply the mask to the trigger word returned by the physics selection + /// - if not, apply the mask to the trigger word built by looking for triggers listed in "fSelectTriggerClass" + void SelectTrigger(Bool_t flag = kTRUE, UInt_t mask = AliVEvent::kMUON) {fSelectTrigger = flag; fTriggerMask = mask;} + /// set the extrapolation mode to get the track parameters and covariances at a given cluster: /// 0 = extrapolate from the closest cluster; 1 = extrapolate from the previous cluster except between stations 2-3-4 void SetExtrapMode(Int_t val) { fExtrapMode = val; } @@ -52,6 +66,14 @@ public: /// return the list of summary canvases TObjArray* GetCanvases() {return fCanvases;} + /// set the flag to show the progression bar + void ShowProgressBar(Bool_t flag = kTRUE) {fShowProgressBar = flag;} + + /// set the flag to print the cluster resolution per chamber/DE + void PrintClusterRes(Bool_t perCh = kTRUE, Bool_t perDE = kFALSE) {fPrintClResPerCh = perCh; fPrintClResPerDE = perDE;} + + void FitResiduals(Bool_t flag = kTRUE); + virtual void UserCreateOutputObjects(); virtual void UserExec(Option_t *); virtual void NotifyRun(); @@ -68,66 +90,78 @@ private: void Zoom(TH1* h, Double_t fractionCut = 0.01); void ZoomLeft(TH1* h, Double_t fractionCut = 0.02); void ZoomRight(TH1* h, Double_t fractionCut = 0.02); - void GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g = 0x0, Int_t i = 0, Double_t x = 0, Bool_t zoom = kTRUE); + void GetMean(TH1* h, Double_t& mean, Double_t& meanErr, TGraphErrors* g = 0x0, Int_t i = 0, Double_t x = 0, Bool_t zoom = kTRUE, Bool_t enableFit = kTRUE); void GetRMS(TH1* h, Double_t& rms, Double_t& rmsErr, TGraphErrors* g = 0x0, Int_t i = 0, Double_t x = 0, Bool_t zoom = kTRUE); void FillSigmaClusterVsP(TH2* hIn, TH2* hOut, TGraphErrors* g, Bool_t zoom = kTRUE); void Cov2CovP(const AliMUONTrackParam ¶m, TMatrixD &covP); + UInt_t BuildTriggerWord(TString& FiredTriggerClasses); private: - enum outputIndices { + enum eResiduals { kResidualPerCh_ClusterIn = 0, ///< cluster-track residual-X/Y distribution per chamber (cluster attached to the track) kResidualPerCh_ClusterOut = 2, ///< cluster-track residual-X/Y distribution per chamber (cluster not attached to the track) kTrackResPerCh = 4, ///< track resolution-X/Y per chamber kMCSPerCh = 6, ///< MCS X/Y-dispersion of extrapolated track per chamber kClusterRes2PerCh = 8, ///< cluster X/Y-resolution per chamber - kResidualInChVsP_ClusterIn = 10, ///< cluster-track residual-X/Y distribution in chamber i versus momentum (cluster attached to the track) - kResidualInChVsP_ClusterOut = 30, ///< cluster-track residual-X/Y distribution in chamber i versus momentum (cluster not attached to the track) - kResidualPerDE_ClusterIn = 50, ///< cluster-track residual-X/Y distribution per DE (cluster attached to the track) - kResidualPerDE_ClusterOut = 52, ///< cluster-track residual-X/Y distribution per DE (cluster not attached to the track) - kTrackResPerDE = 54, ///< track resolution-X/Y per DE - kMCSPerDE = 56, ///< MCS X/Y-dispersion of extrapolated track per DE - kResidualPerHalfCh_ClusterIn = 58, ///< cluster-track residual-X/Y distribution per half chamber (cluster attached to the track) - kResidualPerHalfCh_ClusterOut = 60, ///< cluster-track residual-X/Y distribution per half chamber (cluster not attached to the track) - kTrackResPerHalfCh = 62, ///< track resolution-X/Y per half chamber - kMCSPerHalfCh = 64, ///< MCS X/Y-dispersion of extrapolated track per half chamber - - kLocalChi2PerCh = 100, ///< local chi2-X/Y/total distribution per chamber - kLocalChi2PerDE = 103, ///< local chi2-X/Y/total distribution per DE - kLocalChi2PerChMean = 106, ///< local chi2-X/Y/total per chamber: mean - kLocalChi2PerDEMean = 109, ///< local chi2-X/Y/total per DE: mean - - kResidualPerChMean_ClusterIn = 150, ///< cluster-track residual-X/Y per chamber: mean (cluster in) - kResidualPerChMean_ClusterOut = 152, ///< cluster-track residual-X/Y per chamber: mean (cluster out) - kResidualPerChSigma_ClusterIn = 154, ///< cluster-track residual-X/Y per chamber: sigma (cluster in) - kResidualPerChSigma_ClusterOut = 156, ///< cluster-track residual-X/Y per chamber: sigma (cluster out) - kResidualPerChDispersion_ClusterOut = 158, ///< cluster-track residual-X/Y per chamber: dispersion (cluster out) - kCombinedResidualPerChSigma = 160, ///< combined cluster-track residual-X/Y per chamber - kCombinedResidualSigmaVsP = 162, ///< cluster X/Y-resolution per chamber versus momentum - kTrackResPerChMean = 164, ///< track X/Y-resolution per chamber - kMCSPerChMean = 166, ///< MCS X/Y-dispersion of extrapolated track per chamber - kClusterResPerCh = 168, ///< cluster X/Y-resolution per chamber - kCalcClusterResPerCh = 170, ///< calculated cluster X/Y-resolution per chamber - kResidualPerDEMean_ClusterIn = 172, ///< cluster-track residual-X/Y per DE: mean (cluster in) - kResidualPerDEMean_ClusterOut = 174, ///< cluster-track residual-X/Y per DE: mean (cluster out) - kCombinedResidualPerDESigma = 176, ///< combined cluster-track residual-X/Y per DE - kClusterResPerDE = 178, ///< cluster X/Y-resolution per DE - kResidualPerHalfChMean_ClusterIn = 180, ///< cluster-track residual-X/Y per half chamber: mean (cluster in) - kResidualPerHalfChMean_ClusterOut = 182, ///< cluster-track residual-X/Y per half chamber: mean (cluster out) - kCombinedResidualPerHalfChSigma = 184, ///< combined cluster-track residual-X/Y per half chamber - kClusterResPerHalfCh = 186, ///< cluster X/Y-resolution per half chamber - - kUncorrPRes = 250, ///< muon momentum reconstructed resolution at first cluster vs p - kPRes = 251, ///< muon momentum reconstructed resolution at vertex vs p - kUncorrPtRes = 252, ///< muon transverse momentum reconstructed resolution at first cluster vs p - kPtRes = 253, ///< muon transverse momentum reconstructed resolution at vertex vs p - kUncorrSlopeRes = 254, ///< muon slope-X/Y reconstructed resolution at first cluster vs p - kSlopeRes = 256, ///< muon slope-X/Y reconstructed resolution at vertex vs p - - kResPerCh = 300, ///< summary canvas - kResPerChVsP = 301, ///< summary canvas - kResPerDE = 302, ///< summary canvas - kResPerHalfCh = 303 ///< summary canvas + kResidualPerDE_ClusterIn = 10, ///< cluster-track residual-X/Y distribution per DE (cluster attached to the track) + kResidualPerDE_ClusterOut = 12, ///< cluster-track residual-X/Y distribution per DE (cluster not attached to the track) + kTrackResPerDE = 14, ///< track resolution-X/Y per DE + kMCSPerDE = 16, ///< MCS X/Y-dispersion of extrapolated track per DE + kResidualPerHalfCh_ClusterIn = 18, ///< cluster-track residual-X/Y distribution per half chamber (cluster attached to the track) + kResidualPerHalfCh_ClusterOut = 20, ///< cluster-track residual-X/Y distribution per half chamber (cluster not attached to the track) + kTrackResPerHalfCh = 22, ///< track resolution-X/Y per half chamber + kMCSPerHalfCh = 24, ///< MCS X/Y-dispersion of extrapolated track per half chamber + kLocalChi2PerCh = 26, ///< local chi2-X/Y/total distribution per chamber + kLocalChi2PerDE = 29 ///< local chi2-X/Y/total distribution per DE + }; + + enum eResidualsVsP { + kResidualInChVsP_ClusterIn = 0, ///< cluster-track residual-X/Y distribution in chamber i versus momentum (cluster attached to the track) + kResidualInChVsP_ClusterOut = 20 ///< cluster-track residual-X/Y distribution in chamber i versus momentum (cluster not attached to the track) + }; + + enum eLocalChi2 { + kLocalChi2PerChMean = 0, ///< local chi2-X/Y/total per chamber: mean + kLocalChi2PerDEMean = 3 ///< local chi2-X/Y/total per DE: mean + }; + + enum eChamberRes { + kResidualPerChMean_ClusterIn = 0, ///< cluster-track residual-X/Y per chamber: mean (cluster in) + kResidualPerChMean_ClusterOut = 2, ///< cluster-track residual-X/Y per chamber: mean (cluster out) + kResidualPerChSigma_ClusterIn = 4, ///< cluster-track residual-X/Y per chamber: sigma (cluster in) + kResidualPerChSigma_ClusterOut = 6, ///< cluster-track residual-X/Y per chamber: sigma (cluster out) + kResidualPerChDispersion_ClusterOut = 8, ///< cluster-track residual-X/Y per chamber: dispersion (cluster out) + kCombinedResidualPerChSigma = 10, ///< combined cluster-track residual-X/Y per chamber + kCombinedResidualSigmaVsP = 12, ///< cluster X/Y-resolution per chamber versus momentum + kTrackResPerChMean = 14, ///< track X/Y-resolution per chamber + kMCSPerChMean = 16, ///< MCS X/Y-dispersion of extrapolated track per chamber + kClusterResPerCh = 18, ///< cluster X/Y-resolution per chamber + kCalcClusterResPerCh = 20, ///< calculated cluster X/Y-resolution per chamber + kResidualPerDEMean_ClusterIn = 22, ///< cluster-track residual-X/Y per DE: mean (cluster in) + kResidualPerDEMean_ClusterOut = 24, ///< cluster-track residual-X/Y per DE: mean (cluster out) + kCombinedResidualPerDESigma = 26, ///< combined cluster-track residual-X/Y per DE + kClusterResPerDE = 28, ///< cluster X/Y-resolution per DE + kResidualPerHalfChMean_ClusterIn = 30, ///< cluster-track residual-X/Y per half chamber: mean (cluster in) + kResidualPerHalfChMean_ClusterOut = 32, ///< cluster-track residual-X/Y per half chamber: mean (cluster out) + kCombinedResidualPerHalfChSigma = 34, ///< combined cluster-track residual-X/Y per half chamber + kClusterResPerHalfCh = 36 ///< cluster X/Y-resolution per half chamber + }; + + enum eTrackRes { + kUncorrPRes = 0, ///< muon momentum reconstructed resolution at first cluster vs p + kPRes = 1, ///< muon momentum reconstructed resolution at vertex vs p + kUncorrPtRes = 2, ///< muon transverse momentum reconstructed resolution at first cluster vs p + kPtRes = 3, ///< muon transverse momentum reconstructed resolution at vertex vs p + kUncorrSlopeRes = 4, ///< muon slope-X/Y reconstructed resolution at first cluster vs p + kSlopeRes = 6 ///< muon slope-X/Y reconstructed resolution at vertex vs p + }; + + enum eCanvases { + kResPerCh = 0, ///< summary canvas + kResPerChVsP = 1, ///< summary canvas + kResPerDE = 2, ///< summary canvas + kResPerHalfCh = 3 ///< summary canvas }; static const Int_t fgkMinEntries; ///< minimum number of entries needed to compute resolution @@ -137,15 +171,23 @@ private: TObjArray* fLocalChi2; //!< List of plots related to local chi2 per chamber/DE TObjArray* fChamberRes; //!< List of plots related to chamber/DE resolution TObjArray* fTrackRes; //!< List of plots related to track resolution (p, pT, ...) - TObjArray* fCanvases; //!< List of canvases summarizing te results + TObjArray* fCanvases; //!< List of canvases summarizing the results Double_t fClusterResNB[10]; ///< cluster resolution in non-bending direction Double_t fClusterResB[10]; ///< cluster resolution in bending direction + TString fDefaultStorage; ///< location of the default OCDB storage Int_t fNEvents; //!< number of processed events + Bool_t fShowProgressBar; ///< show the progression bar + Bool_t fPrintClResPerCh; ///< print the cluster resolution per chamber + Bool_t fPrintClResPerDE; ///< print the cluster resolution per DE + TF1* fGaus; ///< gaussian function to fit the residuals Double_t fMinMomentum; ///< use only tracks with momentum higher than this value Bool_t fSelectPhysics; ///< use only tracks passing the physics selection Bool_t fMatchTrig; ///< use only tracks matched with trigger + Bool_t fApplyAccCut; ///< use only tracks passing the acceptance cuts (Rabs, eta) + Bool_t fSelectTrigger; ///< use only tracks passing the trigger selection + UInt_t fTriggerMask; ///< trigger mask to be used when selecting tracks /// extrapolation mode to get the track parameters and covariances at a given cluster: /// 0 = extrapolate from the closest cluster; 1 = extrapolate from the previous cluster except between stations 2-3-4 Int_t fExtrapMode; @@ -160,7 +202,9 @@ private: AliMUONGeometryTransformer* fOldGeoTransformer; //!< geometry transformer used to recontruct the present data AliMUONGeometryTransformer* fNewGeoTransformer; //!< new geometry transformer containing the new alignment to be applied - ClassDef(AliAnalysisTaskMuonResolution, 1); // chamber resolution analysis + TList* fSelectTriggerClass; //!< list of trigger class that can be selected to fill histograms + + ClassDef(AliAnalysisTaskMuonResolution, 2); // chamber resolution analysis }; //________________________________________________________________________ @@ -204,5 +248,15 @@ inline void AliAnalysisTaskMuonResolution::ReAlign(const char* oldAlignStorage, fReAlign = kTRUE; } +//________________________________________________________________________ +inline void AliAnalysisTaskMuonResolution::FitResiduals(Bool_t flag) +{ + /// set gaussian function to fit the residual distribution to extract the mean and the dispersion. + /// if not set: take the mean and the RMS of the distribution + if (fGaus) delete fGaus; + if (flag) fGaus = new TF1("fGaus","gaus"); + else fGaus = NULL; +} + #endif diff --git a/PWG3/muondep/MuonResolution.C b/PWG3/muondep/MuonResolution.C new file mode 100644 index 00000000000..f4d0916d208 --- /dev/null +++ b/PWG3/muondep/MuonResolution.C @@ -0,0 +1,504 @@ +//-------------------------------------------------------------------------- +// Macro compiled and launch by RunMuonResolution.C for submitting muon Resolution analysis locally or on CAF. +// See RunMuonResolution.C for more details +// +// Author: Philippe Pillot - SUBATECH Nantes +//-------------------------------------------------------------------------- + +#if !defined(__CINT__) || defined(__MAKECINT__) +// ROOT includes +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "TAxis.h" +#include "THashList.h" +#include +#include +#include + +// STEER includes +#include "AliLog.h" +#include "AliCDBManager.h" +#include "AliAnalysisManager.h" +#include "AliESDInputHandler.h" +#include "AliTagAnalysis.h" +#include "AliRunTagCuts.h" +#include "AliLHCTagCuts.h" +#include "AliDetectorTagCuts.h" +#include "AliEventTagCuts.h" +#include "AliPhysicsSelectionTask.h" +#include "AliPhysicsSelection.h" +#include "AliBackgroundSelection.h" +#include "AliAnalysisDataContainer.h" +#include "AliAnalysisTaskMuonResolution.h" + +// MUON includes +#include "AliMpCDB.h" +#include "AliMpDetElement.h" +#include "AliMpDDLStore.h" +#include "AliMUONCalibParamND.h" +#include "AliMUON2DMap.h" +#include "AliMUONTrackerData.h" +#include "AliMUONPainterDataRegistry.h" +#include "AliMUONTrackerDataWrapper.h" + +#include "AddTaskPhysicsSelection.C" +#include "AddTaskMuonResolution.C" + +#endif + +enum {kLocal, kInteractif_xml, kInteractif_ESDList, kProof}; + +void LoadAlirootOnProof(TString alirootVersion, TString& extraLibs, Int_t iStep); +AliAnalysisTaskMuonResolution* CreateAnalysisTrain(Int_t mode, Int_t iStep, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig, + Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode, + Double_t clusterResNB[10], Double_t clusterResB[10]); +Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], + Double_t clusterResNBErr[10], Double_t clusterResBErr[10]); +void AddMCHViews(TFile* file); +AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name); +Int_t GetMode(TString smode, TString input); +TChain* CreateChainFromCollection(const char *xmlfile); +TChain* CreateChainFromFile(const char *rootfile); +TChain* CreateChainFromESDList(const char *esdList); +TChain* CreateChain(Int_t mode, TString input); + +//______________________________________________________________________________ +void MuonResolution(TString smode, TString inputFileName, TString alirootVersion, Int_t nSteps, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig, + Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode, Int_t nevents, TString extraLibs) +{ + /// Compute the cluster resolution by studying cluster-track residual, deconvoluting from track resolution + + // timer start... + TStopwatch* localTimer = new TStopwatch; + + // check parameters + nSteps = TMath::Max(nSteps,1); + if (extrapMode != 0 && extrapMode != 1) { + Error("MuonResolution","incorrect extrapolation mode!"); + return; + } + + // Check runing mode + Int_t mode = GetMode(smode, inputFileName); + if(mode < 0){ + Error("MuonResolution","Please provide either an ESD root file, a list of ESDs, a collection of ESDs or a dataset."); + return; + } + + // check for old output file to removed + char remove = '\0'; + if (!gSystem->Exec("ls chamberResolution_step*[0-9].root")) { + cout<<"above files must be removed from the current directory. Delete? [y=yes, n=no] "<>remove; + if (remove == 'y') gSystem->Exec("rm -f chamberResolution_step*[0-9].root"); + else { + Error("MuonResolution","cannot proceed with these files there otherwise results will be mixed up!"); + return; + } + } + + // Create input object + TObject* inputObj = 0x0; + if (mode == kProof) inputObj = new TObjString(inputFileName); + else inputObj = CreateChain(mode, inputFileName); + if (!inputObj) return; + + // set starting chamber resolution (if -1 they will be loaded from recoParam in the task) + Double_t clusterResNB[10] ={-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.}; + Double_t clusterResB[10] ={-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.}; + Double_t clusterResNBErr[10] ={0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + Double_t clusterResBErr[10] ={0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + + // output graphs + TMultiGraph* mgClusterResXVsStep = new TMultiGraph("mgClusterResXVsStep","cluster X-resolution versus step;step;#sigma_{X} (cm)"); + TMultiGraph* mgClusterResYVsStep = new TMultiGraph("mgClusterResYVsStep","cluster Y-resolution versus step;step;#sigma_{Y} (cm)"); + TGraphErrors* clusterResXVsStep[10]; + TGraphErrors* clusterResYVsStep[10]; + for (Int_t i = 0; i < 10; i++) { + clusterResXVsStep[i] = new TGraphErrors(nSteps+1); + clusterResXVsStep[i]->SetName(Form("gResX_ch%d",i+1)); + clusterResXVsStep[i]->SetMarkerStyle(kFullDotMedium); + clusterResXVsStep[i]->SetMarkerColor(i+1+i/9); + mgClusterResXVsStep->Add(clusterResXVsStep[i],"lp"); + + clusterResYVsStep[i] = new TGraphErrors(nSteps+1); + clusterResYVsStep[i]->SetName(Form("gResY_ch%d",i+1)); + clusterResYVsStep[i]->SetMarkerStyle(kFullDotMedium); + clusterResYVsStep[i]->SetMarkerColor(i+1+i/9); + mgClusterResYVsStep->Add(clusterResYVsStep[i],"lp"); + } + + // loop over step + for (Int_t iStep = 0; iStep < nSteps; iStep++) { + cout<<"step "<InitAnalysis()) { + mgr->PrintStatus(); + if (mode == kProof) mgr->StartAnalysis("proof", Form("%s",static_cast(inputObj)->GetName()), nevents); + else mgr->StartAnalysis("local", static_cast(inputObj), nevents); + } + + // save the summary canvases and mchview display + if (muonResolution->GetCanvases()) { + TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"UPDATE"); + if (outFile && outFile->IsOpen()) { + outFile->cd(); + muonResolution->GetCanvases()->Write(); + AddMCHViews(outFile); + outFile->Close(); + } + } + + // fill graph with starting resolutions from the task at first step + if (iStep == 0) { + muonResolution->GetStartingResolution(clusterResNB, clusterResB); + for (Int_t i = 0; i < 10; i++) { + clusterResXVsStep[i]->SetPoint(0, 0, clusterResNB[i]); + clusterResXVsStep[i]->SetPointError(0, 0., clusterResNBErr[i]); + clusterResYVsStep[i]->SetPoint(0, 0, clusterResB[i]); + clusterResYVsStep[i]->SetPointError(0, 0., clusterResBErr[i]); + } + } + + // read the chamber resolution from the output file + if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr)) return; + + // fill graphs with computed resolutions + for (Int_t i = 0; i < 10; i++) { + clusterResXVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResNB[i]); + clusterResXVsStep[i]->SetPointError(iStep+1, 0., clusterResNBErr[i]); + clusterResYVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResB[i]); + clusterResYVsStep[i]->SetPointError(iStep+1, 0., clusterResBErr[i]); + } + + // clean memory + delete mgr; + TObject::SetObjectStat(kFALSE); + + } + + // copy final results in results.root file + gSystem->Exec(Form("cp chamberResolution_step%d.root results.root", nSteps-1)); + + // display convergence + TCanvas* convergence = new TCanvas("convergence","convergence"); + convergence->Divide(1,2); + convergence->cd(1); + mgClusterResXVsStep->Draw("ap"); + convergence->cd(2); + mgClusterResYVsStep->Draw("ap"); + + // save convergence plots + TFile* outFile = TFile::Open("results.root","UPDATE"); + if (!outFile || !outFile->IsOpen()) return; + outFile->cd(); + mgClusterResXVsStep->Write(); + mgClusterResYVsStep->Write(); + convergence->Write(); + outFile->Close(); + + // ...timer stop + localTimer->Stop(); + localTimer->Print(); + +} + +//______________________________________________________________________________ +void LoadAlirootOnProof(TString alirootVersion, TString& extraLibs, Int_t iStep) +{ + /// Load aliroot packages and set environment on Proof + + // set general environment and close previous session + if (iStep == 0) gEnv->SetValue("XSec.GSI.DelegProxy","2"); + else gProof->Close("s"); + + // connect + TString location = "alice-caf.cern.ch"; + TString nWorkers = ""; + if (gSystem->Getenv("alien_API_USER") == NULL) TProof::Open(location.Data(), nWorkers.Data()); + else TProof::Open(Form("%s@%s",gSystem->Getenv("alien_API_USER"), location.Data()), nWorkers.Data()); + if (!gProof) return; + + // set environment and load libraries on workers + TList* list = new TList(); + list->Add(new TNamed("ALIROOT_MODE", "")); + list->Add(new TNamed("ALIROOT_EXTRA_LIBS", extraLibs.Data())); + if (!gSystem->AccessPathName("AliAnalysisTaskMuonResolution.cxx")) + list->Add(new TNamed("ALIROOT_EXTRA_INCLUDES", "MUON:MUON/mapping")); + gProof->EnablePackage(alirootVersion.Data(), list, kTRUE); + + // compile task on workers + if (!gSystem->AccessPathName("AliAnalysisTaskMuonResolution.cxx")) + gProof->Load("AliAnalysisTaskMuonResolution.cxx++g", kTRUE); + +} + +//______________________________________________________________________________ +AliAnalysisTaskMuonResolution* CreateAnalysisTrain(Int_t mode, Int_t iStep, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig, + Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode, + Double_t clusterResNB[10], Double_t clusterResB[10]) +{ + /// create the analysis train and configure it + + // Create the analysis manager + AliAnalysisManager *mgr = new AliAnalysisManager("MuonResolutionAnalysis"); + + // ESD input handler + AliESDInputHandler* esdH = new AliESDInputHandler(); + esdH->SetReadFriends(kFALSE); + esdH->SetInactiveBranches(" FMD PHOS EMCAL Pmd Trd V0s TPC " + "Cascades Kinks CaloClusters ACORDE RawData HLT TZERO ZDC" + " Cells ACORDE Pileup"); + mgr->SetInputEventHandler(esdH); + + // event selection + if (selectPhysics) { + AliPhysicsSelectionTask* physicsSelection = AddTaskPhysicsSelection(); + if (!physicsSelection) { + Error("run","AliPhysicsSelectionTask not created!"); + return 0x0; + } + } + + // Muon Resolution analysis + TString outputFileName = Form("chamberResolution_step%d.root", iStep); + AliAnalysisManager::SetCommonFileName(outputFileName.Data()); + AliAnalysisTaskMuonResolution *muonResolution = AddTaskMuonResolution(selectPhysics, selectTrigger, matchTrig, applyAccCut, minMomentum, correctForSystematics, extrapMode); + if (!muonResolution) { + Error("run","AliAnalysisTaskMuonResolution not created!"); + return 0x0; + } + if (mode == kLocal) muonResolution->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); + if (mode != kProof) muonResolution->ShowProgressBar(); + muonResolution->PrintClusterRes(kTRUE, kTRUE); + muonResolution->SetStartingResolution(clusterResNB, clusterResB); + + return muonResolution; + +} + +//______________________________________________________________________________ +Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10]) +{ + /// read the chamber resolution from the output file + + TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ"); + + if (!outFile || !outFile->IsOpen()) { + Error("GetChamberResolution","output file does not exist!"); + return kFALSE; + } + + TObjArray* summary = static_cast(outFile->FindObjectAny("ChamberRes")); + TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0; + TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast(summary->FindObject("gCombinedResidualYPerChSigma")) : 0x0; + + if (!gCombinedResidualXPerChSigma || !gCombinedResidualYPerChSigma) { + Error("GetChamberResolution","resolution graphs do not exist!"); + return kFALSE; + } + + Double_t dummy; + for (Int_t i = 0; i < 10; i++) { + gCombinedResidualXPerChSigma->GetPoint(i, dummy, clusterResNB[i]); + gCombinedResidualYPerChSigma->GetPoint(i, dummy, clusterResB[i]); + clusterResNBErr[i] = gCombinedResidualXPerChSigma->GetErrorY(i); + clusterResBErr[i] = gCombinedResidualYPerChSigma->GetErrorY(i); + } + + outFile->Close(); + + return kTRUE; +} + +//______________________________________________________________________________ +void AddMCHViews(TFile* file) +{ + /// Get from the file the graphs containing data per DE, convert them into mchview objects and save them + + if ( ! AliMpDDLStore::Instance(false) ) + { + Warning("AddMCHViews","mapping was not loaded. Loading it from $ALICE_ROOT/OCDB"); + AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); + AliCDBManager::Instance()->SetRun(999999999); + } + + AliMpCDB::LoadAll(); + + TObjArray* summary = static_cast(file->FindObjectAny("ChamberRes")); + if (!summary) { + Error("AddMCHViews","resolution graphs do not exist!"); + return; + } + + TGraphErrors* g = 0x0; + g = static_cast(summary->FindObject("gCombinedResidualXPerDESigma")); + if (g) { + file->cd(); + ConvertGraph(*g, "resoX")->Write(); + } + + g = static_cast(summary->FindObject("gCombinedResidualYPerDESigma")); + if (g) { + file->cd(); + ConvertGraph(*g, "resoY")->Write(); + } + + g = static_cast(summary->FindObject("gResidualXPerDEMean_ClusterOut")); + if (g) { + file->cd(); + ConvertGraph(*g, "shiftX")->Write(); + } + + g = static_cast(summary->FindObject("gResidualYPerDEMean_ClusterOut")); + if (g) { + file->cd(); + ConvertGraph(*g, "shiftY")->Write(); + } +} + +//______________________________________________________________________________ +AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name) +{ + /// Convert graph containing data per DE into mchview object + + AliMUON2DMap deValues(kFALSE); + + for ( Int_t i = 0 ; i < g.GetN(); ++i ) + { + double y = g.GetY()[i]; + double ey = g.GetEY()[i]; + int detElemId; + sscanf(g.GetXaxis()->GetBinLabel(i+1),"%d",&detElemId); + + AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId); + + AliMUONVCalibParam* param = new AliMUONCalibParamND(5, 1, detElemId, 0); + + Double_t sumn = 1000.0; + Double_t sumw = sumn*y; + Double_t sumw2 = (sumn-1)*ey*ey+sumw*sumw/sumn; + + param->SetValueAsDouble(0,0,sumw); + param->SetValueAsDouble(0,1,sumw2); + param->SetValueAsDouble(0,2,sumn); + param->SetValueAsDouble(0,3,de->NofChannels()); + param->SetValueAsDouble(0,4,1); + + deValues.Add(param); + } + + AliMUONTrackerData* data = new AliMUONTrackerData(name,name,deValues,1); + data->SetDimensionName(0,name); + + return data; +} + +//______________________________________________________________________________ +Int_t GetMode(TString smode, TString input) +{ + if (smode == "local") { + if ( input.EndsWith(".xml") ) return kInteractif_xml; + else if ( input.EndsWith(".txt") ) return kInteractif_ESDList; + else if ( input.EndsWith(".root") ) return kLocal; + } else if (smode == "proof") return kProof; + return -1; +} + +//______________________________________________________________________________ +TChain* CreateChainFromCollection(const char *xmlfile) +{ + // Create a chain from the collection of tags. + if (!TGrid::Connect("alien://")) return NULL; + + TGridCollection* coll = TAlienCollection::Open(xmlfile); + if (!coll) { + Error("CreateChainFromCollection", "Cannot create the AliEn collection"); + return NULL; + } + + TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE); + AliTagAnalysis *tagAna = new AliTagAnalysis("ESD"); + tagAna->ChainGridTags(tagResult); + + AliRunTagCuts *runCuts = new AliRunTagCuts(); + AliLHCTagCuts *lhcCuts = new AliLHCTagCuts(); + AliDetectorTagCuts *detCuts = new AliDetectorTagCuts(); + AliEventTagCuts *evCuts = new AliEventTagCuts(); + + // Check if the cuts configuration file was provided + if (!gSystem->AccessPathName("ConfigureCuts.C")) + gROOT->ProcessLine(Form(".x ConfigureCuts.C((AliRunTagCuts*)%p, (AliLHCTagCuts*)%p, (AliDetectorTagCuts*)%p," + " (AliEventTagCuts*)%p)", runCuts, lhcCuts, detCuts, evCuts)); + + TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts); + if (!chain || !chain->GetNtrees()) return NULL; + chain->ls(); + return chain; +} + +//______________________________________________________________________________ +TChain* CreateChainFromFile(const char *rootfile) +{ + // Create a chain using the root file. + TChain* chain = new TChain("esdTree"); + chain->Add(rootfile); + if (!chain->GetNtrees()) return NULL; + chain->ls(); + return chain; +} + +//______________________________________________________________________________ +TChain* CreateChainFromESDList(const char *esdList) +{ + // Create a chain using tags from the run list. + TChain* chain = new TChain("esdTree"); + ifstream inFile(esdList); + TString inFileName; + if (inFile.is_open()) { + while (! inFile.eof() ) { + inFileName.ReadLine(inFile,kFALSE); + if(!inFileName.EndsWith(".root")) continue; + chain->Add(inFileName.Data()); + } + } + inFile.close(); + if (!chain->GetNtrees()) return NULL; + chain->ls(); + return chain; +} + +//______________________________________________________________________________ +TChain* CreateChain(Int_t mode, TString input) +{ + printf("*******************************\n"); + printf("*** Getting the Chain ***\n"); + printf("*******************************\n"); + if(mode == kInteractif_xml) return CreateChainFromCollection(input.Data()); + else if (mode == kInteractif_ESDList) return CreateChainFromESDList(input.Data()); + else if (mode == kLocal) return CreateChainFromFile(input.Data()); + else return NULL; +} + diff --git a/PWG3/muondep/RunMuonResolution.C b/PWG3/muondep/RunMuonResolution.C index 336b8b53af9..7725f6e4ea7 100644 --- a/PWG3/muondep/RunMuonResolution.C +++ b/PWG3/muondep/RunMuonResolution.C @@ -25,52 +25,62 @@ // - libCORRFW.so // - libPWG3muondep.so // -// The macro reads ESDs and store outputs in standard output file (AnalysisResults.root) -// It also needs to load magnetic field, mapping, geometry (+alignment), and recoonstruction parameters from the OCDB +// It also needs to load magnetic field, mapping, geometry (+alignment), and reconstruction parameters from the OCDB +// +// The task reads ESDs +// Intermediate results are stored in a file chamberResolution_step<#step>.root +// Final results are stored in the file results.root // // Author: Philippe Pillot - SUBATECH Nantes //-------------------------------------------------------------------------- -TString aliroot="VO_ALICE@AliRoot::v4-19-15-AN"; - -enum {kLocal, kInteractif_xml, kInteractif_ESDList, kProof}; - -// global declaration needed in Proof mode, certainly because of the interpretor -TStopwatch* localTimer; -TMultiGraph* mgClusterResXVsStep; -TMultiGraph* mgClusterResYVsStep; -Double_t clusterResNB[10]; -Double_t clusterResB[10]; -Double_t clusterResNBErr[10]; -Double_t clusterResBErr[10]; +void LoadAlirootLocally(TString& extraLibs); -void RunMuonResolution(TString smode = "local", TString inputFileName = "AliESDs.root", Int_t nSteps = 10, - Bool_t selectPhysics = kFALSE, Bool_t matchTrig = kFALSE, Double_t minMomentum = 0., - Bool_t correctForSystematics = kTRUE, Int_t extrapMode = 1) +//______________________________________________________________________________ +void RunMuonResolution(TString smode = "local", TString inputFileName = "AliESDs.root", + TString alirootVersion = "VO_ALICE@AliRoot::v4-20-12-AN", Int_t nSteps = 5, + Bool_t selectPhysics = kFALSE, Bool_t selectTrigger = kFALSE, Bool_t matchTrig = kTRUE, + Bool_t applyAccCut = kTRUE, Double_t minMomentum = 0., Bool_t correctForSystematics = kTRUE, + Int_t extrapMode = 1, Int_t nevents = 1234567890) { /// Compute the cluster resolution by studying cluster-track residual, deconvoluting from track resolution - /// if extrapMode == 0: extrapolate from the closest cluster - /// if extrapMode == 1: extrapolate from the previous cluster except between stations 2-3-4 - /// if correctForSystematics == kTRUE: the systematic shifts of the residuals is included in the resolution - - // timer start... - localTimer = new TStopwatch; + /// - smode = "local" or "proof" + /// - inputFileName = an ESD root file or a list of ESDs or a collection of ESDs or a dataset in proof mode + /// - alirootVersion = version of aliroot package to enable on AAF (only used in proof mode) + /// - nSteps = number of times to task is run (at each step it starts with the chamber resolution obtained in the previous one) + /// - selectPhysics : apply or not the physics selection + /// - selectTrigger : select only muon trigger events or not (the type of trigger can eventually be changed) + /// - matchTrigger : select only the tracks matching the trigger or not + /// - applyAccCut : select only the tracks passing the Rabs and the eta cut or not + /// - minMomentum : select only the tracks with a total momentum above this value + /// - if correctForSystematics == kTRUE: the systematic shifts of the residuals is included in the resolution + /// - if extrapMode == 0: extrapolate from the closest cluster + /// - if extrapMode == 1: extrapolate from the previous cluster except between stations 2-3-4 + /// - nevents = maximum number of processed events + + // Load libraries locally + TString extraLibs = "RAWDatabase:CDB:STEER:MUONcore:MUONmapping:MUONcalib:MUONgeometry:MUONtrigger:MUONraw:MUONbase:MUONrec:CORRFW:PWG3muondep"; + LoadAlirootLocally(extraLibs); + + // compile analysis macro locally + gROOT->LoadMacro("$ALICE_ROOT/PWG3/muondep/MuonResolution.C++g"); + MuonResolution(smode, inputFileName, alirootVersion, nSteps, selectPhysics, selectTrigger, matchTrig, + applyAccCut, minMomentum, correctForSystematics, extrapMode, nevents, extraLibs); - // check parameters - nSteps = TMath::Max(nSteps,1); - if (extrapMode != 0 && extrapMode != 1) { - Error("RunMuonResolution","incorrect extrapolation mode!"); - return; - } - - // Check runing mode - Int_t mode = GetMode(smode, inputFileName); - if(mode < 0){ - Error("RunMuonResolution","Please provide either an ESD root file a collection of ESDs or a dataset."); - return; - } +} + +//______________________________________________________________________________ +void LoadAlirootLocally(TString& extraLibs) +{ + /// Load libraries locally + // Load common libraries gSystem->Load("libVMC"); + gSystem->Load("libTree.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libMinuit.so"); + gSystem->Load("libXMLParser.so"); + gSystem->Load("libGui.so"); gSystem->Load("libSTEERBase"); gSystem->Load("libESD"); gSystem->Load("libAOD"); @@ -79,326 +89,20 @@ void RunMuonResolution(TString smode = "local", TString inputFileName = "AliESDs // Load additional libraries gSystem->Load("libProofPlayer"); - TString extraLibs="Physics:Minuit:XMLParser:Gui:RAWDatabase:CDB:STEER:MUONcore:MUONmapping:MUONcalib:MUONgeometry:MUONtrigger:MUONraw:MUONbase:MUONrec:CORRFW:PWG3muondep"; TObjArray* libs = extraLibs.Tokenize(":"); for (Int_t i = 0; i < libs->GetEntriesFast(); i++) gSystem->Load(Form("lib%s",static_cast(libs->UncheckedAt(i))->GetName())); + delete libs; - // check for old output file to removed - char remove = ''; - if (!gSystem->Exec("ls chamberResolution_step*[0-9].root")) { - cout<<"above files must be removed from the current directory. Delete? [y=yes, n=no] "<>remove; - if (remove == 'y') gSystem->Exec("rm -f chamberResolution_step*[0-9].root"); - else { - Error("RunMuonResolution","cannot proceed with these files there otherwise results will be mixed up!"); - return; - } - } - - // OCDB access - if(mode == kLocal) AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); - else if (mode != kProof) { - if (!TGrid::Connect("alien://")) return; - if(mode == kInteractif_ESDList || !gSystem->AccessPathName("ConfigureCuts.C")) - AliCDBManager::Instance()->SetDefaultStorage("raw://"); - } - - // Create input object - TObject* inputObj = 0x0; - if (mode == kProof) inputObj = new TObjString(inputFileName); - else inputObj = CreateChain(mode, inputFileName); - if (!inputObj) return; - - // set starting chamber resolution (if -1 they will be loaded from recoParam in the task) - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { - clusterResNB[i] = -1.; - clusterResB[i] = -1.; - clusterResNBErr[i] = 0.; - clusterResBErr[i] = 0.; - } - - // output graphs - mgClusterResXVsStep = new TMultiGraph("mgClusterResXVsStep","cluster X-resolution versus step;step;#sigma_{X} (cm)"); - mgClusterResYVsStep = new TMultiGraph("mgClusterResYVsStep","cluster Y-resolution versus step;step;#sigma_{Y} (cm)"); - TGraphErrors* gClusterResXVsStep[10]; - TGraphErrors* gClusterResYVsStep[10]; - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { - gClusterResXVsStep[i] = new TGraphErrors(nSteps+1); - gClusterResXVsStep[i]->SetName(Form("gResX_ch%d",i+1)); - gClusterResXVsStep[i]->SetMarkerStyle(kFullDotMedium); - gClusterResXVsStep[i]->SetMarkerColor(i+1+i/9); - mgClusterResXVsStep->Add(gClusterResXVsStep[i],"lp"); - - gClusterResYVsStep[i] = new TGraphErrors(nSteps+1); - gClusterResYVsStep[i]->SetName(Form("gResY_ch%d",i+1)); - gClusterResYVsStep[i]->SetMarkerStyle(kFullDotMedium); - gClusterResYVsStep[i]->SetMarkerColor(i+1+i/9); - mgClusterResYVsStep->Add(gClusterResYVsStep[i],"lp"); - } - - // loop over step - for (Int_t iStep = 0; iStep < nSteps; iStep++) { - cout<<"step "<SetValue("XSec.GSI.DelegProxy","2"); - TProof::AddEnvVar("ALIROOT_EXTRA_LIBS", extraLibs.Data()); - } else gProof->Close("s"); - - // connect - if (gSystem->Getenv("alien_API_USER") == NULL) TProof::Open("alice-caf.cern.ch","workers=20"); - else TProof::Open(Form("%s@alice-caf.cern.ch",gSystem->Getenv("alien_API_USER")),"workers=20"); - if (!gProof) return; - - // set environment and compile task on workers - gProof->EnablePackage(aliroot.Data()); - - // prepare OCDB access on workers - gProof->Exec("AliCDBManager::Instance()->SetDefaultStorage(\"raw://\")"); - - } - - // run one step - run(mode, iStep, inputObj, extrapMode, correctForSystematics, minMomentum, matchTrig, selectPhysics, clusterResNB, clusterResB); - - // fill graph with starting resolutions from the task at first step - if (iStep == 0) for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { - gClusterResXVsStep[i]->SetPoint(0, 0, clusterResNB[i]); - gClusterResXVsStep[i]->SetPointError(0, 0., clusterResNBErr[i]); - gClusterResYVsStep[i]->SetPoint(0, 0, clusterResB[i]); - gClusterResYVsStep[i]->SetPointError(0, 0., clusterResBErr[i]); - } - - // read the chamber resolution from the output file - if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr)) return; - - // fill graphs with computed resolutions - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { - gClusterResXVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResNB[i]); - gClusterResXVsStep[i]->SetPointError(iStep+1, 0., clusterResNBErr[i]); - gClusterResYVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResB[i]); - gClusterResYVsStep[i]->SetPointError(iStep+1, 0., clusterResBErr[i]); - } - - } - - // copy final results in results.root file - gSystem->Exec(Form("cp chamberResolution_step%d.root results.root", nSteps-1)); - - // display convergence - TCanvas* convergence = new TCanvas("convergence","convergence"); - convergence->Divide(1,2); - convergence->cd(1); - mgClusterResXVsStep->Draw("ap"); - convergence->cd(2); - mgClusterResYVsStep->Draw("ap"); + // Load lib for final mchview display + gSystem->Load("libMUONgraphics"); - // save convergence plots - TFile* outFile = TFile::Open("results.root","UPDATE"); - if (!outFile || !outFile->IsOpen()) return; - outFile->cd(); - mgClusterResXVsStep->Write(); - mgClusterResYVsStep->Write(); - convergence->Write(); - outFile->Close(); + // Use AliRoot includes and compile our task + gROOT->ProcessLine(".include $ALICE_ROOT/include"); + gROOT->ProcessLine(".include $ALICE_ROOT/MUON"); + gROOT->ProcessLine(".include $ALICE_ROOT/MUON/mapping"); + gROOT->ProcessLine(".include $ALICE_ROOT/ANALYSIS/macros"); + gROOT->ProcessLine(".include $ALICE_ROOT/PWG3/muondep"); - // print results - printf("\nchamber resolution:\n"); - printf(" - non-bending:"); - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",clusterResNB[i]); - printf("\n - bending:"); - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",clusterResB[i]); - printf("\n\n"); - - // ...timer stop - localTimer->Stop(); - localTimer->Print(); - -} - -//______________________________________________________________________________ -void run(Int_t mode, Int_t iStep, TObject* input, Int_t extrapMode, Bool_t correctForSystematics, - Double_t minMomentum, Bool_t matchTrig, Bool_t selectPhysics, Double_t clusterResNB[10], Double_t clusterResB[10]) -{ - /// launch the analysis with these parameters - - // Create the analysis manager - AliAnalysisManager *mgr = new AliAnalysisManager("MuonResolutionAnalysis"); - - // ESD input handler - AliESDInputHandler* esdH = new AliESDInputHandler(); - esdH->SetReadFriends(kFALSE); - mgr->SetInputEventHandler(esdH); - - // event selection - if (selectPhysics) { - gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C"); - AliPhysicsSelectionTask* physicsSelection = AddTaskPhysicsSelection(); - if(!physicsSelection) { - Error("RunMuonResolution","AliPhysicsSelectionTask not created!"); - return; - } - physicsSelection->GetPhysicsSelection()->SetUseMuonTriggers(); - } - - // Muon Resolution analysis - gROOT->LoadMacro("$ALICE_ROOT/PWG3/muondep/AddTaskMuonResolution.C"); - AliAnalysisManager::SetCommonFileName(Form("chamberResolution_step%d.root", iStep)); - AliAnalysisTaskMuonResolution* muonResolution = AddTaskMuonResolution(selectPhysics, matchTrig, minMomentum, correctForSystematics, extrapMode); - if(!muonResolution) { - Error("RunMuonResolution","AliAnalysisTaskMuonResolution not created!"); - return; - } - muonResolution->SetStartingResolution(clusterResNB, clusterResB); - - // Enable debug printouts - //mgr->SetDebugLevel(2); - - // start local analysis - if (mgr->InitAnalysis()) { - mgr->PrintStatus(); - if (mode == kProof) mgr->StartAnalysis("proof", Form("%s",static_cast(input)->GetName())); - else mgr->StartAnalysis("local", static_cast(input)); - } - - // save the summary canvases - if (muonResolution->GetCanvases()) { - TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"UPDATE"); - if (outFile && outFile->IsOpen()) { - muonResolution->GetCanvases()->Write(); - outFile->Close(); - } - } - - // return starting chamber resolution from the task - muonResolution->GetStartingResolution(clusterResNB, clusterResB); - - // clean memory - delete mgr; - TObject::SetObjectStat(kFALSE); -} - -//______________________________________________________________________________ -Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10]) -{ - /// read the chamber resolution from the output file - - TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ"); - - if (!outFile || !outFile->IsOpen()) { - Error("GetChamberResolution","output file does not exist!"); - return kFALSE; - } - - TObjArray* summary = static_cast(outFile->FindObjectAny("ChamberRes")); - TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0; - TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast(summary->FindObject("gCombinedResidualYPerChSigma")) : 0x0; - - if (!gCombinedResidualXPerChSigma || !gCombinedResidualYPerChSigma) { - Error("GetChamberResolution","resolution graphs do not exist!"); - return kFALSE; - } - - Double_t dummy; - for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { - gCombinedResidualXPerChSigma->GetPoint(i, dummy, clusterResNB[i]); - gCombinedResidualYPerChSigma->GetPoint(i, dummy, clusterResB[i]); - clusterResNBErr[i] = gCombinedResidualXPerChSigma->GetErrorY(i); - clusterResBErr[i] = gCombinedResidualYPerChSigma->GetErrorY(i); - } - - outFile->Close(); - - return kTRUE; -} - -//______________________________________________________________________________ -Int_t GetMode(TString smode, TString input) -{ - if (smode == "local") { - if ( input.EndsWith(".xml") ) return kInteractif_xml; - else if ( input.EndsWith(".txt") ) return kInteractif_ESDList; - else if ( input.EndsWith(".root") ) return kLocal; - } else if (smode == "proof") return kProof; - return -1; -} - -//______________________________________________________________________________ -TChain* CreateChainFromCollection(const char *xmlfile) -{ - // Create a chain from the collection of tags. - TAlienCollection* coll = TAlienCollection::Open(xmlfile); - if (!coll) { - Error("CreateChainFromCollection", Form("Cannot create an AliEn collection from %s!", xmlfile)); - return NULL; - } - - TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE); - AliTagAnalysis *tagAna = new AliTagAnalysis("ESD"); - tagAna->ChainGridTags(tagResult); - - AliRunTagCuts *runCuts = new AliRunTagCuts(); - AliLHCTagCuts *lhcCuts = new AliLHCTagCuts(); - AliDetectorTagCuts *detCuts = new AliDetectorTagCuts(); - AliEventTagCuts *evCuts = new AliEventTagCuts(); - - // Check if the cuts configuration file was provided - if (!gSystem->AccessPathName("ConfigureCuts.C")) { - gROOT->LoadMacro("ConfigureCuts.C"); - ConfigureCuts(runCuts, lhcCuts, detCuts, evCuts); - } - - TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts); - if (!chain || !chain->GetNtrees()) return NULL; - chain->ls(); - return chain; -} - -//______________________________________________________________________________ -TChain* CreateChainFromFile(const char *rootfile) -{ - // Create a chain using the root file. - TChain* chain = new TChain("esdTree"); - chain->Add(rootfile); - if (!chain->GetNtrees()) return NULL; - chain->ls(); - return chain; -} - -//______________________________________________________________________________ -TChain* CreateChainFromESDList(const char *esdList) -{ - // Create a chain using tags from the run list. - TChain* chain = new TChain("esdTree"); - ifstream inFile(esdList); - TString inFileName; - if (inFile.is_open()) { - while (! inFile.eof() ) { - inFileName.ReadLine(inFile,kFALSE); - if(!inFileName.EndsWith(".root")) continue; - chain->Add(inFileName.Data()); - } - } - inFile.close(); - if (!chain->GetNtrees()) return NULL; - chain->ls(); - return chain; -} - -//______________________________________________________________________________ -TChain* CreateChain(Int_t mode, TString input) -{ - printf("*******************************\n"); - printf("*** Getting the Chain ***\n"); - printf("*******************************\n"); - if(mode == kInteractif_xml) return CreateChainFromCollection(input.Data()); - else if (mode == kInteractif_ESDList) return CreateChainFromESDList(input.Data()); - else if (mode == kLocal) return CreateChainFromFile(input.Data()); - else return NULL; } -- 2.43.0