[u/mrichter/AliRoot.git] / PWG3 / muondep / MuonResolution.C

//--------------------------------------------------------------------------
// 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 <fstream>
#include <TString.h>
#include <TStopwatch.h>
#include <TMultiGraph.h>
#include <TSystem.h>
#include <TChain.h>
#include <TGraphErrors.h>
#include <TProof.h>
#include <TList.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TGrid.h>
#include <TEnv.h>
#include <TROOT.h>
#include "TAxis.h"
#include "THashList.h"
#include <TAlienCollection.h>
#include <TGridCollection.h>
#include <TGridResult.h>

// 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] "<<flush;
    while (remove != 'y' && remove != 'n') cin>>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 "<<iStep+1<<"/"<<nSteps<<endl;
    
    // Connect to proof if needed and prepare environment
    if (mode == kProof) LoadAlirootOnProof(alirootVersion, extraLibs, iStep);
    
    // create the analysis train
    AliAnalysisTaskMuonResolution *muonResolution = CreateAnalysisTrain(mode, iStep, selectPhysics, selectTrigger, matchTrig,
				      applyAccCut, minMomentum, correctForSystematics, extrapMode, clusterResNB, clusterResB);
    if (!muonResolution) return;
    
    // start analysis
    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
    if (mgr->InitAnalysis()) {
      mgr->PrintStatus();
      if (mode == kProof) mgr->StartAnalysis("proof", Form("%s",static_cast<TObjString*>(inputObj)->GetName()), nevents);
      else mgr->StartAnalysis("local", static_cast<TChain*>(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<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
  TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0;
  TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast<TGraphErrors*>(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<TObjArray*>(file->FindObjectAny("ChamberRes"));
  if (!summary) {
    Error("AddMCHViews","resolution graphs do not exist!");
    return;
  }
  
  TGraphErrors* g = 0x0;
  g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerDESigma"));
  if (g) {
    file->cd();
    ConvertGraph(*g, "resoX")->Write();
  }
  
  g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerDESigma"));
  if (g) {
    file->cd();
    ConvertGraph(*g, "resoY")->Write();
  }
  
  g = static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerDEMean_ClusterOut"));
  if (g) {
    file->cd();
    ConvertGraph(*g, "shiftX")->Write();
  }
  
  g = static_cast<TGraphErrors*>(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;
}
Commit	Line	Data
00612ffc	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
	4	//
	5	// Author: Philippe Pillot - SUBATECH Nantes
	6	//--------------------------------------------------------------------------
	7
	8	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	9	// ROOT includes
	10	#include <fstream>
	11	#include <TString.h>
	12	#include <TStopwatch.h>
	13	#include <TMultiGraph.h>
	14	#include <TSystem.h>
	15	#include <TChain.h>
	16	#include <TGraphErrors.h>
	17	#include <TProof.h>
	18	#include <TList.h>
	19	#include <TCanvas.h>
	20	#include <TFile.h>
	21	#include <TGrid.h>
	22	#include <TEnv.h>
	23	#include <TROOT.h>
	24	#include "TAxis.h"
	25	#include "THashList.h"
	26	#include <TAlienCollection.h>
	27	#include <TGridCollection.h>
	28	#include <TGridResult.h>
	29
	30	// STEER includes
	31	#include "AliLog.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"
	45
	46	// MUON includes
	47	#include "AliMpCDB.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"
	55
	56	#include "AddTaskPhysicsSelection.C"
	57	#include "AddTaskMuonResolution.C"
	58
	59	#endif
	60
	61	enum {kLocal, kInteractif_xml, kInteractif_ESDList, kProof};
	62
	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);
76
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)
80	{
81	/// Compute the cluster resolution by studying cluster-track residual, deconvoluting from track resolution
82
83	// timer start...
84	TStopwatch* localTimer = new TStopwatch;
85
86	// check parameters
87	nSteps = TMath::Max(nSteps,1);
88	if (extrapMode != 0 && extrapMode != 1) {
89	Error("MuonResolution","incorrect extrapolation mode!");
90	return;
91	}
92
93	// Check runing mode
94	Int_t mode = GetMode(smode, inputFileName);
95	if(mode < 0){
96	Error("MuonResolution","Please provide either an ESD root file, a list of ESDs, a collection of ESDs or a dataset.");
97	return;
98	}
99
100	// check for old output file to removed
101	char remove = '\0';
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");
106	else {
107	Error("MuonResolution","cannot proceed with these files there otherwise results will be mixed up!");
108	return;
109	}
110	}
111
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;
117
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.};
123
124	// output graphs
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");
135
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");
141	}
142
143	// loop over step
144	for (Int_t iStep = 0; iStep < nSteps; iStep++) {
145	cout<<"step "<<iStep+1<<"/"<<nSteps<<endl;
146
147	// Connect to proof if needed and prepare environment
148	if (mode == kProof) LoadAlirootOnProof(alirootVersion, extraLibs, iStep);
149
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;
154
155	// start analysis
156	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
157	if (mgr->InitAnalysis()) {
158	mgr->PrintStatus();
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);
161	}
162
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()) {
167	outFile->cd();
168	muonResolution->GetCanvases()->Write();
169	AddMCHViews(outFile);
170	outFile->Close();
171	}
172	}
173
174	// fill graph with starting resolutions from the task at first step
175	if (iStep == 0) {
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]);
182	}
183	}
184
185	// read the chamber resolution from the output file
186	if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr)) return;
187
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]);
194	}
195
196	// clean memory
197	delete mgr;
198	TObject::SetObjectStat(kFALSE);
199
200	}
201
202	// copy final results in results.root file
203	gSystem->Exec(Form("cp chamberResolution_step%d.root results.root", nSteps-1));
204
205	// display convergence
206	TCanvas* convergence = new TCanvas("convergence","convergence");
207	convergence->Divide(1,2);
208	convergence->cd(1);
209	mgClusterResXVsStep->Draw("ap");
210	convergence->cd(2);
211	mgClusterResYVsStep->Draw("ap");
212
213	// save convergence plots
214	TFile* outFile = TFile::Open("results.root","UPDATE");
215	if (!outFile \|\| !outFile->IsOpen()) return;
216	outFile->cd();
217	mgClusterResXVsStep->Write();
218	mgClusterResYVsStep->Write();
219	convergence->Write();
220	outFile->Close();
221
222	// ...timer stop
223	localTimer->Stop();
224	localTimer->Print();
225
226	}
227
228	//______________________________________________________________________________
229	void LoadAlirootOnProof(TString alirootVersion, TString& extraLibs, Int_t iStep)
230	{
231	/// Load aliroot packages and set environment on Proof
232
233	// set general environment and close previous session
234	if (iStep == 0) gEnv->SetValue("XSec.GSI.DelegProxy","2");
235	else gProof->Close("s");
236
237	// connect
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());
242	if (!gProof) return;
243
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);
251
252	// compile task on workers
253	if (!gSystem->AccessPathName("AliAnalysisTaskMuonResolution.cxx"))
254	gProof->Load("AliAnalysisTaskMuonResolution.cxx++g", kTRUE);
255
256	}
257
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])
262	{
263	/// create the analysis train and configure it
264
265	// Create the analysis manager
266	AliAnalysisManager *mgr = new AliAnalysisManager("MuonResolutionAnalysis");
267
268	// ESD input handler
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);
275
276	// event selection
277	if (selectPhysics) {
278	AliPhysicsSelectionTask* physicsSelection = AddTaskPhysicsSelection();
279	if (!physicsSelection) {
280	Error("run","AliPhysicsSelectionTask not created!");
281	return 0x0;
282	}
283	}
284
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!");
291	return 0x0;
292	}
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);
297
298	return muonResolution;
299
300	}
301
302	//______________________________________________________________________________
303	Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10])
304	{
305	/// read the chamber resolution from the output file
306
307	TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
308
309	if (!outFile \|\| !outFile->IsOpen()) {
310	Error("GetChamberResolution","output file does not exist!");
311	return kFALSE;
312	}
313
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;
317
318	if (!gCombinedResidualXPerChSigma \|\| !gCombinedResidualYPerChSigma) {
319	Error("GetChamberResolution","resolution graphs do not exist!");
320	return kFALSE;
321	}
322
323	Double_t dummy;
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);
329	}
330
331	outFile->Close();
332
333	return kTRUE;
334	}
335
336	//______________________________________________________________________________
337	void AddMCHViews(TFile* file)
338	{
339	/// Get from the file the graphs containing data per DE, convert them into mchview objects and save them
340
341	if ( ! AliMpDDLStore::Instance(false) )
342	{
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);
346	}
347
348	AliMpCDB::LoadAll();
349
350	TObjArray* summary = static_cast<TObjArray*>(file->FindObjectAny("ChamberRes"));
351	if (!summary) {
352	Error("AddMCHViews","resolution graphs do not exist!");
353	return;
354	}
355
356	TGraphErrors* g = 0x0;
357	g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerDESigma"));
358	if (g) {
359	file->cd();
360	ConvertGraph(*g, "resoX")->Write();
361	}
362
363	g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerDESigma"));
364	if (g) {
365	file->cd();
366	ConvertGraph(*g, "resoY")->Write();
367	}
368
369	g = static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerDEMean_ClusterOut"));
370	if (g) {
371	file->cd();
372	ConvertGraph(*g, "shiftX")->Write();
373	}
374
375	g = static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerDEMean_ClusterOut"));
376	if (g) {
377	file->cd();
378	ConvertGraph(*g, "shiftY")->Write();
379	}
380	}
381
382	//______________________________________________________________________________
383	AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name)
384	{
385	/// Convert graph containing data per DE into mchview object
386
387	AliMUON2DMap deValues(kFALSE);
388
389	for ( Int_t i = 0 ; i < g.GetN(); ++i )
390	{
391	double y = g.GetY()[i];
392	double ey = g.GetEY()[i];
393	int detElemId;
394	sscanf(g.GetXaxis()->GetBinLabel(i+1),"%d",&detElemId);
395
396	AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
397
398	AliMUONVCalibParam* param = new AliMUONCalibParamND(5, 1, detElemId, 0);
399
400	Double_t sumn = 1000.0;
401	Double_t sumw = sumn*y;
402	Double_t sumw2 = (sumn-1)eyey+sumw*sumw/sumn;
403
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);
409
410	deValues.Add(param);
411	}
412
413	AliMUONTrackerData* data = new AliMUONTrackerData(name,name,deValues,1);
414	data->SetDimensionName(0,name);
415
416	return data;
417	}
418
419	//______________________________________________________________________________
420	Int_t GetMode(TString smode, TString input)
421	{
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;
427	return -1;
428	}
429
430	//______________________________________________________________________________
431	TChain* CreateChainFromCollection(const char *xmlfile)
432	{
433	// Create a chain from the collection of tags.
434	if (!TGrid::Connect("alien://")) return NULL;
435
436	TGridCollection* coll = TAlienCollection::Open(xmlfile);
437	if (!coll) {
438	Error("CreateChainFromCollection", "Cannot create the AliEn collection");
439	return NULL;
440	}
441
442	TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE);
443	AliTagAnalysis *tagAna = new AliTagAnalysis("ESD");
444	tagAna->ChainGridTags(tagResult);
445
446	AliRunTagCuts *runCuts = new AliRunTagCuts();
447	AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();
448	AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
449	AliEventTagCuts *evCuts = new AliEventTagCuts();
450
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));
455
456	TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts);
457	if (!chain \|\| !chain->GetNtrees()) return NULL;
458	chain->ls();
459	return chain;
460	}
461
462	//______________________________________________________________________________
463	TChain* CreateChainFromFile(const char *rootfile)
464	{
465	// Create a chain using the root file.
466	TChain* chain = new TChain("esdTree");
467	chain->Add(rootfile);
468	if (!chain->GetNtrees()) return NULL;
469	chain->ls();
470	return chain;
471	}
472
473	//______________________________________________________________________________
474	TChain* CreateChainFromESDList(const char *esdList)
475	{
476	// Create a chain using tags from the run list.
477	TChain* chain = new TChain("esdTree");
478	ifstream inFile(esdList);
479	TString inFileName;
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());
485	}
486	}
487	inFile.close();
488	if (!chain->GetNtrees()) return NULL;
489	chain->ls();
490	return chain;
491	}
492
493	//______________________________________________________________________________
494	TChain* CreateChain(Int_t mode, TString input)
495	{
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());
502	else return NULL;
503	}
504