PWG3/muondep/MuonResolution.C

   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)*ey*ey+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