[u/mrichter/AliRoot.git] / PWGPP / MUON / dep / 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 "AliCentralitySelectionTask.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 "AddTaskCentrality.C"
#include "AddTaskMuonResolution.C"

#endif

enum {kLocal, kInteractif_xml, kInteractif_ESDList, kProof};

void    LoadAlirootOnProof(TString& aaf, 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(smode, 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& aaf, 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 = (aaf == "caf") ? "alice-caf.cern.ch" : "nansafmaster.in2p3.fr";
  //TString location = (aaf == "caf") ? "alice-caf.cern.ch" : "localhost:1093";
  TString nWorkers = (aaf == "caf") ? "workers=80" : "";
  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("*");
  esdH->SetActiveBranches("MuonTracks AliESDRun. AliESDHeader. AliMultiplicity. AliESDFMD. AliESDVZERO. SPDVertex. PrimaryVertex. AliESDZDC.");
  mgr->SetInputEventHandler(esdH);
  
  // event selection
  if (selectPhysics) {
    AliPhysicsSelectionTask* physicsSelection = AddTaskPhysicsSelection();
    if (!physicsSelection) {
      Error("CreateAnalysisTrain","AliPhysicsSelectionTask not created!");
      return 0x0;
    }
  }
  
  // centrality selection
  AliCentralitySelectionTask* centralityTask = AddTaskCentrality();
  if (!centralityTask) {
    Error("CreateAnalysisTrain","AliCentralitySelectionTask not created!");
    return 0x0;
  }
  centralityTask->SetPass(1);
  
  // 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("CreateAnalysisTrain","AliAnalysisTaskMuonResolution not created!");
    return 0x0;
  }
  //if (mode == kLocal) muonResolution->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  muonResolution->SetDefaultStorage("alien://folder=/alice/data/2011/OCDB");
  if (mode != kProof) muonResolution->ShowProgressBar();
  muonResolution->PrintClusterRes(kTRUE, kTRUE);
  muonResolution->SetStartingResolution(clusterResNB, clusterResB);
  //muonResolution->RemoveMonoCathodClusters(kTRUE, kFALSE);
//  muonResolution->FitResiduals(kFALSE);
//  muonResolution->ReAlign("alien://folder=/alice/cern.ch/user/p/ppillot/OCDB2011","");
//  muonResolution->ReAlign("alien://folder=/alice/cern.ch/user/p/ppillot/OCDB2011", "alien://folder=/alice/cern.ch/user/p/ppillot/OCDB2011_Align2");
  
  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 == "caf" || smode == "saf") 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"
4060971a	43	#include "AliCentralitySelectionTask.h"
00612ffc	44	#include "AliAnalysisDataContainer.h"
	45	#include "AliAnalysisTaskMuonResolution.h"
	46
	47	// MUON includes
	48	#include "AliMpCDB.h"
	49	#include "AliMpDetElement.h"
	50	#include "AliMpDDLStore.h"
	51	#include "AliMUONCalibParamND.h"
	52	#include "AliMUON2DMap.h"
	53	#include "AliMUONTrackerData.h"
	54	#include "AliMUONPainterDataRegistry.h"
	55	#include "AliMUONTrackerDataWrapper.h"
	56
	57	#include "AddTaskPhysicsSelection.C"
4060971a	58	#include "AddTaskCentrality.C"
00612ffc	59	#include "AddTaskMuonResolution.C"
	60
	61	#endif
	62
	63	enum {kLocal, kInteractif_xml, kInteractif_ESDList, kProof};
	64
