1 //--------------------------------------------------------------------------
2 // Macro compiled and launch by RunMuonResolution.C for submitting muon Resolution analysis locally or on CAF.
3 // See RunMuonResolution.C for more details
5 // Author: Philippe Pillot - SUBATECH Nantes
6 //--------------------------------------------------------------------------
8 #if !defined(__CINT__) || defined(__MAKECINT__)
12 #include <TStopwatch.h>
13 #include <TMultiGraph.h>
16 #include <TGraphErrors.h>
25 #include "THashList.h"
26 #include <TAlienCollection.h>
27 #include <TGridCollection.h>
28 #include <TGridResult.h>
32 #include "AliCDBManager.h"
33 #include "AliAnalysisManager.h"
34 #include "AliESDInputHandler.h"
35 #include "AliTagAnalysis.h"
36 #include "AliRunTagCuts.h"
37 #include "AliLHCTagCuts.h"
38 #include "AliDetectorTagCuts.h"
39 #include "AliEventTagCuts.h"
40 #include "AliPhysicsSelectionTask.h"
41 #include "AliPhysicsSelection.h"
42 #include "AliBackgroundSelection.h"
43 #include "AliAnalysisDataContainer.h"
44 #include "AliAnalysisTaskMuonResolution.h"
48 #include "AliMpDetElement.h"
49 #include "AliMpDDLStore.h"
50 #include "AliMUONCalibParamND.h"
51 #include "AliMUON2DMap.h"
52 #include "AliMUONTrackerData.h"
53 #include "AliMUONPainterDataRegistry.h"
54 #include "AliMUONTrackerDataWrapper.h"
56 #include "AddTaskPhysicsSelection.C"
57 #include "AddTaskMuonResolution.C"
61 enum {kLocal, kInteractif_xml, kInteractif_ESDList, kProof};
63 void LoadAlirootOnProof(TString alirootVersion, TString& extraLibs, Int_t iStep);
64 AliAnalysisTaskMuonResolution* CreateAnalysisTrain(Int_t mode, Int_t iStep, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig,
65 Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode,
66 Double_t clusterResNB[10], Double_t clusterResB[10]);
67 Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10],
68 Double_t clusterResNBErr[10], Double_t clusterResBErr[10]);
69 void AddMCHViews(TFile* file);
70 AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name);
71 Int_t GetMode(TString smode, TString input);
72 TChain* CreateChainFromCollection(const char *xmlfile);
73 TChain* CreateChainFromFile(const char *rootfile);
74 TChain* CreateChainFromESDList(const char *esdList);
75 TChain* CreateChain(Int_t mode, TString input);
77 //______________________________________________________________________________
78 void MuonResolution(TString smode, TString inputFileName, TString alirootVersion, Int_t nSteps, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig,
79 Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode, Int_t nevents, TString extraLibs)
81 /// Compute the cluster resolution by studying cluster-track residual, deconvoluting from track resolution
84 TStopwatch* localTimer = new TStopwatch;
87 nSteps = TMath::Max(nSteps,1);
88 if (extrapMode != 0 && extrapMode != 1) {
89 Error("MuonResolution","incorrect extrapolation mode!");
94 Int_t mode = GetMode(smode, inputFileName);
96 Error("MuonResolution","Please provide either an ESD root file, a list of ESDs, a collection of ESDs or a dataset.");
100 // check for old output file to removed
102 if (!gSystem->Exec("ls chamberResolution_step*[0-9].root")) {
103 cout<<"above files must be removed from the current directory. Delete? [y=yes, n=no] "<<flush;
104 while (remove != 'y' && remove != 'n') cin>>remove;
105 if (remove == 'y') gSystem->Exec("rm -f chamberResolution_step*[0-9].root");
107 Error("MuonResolution","cannot proceed with these files there otherwise results will be mixed up!");
112 // Create input object
113 TObject* inputObj = 0x0;
114 if (mode == kProof) inputObj = new TObjString(inputFileName);
115 else inputObj = CreateChain(mode, inputFileName);
116 if (!inputObj) return;
118 // set starting chamber resolution (if -1 they will be loaded from recoParam in the task)
119 Double_t clusterResNB[10] ={-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.};
120 Double_t clusterResB[10] ={-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.};
121 Double_t clusterResNBErr[10] ={0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
122 Double_t clusterResBErr[10] ={0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
125 TMultiGraph* mgClusterResXVsStep = new TMultiGraph("mgClusterResXVsStep","cluster X-resolution versus step;step;#sigma_{X} (cm)");
126 TMultiGraph* mgClusterResYVsStep = new TMultiGraph("mgClusterResYVsStep","cluster Y-resolution versus step;step;#sigma_{Y} (cm)");
127 TGraphErrors* clusterResXVsStep[10];
128 TGraphErrors* clusterResYVsStep[10];
129 for (Int_t i = 0; i < 10; i++) {
130 clusterResXVsStep[i] = new TGraphErrors(nSteps+1);
131 clusterResXVsStep[i]->SetName(Form("gResX_ch%d",i+1));
132 clusterResXVsStep[i]->SetMarkerStyle(kFullDotMedium);
133 clusterResXVsStep[i]->SetMarkerColor(i+1+i/9);
134 mgClusterResXVsStep->Add(clusterResXVsStep[i],"lp");
136 clusterResYVsStep[i] = new TGraphErrors(nSteps+1);
137 clusterResYVsStep[i]->SetName(Form("gResY_ch%d",i+1));
138 clusterResYVsStep[i]->SetMarkerStyle(kFullDotMedium);
139 clusterResYVsStep[i]->SetMarkerColor(i+1+i/9);
140 mgClusterResYVsStep->Add(clusterResYVsStep[i],"lp");
144 for (Int_t iStep = 0; iStep < nSteps; iStep++) {
145 cout<<"step "<<iStep+1<<"/"<<nSteps<<endl;
147 // Connect to proof if needed and prepare environment
148 if (mode == kProof) LoadAlirootOnProof(alirootVersion, extraLibs, iStep);
150 // create the analysis train
151 AliAnalysisTaskMuonResolution *muonResolution = CreateAnalysisTrain(mode, iStep, selectPhysics, selectTrigger, matchTrig,
152 applyAccCut, minMomentum, correctForSystematics, extrapMode, clusterResNB, clusterResB);
153 if (!muonResolution) return;
156 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
157 if (mgr->InitAnalysis()) {
159 if (mode == kProof) mgr->StartAnalysis("proof", Form("%s",static_cast<TObjString*>(inputObj)->GetName()), nevents);
160 else mgr->StartAnalysis("local", static_cast<TChain*>(inputObj), nevents);
163 // save the summary canvases and mchview display
164 if (muonResolution->GetCanvases()) {
165 TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"UPDATE");
166 if (outFile && outFile->IsOpen()) {
168 muonResolution->GetCanvases()->Write();
169 AddMCHViews(outFile);
174 // fill graph with starting resolutions from the task at first step
176 muonResolution->GetStartingResolution(clusterResNB, clusterResB);
177 for (Int_t i = 0; i < 10; i++) {
178 clusterResXVsStep[i]->SetPoint(0, 0, clusterResNB[i]);
179 clusterResXVsStep[i]->SetPointError(0, 0., clusterResNBErr[i]);
180 clusterResYVsStep[i]->SetPoint(0, 0, clusterResB[i]);
181 clusterResYVsStep[i]->SetPointError(0, 0., clusterResBErr[i]);
185 // read the chamber resolution from the output file
186 if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr)) return;
188 // fill graphs with computed resolutions
189 for (Int_t i = 0; i < 10; i++) {
190 clusterResXVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResNB[i]);
191 clusterResXVsStep[i]->SetPointError(iStep+1, 0., clusterResNBErr[i]);
192 clusterResYVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResB[i]);
193 clusterResYVsStep[i]->SetPointError(iStep+1, 0., clusterResBErr[i]);
198 TObject::SetObjectStat(kFALSE);
202 // copy final results in results.root file
203 gSystem->Exec(Form("cp chamberResolution_step%d.root results.root", nSteps-1));
205 // display convergence
206 TCanvas* convergence = new TCanvas("convergence","convergence");
207 convergence->Divide(1,2);
209 mgClusterResXVsStep->Draw("ap");
211 mgClusterResYVsStep->Draw("ap");
213 // save convergence plots
214 TFile* outFile = TFile::Open("results.root","UPDATE");
215 if (!outFile || !outFile->IsOpen()) return;
217 mgClusterResXVsStep->Write();
218 mgClusterResYVsStep->Write();
219 convergence->Write();
228 //______________________________________________________________________________
229 void LoadAlirootOnProof(TString alirootVersion, TString& extraLibs, Int_t iStep)
231 /// Load aliroot packages and set environment on Proof
233 // set general environment and close previous session
234 if (iStep == 0) gEnv->SetValue("XSec.GSI.DelegProxy","2");
235 else gProof->Close("s");
238 TString location = "alice-caf.cern.ch";
239 TString nWorkers = "";
240 if (gSystem->Getenv("alien_API_USER") == NULL) TProof::Open(location.Data(), nWorkers.Data());
241 else TProof::Open(Form("%s@%s",gSystem->Getenv("alien_API_USER"), location.Data()), nWorkers.Data());
244 // set environment and load libraries on workers
245 TList* list = new TList();
246 list->Add(new TNamed("ALIROOT_MODE", ""));
247 list->Add(new TNamed("ALIROOT_EXTRA_LIBS", extraLibs.Data()));
248 if (!gSystem->AccessPathName("AliAnalysisTaskMuonResolution.cxx"))
249 list->Add(new TNamed("ALIROOT_EXTRA_INCLUDES", "MUON:MUON/mapping"));
250 gProof->EnablePackage(alirootVersion.Data(), list, kTRUE);
252 // compile task on workers
253 if (!gSystem->AccessPathName("AliAnalysisTaskMuonResolution.cxx"))
254 gProof->Load("AliAnalysisTaskMuonResolution.cxx++g", kTRUE);
258 //______________________________________________________________________________
259 AliAnalysisTaskMuonResolution* CreateAnalysisTrain(Int_t mode, Int_t iStep, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig,
260 Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode,
261 Double_t clusterResNB[10], Double_t clusterResB[10])
263 /// create the analysis train and configure it
265 // Create the analysis manager
266 AliAnalysisManager *mgr = new AliAnalysisManager("MuonResolutionAnalysis");
269 AliESDInputHandler* esdH = new AliESDInputHandler();
270 esdH->SetReadFriends(kFALSE);
271 esdH->SetInactiveBranches(" FMD PHOS EMCAL Pmd Trd V0s TPC "
272 "Cascades Kinks CaloClusters ACORDE RawData HLT TZERO ZDC"
273 " Cells ACORDE Pileup");
274 mgr->SetInputEventHandler(esdH);
278 AliPhysicsSelectionTask* physicsSelection = AddTaskPhysicsSelection();
279 if (!physicsSelection) {
280 Error("run","AliPhysicsSelectionTask not created!");
285 // Muon Resolution analysis
286 TString outputFileName = Form("chamberResolution_step%d.root", iStep);
287 AliAnalysisManager::SetCommonFileName(outputFileName.Data());
288 AliAnalysisTaskMuonResolution *muonResolution = AddTaskMuonResolution(selectPhysics, selectTrigger, matchTrig, applyAccCut, minMomentum, correctForSystematics, extrapMode);
289 if (!muonResolution) {
290 Error("run","AliAnalysisTaskMuonResolution not created!");
293 if (mode == kLocal) muonResolution->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
294 if (mode != kProof) muonResolution->ShowProgressBar();
295 muonResolution->PrintClusterRes(kTRUE, kTRUE);
296 muonResolution->SetStartingResolution(clusterResNB, clusterResB);
298 return muonResolution;
302 //______________________________________________________________________________
303 Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10])
305 /// read the chamber resolution from the output file
307 TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
309 if (!outFile || !outFile->IsOpen()) {
310 Error("GetChamberResolution","output file does not exist!");
314 TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
315 TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0;
316 TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerChSigma")) : 0x0;
318 if (!gCombinedResidualXPerChSigma || !gCombinedResidualYPerChSigma) {
319 Error("GetChamberResolution","resolution graphs do not exist!");
324 for (Int_t i = 0; i < 10; i++) {
325 gCombinedResidualXPerChSigma->GetPoint(i, dummy, clusterResNB[i]);
326 gCombinedResidualYPerChSigma->GetPoint(i, dummy, clusterResB[i]);
327 clusterResNBErr[i] = gCombinedResidualXPerChSigma->GetErrorY(i);
328 clusterResBErr[i] = gCombinedResidualYPerChSigma->GetErrorY(i);
336 //______________________________________________________________________________
337 void AddMCHViews(TFile* file)
339 /// Get from the file the graphs containing data per DE, convert them into mchview objects and save them
341 if ( ! AliMpDDLStore::Instance(false) )
343 Warning("AddMCHViews","mapping was not loaded. Loading it from $ALICE_ROOT/OCDB");
344 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
345 AliCDBManager::Instance()->SetRun(999999999);
350 TObjArray* summary = static_cast<TObjArray*>(file->FindObjectAny("ChamberRes"));
352 Error("AddMCHViews","resolution graphs do not exist!");
356 TGraphErrors* g = 0x0;
357 g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerDESigma"));
360 ConvertGraph(*g, "resoX")->Write();
363 g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerDESigma"));
366 ConvertGraph(*g, "resoY")->Write();
369 g = static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerDEMean_ClusterOut"));
372 ConvertGraph(*g, "shiftX")->Write();
375 g = static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerDEMean_ClusterOut"));
378 ConvertGraph(*g, "shiftY")->Write();
382 //______________________________________________________________________________
383 AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name)
385 /// Convert graph containing data per DE into mchview object
387 AliMUON2DMap deValues(kFALSE);
389 for ( Int_t i = 0 ; i < g.GetN(); ++i )
391 double y = g.GetY()[i];
392 double ey = g.GetEY()[i];
394 sscanf(g.GetXaxis()->GetBinLabel(i+1),"%d",&detElemId);
396 AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
398 AliMUONVCalibParam* param = new AliMUONCalibParamND(5, 1, detElemId, 0);
400 Double_t sumn = 1000.0;
401 Double_t sumw = sumn*y;
402 Double_t sumw2 = (sumn-1)*ey*ey+sumw*sumw/sumn;
404 param->SetValueAsDouble(0,0,sumw);
405 param->SetValueAsDouble(0,1,sumw2);
406 param->SetValueAsDouble(0,2,sumn);
407 param->SetValueAsDouble(0,3,de->NofChannels());
408 param->SetValueAsDouble(0,4,1);
413 AliMUONTrackerData* data = new AliMUONTrackerData(name,name,deValues,1);
414 data->SetDimensionName(0,name);
419 //______________________________________________________________________________
420 Int_t GetMode(TString smode, TString input)
422 if (smode == "local") {
423 if ( input.EndsWith(".xml") ) return kInteractif_xml;
424 else if ( input.EndsWith(".txt") ) return kInteractif_ESDList;
425 else if ( input.EndsWith(".root") ) return kLocal;
426 } else if (smode == "proof") return kProof;
430 //______________________________________________________________________________
431 TChain* CreateChainFromCollection(const char *xmlfile)
433 // Create a chain from the collection of tags.
434 if (!TGrid::Connect("alien://")) return NULL;
436 TGridCollection* coll = TAlienCollection::Open(xmlfile);
438 Error("CreateChainFromCollection", "Cannot create the AliEn collection");
442 TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE);
443 AliTagAnalysis *tagAna = new AliTagAnalysis("ESD");
444 tagAna->ChainGridTags(tagResult);
446 AliRunTagCuts *runCuts = new AliRunTagCuts();
447 AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();
448 AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
449 AliEventTagCuts *evCuts = new AliEventTagCuts();
451 // Check if the cuts configuration file was provided
452 if (!gSystem->AccessPathName("ConfigureCuts.C"))
453 gROOT->ProcessLine(Form(".x ConfigureCuts.C((AliRunTagCuts*)%p, (AliLHCTagCuts*)%p, (AliDetectorTagCuts*)%p,"
454 " (AliEventTagCuts*)%p)", runCuts, lhcCuts, detCuts, evCuts));
456 TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts);
457 if (!chain || !chain->GetNtrees()) return NULL;
462 //______________________________________________________________________________
463 TChain* CreateChainFromFile(const char *rootfile)
465 // Create a chain using the root file.
466 TChain* chain = new TChain("esdTree");
467 chain->Add(rootfile);
468 if (!chain->GetNtrees()) return NULL;
473 //______________________________________________________________________________
474 TChain* CreateChainFromESDList(const char *esdList)
476 // Create a chain using tags from the run list.
477 TChain* chain = new TChain("esdTree");
478 ifstream inFile(esdList);
480 if (inFile.is_open()) {
481 while (! inFile.eof() ) {
482 inFileName.ReadLine(inFile,kFALSE);
483 if(!inFileName.EndsWith(".root")) continue;
484 chain->Add(inFileName.Data());
488 if (!chain->GetNtrees()) return NULL;
493 //______________________________________________________________________________
494 TChain* CreateChain(Int_t mode, TString input)
496 printf("*******************************\n");
497 printf("*** Getting the Chain ***\n");
498 printf("*******************************\n");
499 if(mode == kInteractif_xml) return CreateChainFromCollection(input.Data());
500 else if (mode == kInteractif_ESDList) return CreateChainFromESDList(input.Data());
501 else if (mode == kLocal) return CreateChainFromFile(input.Data());