4060971a	65	void LoadAlirootOnProof(TString& aaf, TString alirootVersion, TString& extraLibs, Int_t iStep);
00612ffc	66	AliAnalysisTaskMuonResolution* CreateAnalysisTrain(Int_t mode, Int_t iStep, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig,
	67	Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode,
	68	Double_t clusterResNB[10], Double_t clusterResB[10]);
	69	Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10],
	70	Double_t clusterResNBErr[10], Double_t clusterResBErr[10]);
	71	void AddMCHViews(TFile* file);
	72	AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name);
	73	Int_t GetMode(TString smode, TString input);
	74	TChain* CreateChainFromCollection(const char *xmlfile);
	75	TChain* CreateChainFromFile(const char *rootfile);
	76	TChain* CreateChainFromESDList(const char *esdList);
	77	TChain* CreateChain(Int_t mode, TString input);
	78
	79	//______________________________________________________________________________
	80	void MuonResolution(TString smode, TString inputFileName, TString alirootVersion, Int_t nSteps, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig,
	81	Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode, Int_t nevents, TString extraLibs)
	82	{
	83	/// Compute the cluster resolution by studying cluster-track residual, deconvoluting from track resolution
	84
	85	// timer start...
	86	TStopwatch* localTimer = new TStopwatch;
	87
	88	// check parameters
	89	nSteps = TMath::Max(nSteps,1);
	90	if (extrapMode != 0 && extrapMode != 1) {
	91	Error("MuonResolution","incorrect extrapolation mode!");
	92	return;
	93	}
	94
	95	// Check runing mode
	96	Int_t mode = GetMode(smode, inputFileName);
	97	if(mode < 0){
	98	Error("MuonResolution","Please provide either an ESD root file, a list of ESDs, a collection of ESDs or a dataset.");
	99	return;
	100	}
	101
	102	// check for old output file to removed
	103	char remove = '\0';
	104	if (!gSystem->Exec("ls chamberResolution_step*[0-9].root")) {
	105	cout<<"above files must be removed from the current directory. Delete? [y=yes, n=no] "<<flush;
	106	while (remove != 'y' && remove != 'n') cin>>remove;
	107	if (remove == 'y') gSystem->Exec("rm -f chamberResolution_step*[0-9].root");
	108	else {
	109	Error("MuonResolution","cannot proceed with these files there otherwise results will be mixed up!");
	110	return;
	111	}
	112	}
	113
	114	// Create input object
	115	TObject* inputObj = 0x0;
	116	if (mode == kProof) inputObj = new TObjString(inputFileName);
	117	else inputObj = CreateChain(mode, inputFileName);
	118	if (!inputObj) return;
	119
	120	// set starting chamber resolution (if -1 they will be loaded from recoParam in the task)
	121	Double_t clusterResNB[10] ={-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.};
	122	Double_t clusterResB[10] ={-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.};
	123	Double_t clusterResNBErr[10] ={0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
	124	Double_t clusterResBErr[10] ={0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
	125
	126	// output graphs
	127	TMultiGraph* mgClusterResXVsStep = new TMultiGraph("mgClusterResXVsStep","cluster X-resolution versus step;step;#sigma_{X} (cm)");
	128	TMultiGraph* mgClusterResYVsStep = new TMultiGraph("mgClusterResYVsStep","cluster Y-resolution versus step;step;#sigma_{Y} (cm)");
	129	TGraphErrors* clusterResXVsStep[10];
130	TGraphErrors* clusterResYVsStep[10];
131	for (Int_t i = 0; i < 10; i++) {
132	clusterResXVsStep[i] = new TGraphErrors(nSteps+1);
133	clusterResXVsStep[i]->SetName(Form("gResX_ch%d",i+1));
134	clusterResXVsStep[i]->SetMarkerStyle(kFullDotMedium);
135	clusterResXVsStep[i]->SetMarkerColor(i+1+i/9);
136	mgClusterResXVsStep->Add(clusterResXVsStep[i],"lp");
137
138	clusterResYVsStep[i] = new TGraphErrors(nSteps+1);
139	clusterResYVsStep[i]->SetName(Form("gResY_ch%d",i+1));
140	clusterResYVsStep[i]->SetMarkerStyle(kFullDotMedium);
141	clusterResYVsStep[i]->SetMarkerColor(i+1+i/9);
142	mgClusterResYVsStep->Add(clusterResYVsStep[i],"lp");
143	}
144
145	// loop over step
146	for (Int_t iStep = 0; iStep < nSteps; iStep++) {
147	cout<<"step "<<iStep+1<<"/"<<nSteps<<endl;
148
149	// Connect to proof if needed and prepare environment
4060971a	150	if (mode == kProof) LoadAlirootOnProof(smode, alirootVersion, extraLibs, iStep);
00612ffc	151
	152	// create the analysis train
	153	AliAnalysisTaskMuonResolution *muonResolution = CreateAnalysisTrain(mode, iStep, selectPhysics, selectTrigger, matchTrig,
	154	applyAccCut, minMomentum, correctForSystematics, extrapMode, clusterResNB, clusterResB);
	155	if (!muonResolution) return;
	156
	157	// start analysis
	158	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	159	if (mgr->InitAnalysis()) {
	160	mgr->PrintStatus();
	161	if (mode == kProof) mgr->StartAnalysis("proof", Form("%s",static_cast<TObjString*>(inputObj)->GetName()), nevents);
	162	else mgr->StartAnalysis("local", static_cast<TChain*>(inputObj), nevents);
	163	}
	164
	165	// save the summary canvases and mchview display
	166	if (muonResolution->GetCanvases()) {
	167	TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"UPDATE");
	168	if (outFile && outFile->IsOpen()) {
	169	outFile->cd();
	170	muonResolution->GetCanvases()->Write();
	171	AddMCHViews(outFile);
	172	outFile->Close();
	173	}
	174	}
	175
	176	// fill graph with starting resolutions from the task at first step
	177	if (iStep == 0) {
	178	muonResolution->GetStartingResolution(clusterResNB, clusterResB);
	179	for (Int_t i = 0; i < 10; i++) {
	180	clusterResXVsStep[i]->SetPoint(0, 0, clusterResNB[i]);
	181	clusterResXVsStep[i]->SetPointError(0, 0., clusterResNBErr[i]);
	182	clusterResYVsStep[i]->SetPoint(0, 0, clusterResB[i]);
	183	clusterResYVsStep[i]->SetPointError(0, 0., clusterResBErr[i]);
	184	}
	185	}
	186
	187	// read the chamber resolution from the output file
	188	if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr)) return;
	189
	190	// fill graphs with computed resolutions
	191	for (Int_t i = 0; i < 10; i++) {
	192	clusterResXVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResNB[i]);
	193	clusterResXVsStep[i]->SetPointError(iStep+1, 0., clusterResNBErr[i]);
	194	clusterResYVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResB[i]);
	195	clusterResYVsStep[i]->SetPointError(iStep+1, 0., clusterResBErr[i]);
	196	}
	197
	198	// clean memory
	199	delete mgr;
	200	TObject::SetObjectStat(kFALSE);
	201
	202	}
	203
	204	// copy final results in results.root file
	205	gSystem->Exec(Form("cp chamberResolution_step%d.root results.root", nSteps-1));
	206
	207	// display convergence
	208	TCanvas* convergence = new TCanvas("convergence","convergence");
	209	convergence->Divide(1,2);
	210	convergence->cd(1);
	211	mgClusterResXVsStep->Draw("ap");
	212	convergence->cd(2);
	213	mgClusterResYVsStep->Draw("ap");
	214
215	// save convergence plots
216	TFile* outFile = TFile::Open("results.root","UPDATE");
217	if (!outFile \|\| !outFile->IsOpen()) return;
218	outFile->cd();
219	mgClusterResXVsStep->Write();
220	mgClusterResYVsStep->Write();
221	convergence->Write();
222	outFile->Close();
223
224	// ...timer stop
225	localTimer->Stop();
226	localTimer->Print();
227
228	}
229
230	//______________________________________________________________________________
4060971a	231	void LoadAlirootOnProof(TString& aaf, TString alirootVersion, TString& extraLibs, Int_t iStep)
00612ffc	232	{
	233	/// Load aliroot packages and set environment on Proof
	234
	235	// set general environment and close previous session
	236	if (iStep == 0) gEnv->SetValue("XSec.GSI.DelegProxy","2");
	237	else gProof->Close("s");
	238
	239	// connect
4060971a	240	TString location = (aaf == "caf") ? "alice-caf.cern.ch" : "nansafmaster.in2p3.fr";
	241	//TString location = (aaf == "caf") ? "alice-caf.cern.ch" : "localhost:1093";
	242	TString nWorkers = (aaf == "caf") ? "workers=80" : "";
00612ffc	243	if (gSystem->Getenv("alien_API_USER") == NULL) TProof::Open(location.Data(), nWorkers.Data());
	244	else TProof::Open(Form("%s@%s",gSystem->Getenv("alien_API_USER"), location.Data()), nWorkers.Data());
	245	if (!gProof) return;
	246
	247	// set environment and load libraries on workers
	248	TList* list = new TList();
	249	list->Add(new TNamed("ALIROOT_MODE", ""));
	250	list->Add(new TNamed("ALIROOT_EXTRA_LIBS", extraLibs.Data()));
	251	if (!gSystem->AccessPathName("AliAnalysisTaskMuonResolution.cxx"))
	252	list->Add(new TNamed("ALIROOT_EXTRA_INCLUDES", "MUON:MUON/mapping"));
	253	gProof->EnablePackage(alirootVersion.Data(), list, kTRUE);
	254
	255	// compile task on workers
	256	if (!gSystem->AccessPathName("AliAnalysisTaskMuonResolution.cxx"))
	257	gProof->Load("AliAnalysisTaskMuonResolution.cxx++g", kTRUE);
	258
	259	}
	260
	261	//______________________________________________________________________________
	262	AliAnalysisTaskMuonResolution* CreateAnalysisTrain(Int_t mode, Int_t iStep, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig,
	263	Bool_t applyAccCut, Double_t minMomentum, Bool_t correctForSystematics, Int_t extrapMode,
	264	Double_t clusterResNB[10], Double_t clusterResB[10])
	265	{
	266	/// create the analysis train and configure it
	267
	268	// Create the analysis manager
	269	AliAnalysisManager *mgr = new AliAnalysisManager("MuonResolutionAnalysis");
	270
	271	// ESD input handler
	272	AliESDInputHandler* esdH = new AliESDInputHandler();
	273	esdH->SetReadFriends(kFALSE);
4060971a	274	esdH->SetInactiveBranches("*");
4060971a	275	esdH->SetActiveBranches("MuonTracks AliESDRun. AliESDHeader. AliMultiplicity. AliESDFMD. AliESDVZERO. SPDVertex. PrimaryVertex. AliESDZDC.");
00612ffc	276	mgr->SetInputEventHandler(esdH);
	277
	278	// event selection
	279	if (selectPhysics) {
	280	AliPhysicsSelectionTask* physicsSelection = AddTaskPhysicsSelection();
	281	if (!physicsSelection) {
4060971a	282	Error("CreateAnalysisTrain","AliPhysicsSelectionTask not created!");
00612ffc	283	return 0x0;
	284	}
	285	}
	286
4060971a	287	// centrality selection
	288	AliCentralitySelectionTask* centralityTask = AddTaskCentrality();
	289	if (!centralityTask) {
	290	Error("CreateAnalysisTrain","AliCentralitySelectionTask not created!");
	291	return 0x0;
	292	}
	293	centralityTask->SetPass(1);
	294
00612ffc	295	// Muon Resolution analysis
	296	TString outputFileName = Form("chamberResolution_step%d.root", iStep);
	297	AliAnalysisManager::SetCommonFileName(outputFileName.Data());
	298	AliAnalysisTaskMuonResolution *muonResolution = AddTaskMuonResolution(selectPhysics, selectTrigger, matchTrig, applyAccCut, minMomentum, correctForSystematics, extrapMode);
	299	if (!muonResolution) {
4060971a	300	Error("CreateAnalysisTrain","AliAnalysisTaskMuonResolution not created!");
00612ffc	301	return 0x0;
00612ffc	302	}
4060971a	303	//if (mode == kLocal) muonResolution->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
4060971a	304	muonResolution->SetDefaultStorage("alien://folder=/alice/data/2011/OCDB");
00612ffc	305	if (mode != kProof) muonResolution->ShowProgressBar();
	306	muonResolution->PrintClusterRes(kTRUE, kTRUE);
	307	muonResolution->SetStartingResolution(clusterResNB, clusterResB);
4060971a	308	//muonResolution->RemoveMonoCathodClusters(kTRUE, kFALSE);
	309	// muonResolution->FitResiduals(kFALSE);
	310	// muonResolution->ReAlign("alien://folder=/alice/cern.ch/user/p/ppillot/OCDB2011","");
	311	// muonResolution->ReAlign("alien://folder=/alice/cern.ch/user/p/ppillot/OCDB2011", "alien://folder=/alice/cern.ch/user/p/ppillot/OCDB2011_Align2");
00612ffc	312
	313	return muonResolution;
	314
	315	}
	316
	317	//______________________________________________________________________________
	318	Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10])
	319	{
	320	/// read the chamber resolution from the output file
	321
	322	TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
	323
	324	if (!outFile \|\| !outFile->IsOpen()) {
	325	Error("GetChamberResolution","output file does not exist!");
	326	return kFALSE;
	327	}
	328
	329	TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
	330	TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0;
	331	TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerChSigma")) : 0x0;
	332
	333	if (!gCombinedResidualXPerChSigma \|\| !gCombinedResidualYPerChSigma) {
	334	Error("GetChamberResolution","resolution graphs do not exist!");
	335	return kFALSE;
	336	}
	337
	338	Double_t dummy;
	339	for (Int_t i = 0; i < 10; i++) {
	340	gCombinedResidualXPerChSigma->GetPoint(i, dummy, clusterResNB[i]);
	341	gCombinedResidualYPerChSigma->GetPoint(i, dummy, clusterResB[i]);
	342	clusterResNBErr[i] = gCombinedResidualXPerChSigma->GetErrorY(i);
	343	clusterResBErr[i] = gCombinedResidualYPerChSigma->GetErrorY(i);
	344	}
	345
	346	outFile->Close();
	347
	348	return kTRUE;
	349	}
	350
	351	//______________________________________________________________________________
	352	void AddMCHViews(TFile* file)
	353	{
	354	/// Get from the file the graphs containing data per DE, convert them into mchview objects and save them
	355
	356	if ( ! AliMpDDLStore::Instance(false) )
	357	{
	358	Warning("AddMCHViews","mapping was not loaded. Loading it from $ALICE_ROOT/OCDB");
	359	AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
	360	AliCDBManager::Instance()->SetRun(999999999);
	361	}
	362
	363	AliMpCDB::LoadAll();
	364
	365	TObjArray* summary = static_cast<TObjArray*>(file->FindObjectAny("ChamberRes"));
	366	if (!summary) {
	367	Error("AddMCHViews","resolution graphs do not exist!");
	368	return;
	369	}
	370
	371	TGraphErrors* g = 0x0;
	372	g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerDESigma"));
	373	if (g) {
	374	file->cd();
	375	ConvertGraph(*g, "resoX")->Write();
376	}
377
378	g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerDESigma"));
379	if (g) {
380	file->cd();
381	ConvertGraph(*g, "resoY")->Write();
382	}
383
384	g = static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerDEMean_ClusterOut"));
385	if (g) {
386	file->cd();
387	ConvertGraph(*g, "shiftX")->Write();
388	}
389
390	g = static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerDEMean_ClusterOut"));
391	if (g) {
392	file->cd();
393	ConvertGraph(*g, "shiftY")->Write();
394	}
395	}
396
397	//______________________________________________________________________________
398	AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name)
399	{
400	/// Convert graph containing data per DE into mchview object
401
402	AliMUON2DMap deValues(kFALSE);
403
404	for ( Int_t i = 0 ; i < g.GetN(); ++i )
405	{
406	double y = g.GetY()[i];
407	double ey = g.GetEY()[i];
408	int detElemId;
409	sscanf(g.GetXaxis()->GetBinLabel(i+1),"%d",&detElemId);
410
411	AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
412
413	AliMUONVCalibParam* param = new AliMUONCalibParamND(5, 1, detElemId, 0);
414
415	Double_t sumn = 1000.0;
416	Double_t sumw = sumn*y;
417	Double_t sumw2 = (sumn-1)eyey+sumw*sumw/sumn;
418
419	param->SetValueAsDouble(0,0,sumw);
420	param->SetValueAsDouble(0,1,sumw2);
421	param->SetValueAsDouble(0,2,sumn);
422	param->SetValueAsDouble(0,3,de->NofChannels());
423	param->SetValueAsDouble(0,4,1);
424
425	deValues.Add(param);
426	}
427
428	AliMUONTrackerData* data = new AliMUONTrackerData(name,name,deValues,1);
429	data->SetDimensionName(0,name);
430
431	return data;
432	}
433
434	//______________________________________________________________________________
435	Int_t GetMode(TString smode, TString input)
436	{
437	if (smode == "local") {
438	if ( input.EndsWith(".xml") ) return kInteractif_xml;
439	else if ( input.EndsWith(".txt") ) return kInteractif_ESDList;
440	else if ( input.EndsWith(".root") ) return kLocal;
4060971a	441	} else if (smode == "caf" \|\| smode == "saf") return kProof;
00612ffc	442	return -1;
	443	}
	444
	445	//______________________________________________________________________________
	446	TChain* CreateChainFromCollection(const char *xmlfile)
	447	{
	448	// Create a chain from the collection of tags.
	449	if (!TGrid::Connect("alien://")) return NULL;
	450
	451	TGridCollection* coll = TAlienCollection::Open(xmlfile);
	452	if (!coll) {
	453	Error("CreateChainFromCollection", "Cannot create the AliEn collection");
	454	return NULL;
	455	}
	456
	457	TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE);
	458	AliTagAnalysis *tagAna = new AliTagAnalysis("ESD");
	459	tagAna->ChainGridTags(tagResult);
	460
	461	AliRunTagCuts *runCuts = new AliRunTagCuts();
	462	AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();
	463	AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
	464	AliEventTagCuts *evCuts = new AliEventTagCuts();
	465
	466	// Check if the cuts configuration file was provided
	467	if (!gSystem->AccessPathName("ConfigureCuts.C"))
	468	gROOT->ProcessLine(Form(".x ConfigureCuts.C((AliRunTagCuts)%p, (AliLHCTagCuts)%p, (AliDetectorTagCuts*)%p,"
	469	" (AliEventTagCuts*)%p)", runCuts, lhcCuts, detCuts, evCuts));
	470
	471	TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts);
	472	if (!chain \|\| !chain->GetNtrees()) return NULL;
	473	chain->ls();
	474	return chain;
	475	}
	476
	477	//______________________________________________________________________________
	478	TChain* CreateChainFromFile(const char *rootfile)
	479	{
	480	// Create a chain using the root file.
	481	TChain* chain = new TChain("esdTree");
	482	chain->Add(rootfile);
	483	if (!chain->GetNtrees()) return NULL;
	484	chain->ls();
	485	return chain;
	486	}
	487
	488	//______________________________________________________________________________
	489	TChain* CreateChainFromESDList(const char *esdList)
	490	{
	491	// Create a chain using tags from the run list.
	492	TChain* chain = new TChain("esdTree");
	493	ifstream inFile(esdList);
	494	TString inFileName;
	495	if (inFile.is_open()) {
	496	while (! inFile.eof() ) {
	497	inFileName.ReadLine(inFile,kFALSE);
	498	if(!inFileName.EndsWith(".root")) continue;
	499	chain->Add(inFileName.Data());
	500	}
	501	}
	502	inFile.close();
	503	if (!chain->GetNtrees()) return NULL;
	504	chain->ls();
	505	return chain;
506	}
507
508	//______________________________________________________________________________
509	TChain* CreateChain(Int_t mode, TString input)
510	{
511	printf("*******************************\n");
512	printf("* Getting the Chain *\n");
513	printf("*******************************\n");
514	if(mode == kInteractif_xml) return CreateChainFromCollection(input.Data());
515	else if (mode == kInteractif_ESDList) return CreateChainFromESDList(input.Data());
516	else if (mode == kLocal) return CreateChainFromFile(input.Data());
517	else return NULL;
518	}
519