AliPhysics  ec7afe5 (ec7afe5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MuonResolution.C
Go to the documentation of this file.
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 <TFileCollection.h>
27 #include <TAlienCollection.h>
28 #include <TGridCollection.h>
29 #include <TGridResult.h>
30 
31 // STEER includes
32 #include "AliLog.h"
33 #include "AliCDBManager.h"
34 #include "AliAnalysisManager.h"
35 #include "AliESDInputHandler.h"
36 #include "AliTagAnalysis.h"
37 #include "AliRunTagCuts.h"
38 #include "AliLHCTagCuts.h"
39 #include "AliDetectorTagCuts.h"
40 #include "AliEventTagCuts.h"
41 #include "AliAnalysisDataContainer.h"
42 
43 // PHYSICS includes
44 #include "AliPhysicsSelectionTask.h"
45 #include "AliPhysicsSelection.h"
46 #include "AliBackgroundSelection.h"
47 #include "AliCentralitySelectionTask.h"
49 
50 // MUON includes
51 #include "AliMpCDB.h"
52 #include "AliMpDetElement.h"
53 #include "AliMpDDLStore.h"
54 #include "AliMUONCalibParamND.h"
55 #include "AliMUON2DMap.h"
56 #include "AliMUONTrackerData.h"
57 #include "AliMUONPainterDataRegistry.h"
58 #include "AliMUONTrackerDataWrapper.h"
59 
60 #include "AliMuonEventCuts.h"
61 #include "AliMuonTrackCuts.h"
62 
63 #include "AddTaskMuonResolution.C"
64 
65 #endif
66 
68 Int_t nDE = 200;
69 
70 Bool_t Resume(Int_t &firstStep, Double_t clusterResNB[10], Double_t clusterResB[10],
71  Double_t clusterResNBErr[10], Double_t clusterResBErr[10],
72  Bool_t shiftHalfCh, Double_t halfChShiftNB[20], Double_t halfChShiftB[20],
73  Double_t halfChShiftNBErr[20], Double_t halfChShiftBErr[20],
74  Bool_t shiftDE, Double_t deShiftNB[200], Double_t deShiftB[200],
75  TGraphErrors* clusterResXVsStep[10], TGraphErrors* clusterResYVsStep[10],
76  TGraphErrors* halfChShiftXVsStep[20], TGraphErrors* halfChShiftYVsStep[20]);
77 void LoadAlirootOnProof(TString& aaf, TString rootVersion, TString aliphysicsVersion, Int_t iStep);
78 AliAnalysisTaskMuonResolution* CreateAnalysisTrain(Int_t mode, Int_t iStep, Bool_t selectPhysics, Bool_t selectTrigger,
79  Bool_t matchTrig, Bool_t applyAccCut, Bool_t applyPDCACut,
80  Double_t minMomentum, Double_t minPt, Bool_t isMC, Bool_t correctForSystematics,
81  Int_t extrapMode, Double_t clusterResNB[10], Double_t clusterResB[10],
82  Bool_t shiftHalfCh, Double_t halfChShiftNB[20], Double_t halfChShiftB[20],
83  Bool_t shiftDE, Double_t deShiftNB[200], Double_t deShiftB[200]);
84 Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10],
85  Double_t clusterResNBErr[10], Double_t clusterResBErr[10]);
86 Bool_t AddHalfChShift(Int_t iStep, Double_t halfChShiftNB[20], Double_t halfChShiftB[20], Double_t halfChShiftNBErr[20], Double_t halfChShiftBErr[20]);
87 Bool_t AddDEShift(Int_t iStep, Double_t deShiftNB[200], Double_t deShiftB[200]);
88 void AddMCHViews(TString smode, TFile* file);
89 AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name);
90 Int_t GetMode(TString smode, TString input);
91 TChain* CreateChainFromCollection(const char *xmlfile);
92 TChain* CreateChainFromFile(const char *rootfile);
93 TChain* CreateChainFromESDList(const char *esdList);
95 
96 //______________________________________________________________________________
97 void MuonResolution(TString smode, TString inputFileName, TString rootVersion, TString aliphysicsVersion, Int_t nSteps,
98  Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig, Bool_t applyAccCut, Bool_t applyPDCACut,
99  Double_t minMomentum, Double_t minPt, Bool_t isMC, Bool_t correctForSystematics, Int_t extrapMode,
100  Bool_t shiftHalfCh, Bool_t shiftDE, Int_t nevents)
101 {
103 
104  // timer start...
105  TStopwatch* localTimer = new TStopwatch;
106 
107  // check parameters
108  nSteps = TMath::Max(nSteps,1);
109  if (extrapMode != 0 && extrapMode != 1) {
110  Error("MuonResolution","incorrect extrapolation mode!");
111  return;
112  }
113 
114  // Check runing mode
115  Int_t mode = GetMode(smode, inputFileName);
116  if(mode < 0){
117  Error("MuonResolution","Please provide either an ESD root file, a list of ESDs, a collection of ESDs or a dataset.");
118  return;
119  }
120 
121  // set starting chamber resolution (if -1 they will be loaded from recoParam in the task)
122  Double_t clusterResNB[10] = {-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.};
123  Double_t clusterResB[10] = {-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.};
124  Double_t clusterResNBErr[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
125  Double_t clusterResBErr[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
126  Double_t halfChShiftNB[20] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
127  Double_t halfChShiftB[20] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
128  Double_t halfChShiftNBErr[20] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
129  Double_t halfChShiftBErr[20] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
130  Double_t deShiftNB[200];
131  Double_t deShiftB[200];
132  for (Int_t i=0; i<200; i++) {
133  deShiftNB[i] = 0.;
134  deShiftB[i] = 0.;
135  }
136 
137  // output graphs
138  TMultiGraph* mgClusterResXVsStep = new TMultiGraph("mgClusterResXVsStep","cluster X-resolution versus step;step;#sigma_{X} (cm)");
139  TMultiGraph* mgClusterResYVsStep = new TMultiGraph("mgClusterResYVsStep","cluster Y-resolution versus step;step;#sigma_{Y} (cm)");
140  TGraphErrors* clusterResXVsStep[10];
141  TGraphErrors* clusterResYVsStep[10];
142  for (Int_t i = 0; i < 10; i++) {
143  clusterResXVsStep[i] = new TGraphErrors(nSteps+1);
144  clusterResXVsStep[i]->SetName(Form("gResX_ch%d",i+1));
145  clusterResXVsStep[i]->SetMarkerStyle(kFullDotMedium);
146  clusterResXVsStep[i]->SetMarkerColor(i+1+i/9);
147  mgClusterResXVsStep->Add(clusterResXVsStep[i],"lp");
148 
149  clusterResYVsStep[i] = new TGraphErrors(nSteps+1);
150  clusterResYVsStep[i]->SetName(Form("gResY_ch%d",i+1));
151  clusterResYVsStep[i]->SetMarkerStyle(kFullDotMedium);
152  clusterResYVsStep[i]->SetMarkerColor(i+1+i/9);
153  mgClusterResYVsStep->Add(clusterResYVsStep[i],"lp");
154  }
155  TMultiGraph* mgHalfChShiftXVsStep = new TMultiGraph("mgHalfChShiftXVsStep","half-chamber displacement in X direction versus step;step;#Delta_{X} (cm)");
156  TMultiGraph* mgHalfChShiftYVsStep = new TMultiGraph("mgHalfChShiftYVsStep","half-chamber displacement in Y direction versus step;step;#Delta_{Y} (cm)");
157  TGraphErrors* halfChShiftXVsStep[20];
158  TGraphErrors* halfChShiftYVsStep[20];
159  for (Int_t i = 0; i < 20; i++) {
160  halfChShiftXVsStep[i] = new TGraphErrors(nSteps+1);
161  halfChShiftXVsStep[i]->SetName(Form("gShiftX_hch%d",i+1));
162  halfChShiftXVsStep[i]->SetMarkerStyle(kFullDotMedium);
163  halfChShiftXVsStep[i]->SetMarkerColor(i+1+i/9+i/18);
164  mgHalfChShiftXVsStep->Add(halfChShiftXVsStep[i],"lp");
165  halfChShiftXVsStep[i]->SetPoint(0, 0, halfChShiftNB[i]);
166  halfChShiftXVsStep[i]->SetPointError(0, 0., halfChShiftNBErr[i]);
167 
168  halfChShiftYVsStep[i] = new TGraphErrors(nSteps+1);
169  halfChShiftYVsStep[i]->SetName(Form("gShiftY_hch%d",i+1));
170  halfChShiftYVsStep[i]->SetMarkerStyle(kFullDotMedium);
171  halfChShiftYVsStep[i]->SetMarkerColor(i+1+i/9+i/18);
172  mgHalfChShiftYVsStep->Add(halfChShiftYVsStep[i],"lp");
173  halfChShiftYVsStep[i]->SetPoint(0, 0, halfChShiftB[i]);
174  halfChShiftYVsStep[i]->SetPointError(0, 0., halfChShiftBErr[i]);
175  }
176 
177  // check for old output files
178  Int_t firstStep = 0;
179  char remove = '\0';
180  if (!gSystem->Exec("ls chamberResolution_step*[0-9].root")) {
181  cout<<"Above files already exist in the current directory. [d=delete, r=resume, e=exit] "<<flush;
182  while (remove != 'd' && remove != 'r' && remove != 'e') cin>>remove;
183  if (remove == 'y') gSystem->Exec("rm -f chamberResolution_step*[0-9].root");
184  else if (remove == 'r' && !Resume(firstStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr,
185  shiftHalfCh, halfChShiftNB, halfChShiftB, halfChShiftNBErr, halfChShiftBErr,
186  shiftDE, deShiftNB, deShiftB, clusterResXVsStep, clusterResYVsStep,
187  halfChShiftXVsStep, halfChShiftYVsStep)) return;
188  else if (remove == 'e') return;
189  }
190 
191  // Create input object
192  TObject* inputObj = 0x0;
193  if (mode == kProof) {
194  if (inputFileName.EndsWith(".root")) {
195  TFile *inFile = TFile::Open(inputFileName.Data(),"READ");
196  if (inFile && inFile->IsOpen()) {
197  inputObj = dynamic_cast<TFileCollection*>(inFile->FindObjectAny("dataset"));
198  inFile->Close();
199  }
200  } else inputObj = new TObjString(inputFileName);
201  } else inputObj = CreateChain(mode, inputFileName);
202  if (!inputObj) return;
203 
204  // loop over step
205  for (Int_t iStep = firstStep; iStep < nSteps; iStep++) {
206  cout<<"step "<<iStep+1<<"/"<<nSteps<<endl;
207 
208  // Connect to proof if needed and prepare environment
209  if (mode == kProof) LoadAlirootOnProof(smode, rootVersion, aliphysicsVersion, iStep-firstStep);
210 
211  // create the analysis train
212  AliAnalysisTaskMuonResolution *muonResolution = CreateAnalysisTrain(mode, iStep, selectPhysics, selectTrigger, matchTrig,
213  applyAccCut, applyPDCACut, minMomentum, minPt, isMC,
214  correctForSystematics, extrapMode, clusterResNB, clusterResB,
215  shiftHalfCh, halfChShiftNB, halfChShiftB, shiftDE, deShiftNB, deShiftB);
216  if (!muonResolution) return;
217 
218  // start analysis
219  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
220  if (mgr->InitAnalysis()) {
221  mgr->PrintStatus();
222  if (mode == kProof) {
223  if (inputObj->IsA() == TFileCollection::Class())
224  mgr->StartAnalysis("proof", static_cast<TFileCollection*>(inputObj), nevents);
225  else mgr->StartAnalysis("proof", static_cast<TObjString*>(inputObj)->GetName(), nevents);
226  } else mgr->StartAnalysis("local", static_cast<TChain*>(inputObj), nevents);
227  }
228 
229  // save the summary canvases and mchview display
230  if (muonResolution->GetCanvases()) {
231  TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"UPDATE");
232  if (outFile && outFile->IsOpen()) {
233  outFile->cd();
234  muonResolution->GetCanvases()->Write();
235  AddMCHViews(smode, outFile);
236  outFile->Close();
237  delete outFile;
238  }
239  }
240 
241  // fill graphs with starting resolutions from the task at very first step
242  if (iStep == 0) {
243  muonResolution->GetStartingResolution(clusterResNB, clusterResB);
244  for (Int_t i = 0; i < 10; i++) {
245  clusterResXVsStep[i]->SetPoint(0, 0, clusterResNB[i]);
246  clusterResXVsStep[i]->SetPointError(0, 0., clusterResNBErr[i]);
247  clusterResYVsStep[i]->SetPoint(0, 0, clusterResB[i]);
248  clusterResYVsStep[i]->SetPointError(0, 0., clusterResBErr[i]);
249  }
250  }
251 
252  // read the chamber resolution from the output file
253  if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr)) return;
254 
255  // fill graphs with computed resolutions
256  for (Int_t i = 0; i < 10; i++) {
257  clusterResXVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResNB[i]);
258  clusterResXVsStep[i]->SetPointError(iStep+1, 0., clusterResNBErr[i]);
259  clusterResYVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResB[i]);
260  clusterResYVsStep[i]->SetPointError(iStep+1, 0., clusterResBErr[i]);
261  }
262 
263  // get the half-chamber displacements currently used and add the new measurements from the output file
264  muonResolution->GetHalfChShift(halfChShiftNB, halfChShiftB);
265  if (!AddHalfChShift(iStep, halfChShiftNB, halfChShiftB, halfChShiftNBErr, halfChShiftBErr)) return;
266 
267  // fill graphs with computed displacements
268  for (Int_t i = 0; i < 20; i++) {
269  halfChShiftXVsStep[i]->SetPoint(iStep+1, iStep+1, halfChShiftNB[i]);
270  halfChShiftXVsStep[i]->SetPointError(iStep+1, 0., halfChShiftNBErr[i]);
271  halfChShiftYVsStep[i]->SetPoint(iStep+1, iStep+1, halfChShiftB[i]);
272  halfChShiftYVsStep[i]->SetPointError(iStep+1, 0., halfChShiftBErr[i]);
273  }
274 
275  // get the DE displacements currently used and add the new measurements from the output file
276  muonResolution->GetDEShift(deShiftNB, deShiftB);
277  if (!AddDEShift(iStep, deShiftNB, deShiftB)) return;
278 
279  // clean memory
280  mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis); // to make sure that all output containers are deleted
281  delete mgr;
282  TObject::SetObjectStat(kFALSE);
283 
284  }
285 
286  // copy final results in results.root file
287  gSystem->Exec(Form("cp chamberResolution_step%d.root results.root", nSteps-1));
288 
289  // display convergence of cluster resolution
290  TCanvas* convergence1 = new TCanvas("convergenceRes","convergence of cluster resolution");
291  convergence1->Divide(1,2);
292  convergence1->cd(1);
293  mgClusterResXVsStep->Draw("ap");
294  convergence1->cd(2);
295  mgClusterResYVsStep->Draw("ap");
296 
297  // display convergence of half-chamber displacements
298  TCanvas* convergence2 = new TCanvas("convergenceShift","convergence of half-chamber displacements");
299  convergence2->Divide(1,2);
300  convergence2->cd(1);
301  mgHalfChShiftXVsStep->Draw("ap");
302  convergence2->cd(2);
303  mgHalfChShiftYVsStep->Draw("ap");
304 
305  // save convergence plots
306  TFile* outFile = TFile::Open("results.root","UPDATE");
307  if (!outFile || !outFile->IsOpen()) return;
308  outFile->cd();
309  mgClusterResXVsStep->Write();
310  mgClusterResYVsStep->Write();
311  convergence1->Write();
312  mgHalfChShiftXVsStep->Write();
313  mgHalfChShiftYVsStep->Write();
314  convergence2->Write();
315  outFile->Close();
316  delete outFile;
317 
318  // print final half-chamber displacements
319  printf("\nhalf-chamber total displacements:\n");
320  printf(" - non-bending:");
321  for (Int_t i = 0; i < 20; i++) printf((i==0)?" %6.4f":", %6.4f", halfChShiftNB[i]);
322  printf("\n - bending:");
323  for (Int_t i = 0; i < 20; i++) printf((i==0)?" %6.4f":", %6.4f", halfChShiftB[i]);
324  printf("\n\n");
325 
326  // print final DE displacements
327  printf("\nDE total displacements:\n");
328  printf(" - non-bending:");
329  for (Int_t i = 0; i < nDE; i++) printf((i==0)?" %6.4f":", %6.4f", deShiftNB[i]);
330  printf("\n - bending:");
331  for (Int_t i = 0; i < nDE; i++) printf((i==0)?" %6.4f":", %6.4f", deShiftB[i]);
332  printf("\n\n");
333 
334  // ...timer stop
335  localTimer->Stop();
336  localTimer->Print();
337  delete localTimer;
338 
339 }
340 
341 //______________________________________________________________________________
342 Bool_t Resume(Int_t &firstStep, Double_t clusterResNB[10], Double_t clusterResB[10],
343  Double_t clusterResNBErr[10], Double_t clusterResBErr[10],
344  Bool_t shiftHalfCh, Double_t halfChShiftNB[20], Double_t halfChShiftB[20],
345  Double_t halfChShiftNBErr[20], Double_t halfChShiftBErr[20],
346  Bool_t shiftDE, Double_t deShiftNB[200], Double_t deShiftB[200],
347  TGraphErrors* clusterResXVsStep[10], TGraphErrors* clusterResYVsStep[10],
348  TGraphErrors* halfChShiftXVsStep[20], TGraphErrors* halfChShiftYVsStep[20])
349 {
351 
352  while (kTRUE) {
353 
354  // Get the step to restart from
355  cout<<"From which step (included) you want to resume? [#, e=exit] "<<flush;
356  TString step = "";
357  do {step.Gets(stdin,kTRUE);} while (!step.IsDigit() && step != "e");
358  if (step == "e") return kFALSE;
359  firstStep = step.Atoi();
360 
361  // restart from scratch if requested
362  if (firstStep == 0) {
363  gSystem->Exec("rm -f chamberResolution_step*[0-9].root");
364  return kTRUE;
365  }
366 
367  // look for results from the previous step
368  if (gSystem->AccessPathName(Form("chamberResolution_step%d.root", firstStep-1))) {
369  cout<<"No result found from the previous step ("<<firstStep-1<<"). Unable to resume from step "<<firstStep<<endl;
370  continue;
371  }
372 
373  // fill graph with starting resolutions
374  for (Int_t i = 0; i < 10; i++) {
375  clusterResXVsStep[i]->SetPoint(0, 0, clusterResNB[i]);
376  clusterResXVsStep[i]->SetPointError(0, 0., clusterResNBErr[i]);
377  clusterResYVsStep[i]->SetPoint(0, 0, clusterResB[i]);
378  clusterResYVsStep[i]->SetPointError(0, 0., clusterResBErr[i]);
379  }
380 
381  // loop over previous steps
382  Bool_t missingInfo = kFALSE;
383  for (Int_t iStep = 0; iStep < firstStep; iStep++) {
384 
385  // read the chamber resolution from the output file
386  if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr) && iStep == firstStep-1) {
387  missingInfo = kTRUE;
388  break;
389  }
390 
391  // fill graphs with computed resolutions
392  for (Int_t i = 0; i < 10; i++) {
393  clusterResXVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResNB[i]);
394  clusterResXVsStep[i]->SetPointError(iStep+1, 0., clusterResNBErr[i]);
395  clusterResYVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResB[i]);
396  clusterResYVsStep[i]->SetPointError(iStep+1, 0., clusterResBErr[i]);
397  }
398 
399  // reset the half-chamber displacements if not used and add the new measurements from the output file
400  if (!shiftHalfCh) for (Int_t i=0; i<20; i++) {
401  halfChShiftNB[i] = 0.; halfChShiftB[i] = 0.;
402  halfChShiftNBErr[i] = 0.; halfChShiftBErr[i] = 0.;
403  }
404  if (!AddHalfChShift(iStep, halfChShiftNB, halfChShiftB, halfChShiftNBErr, halfChShiftBErr) && shiftHalfCh) {
405  missingInfo = kTRUE;
406  break;
407  }
408 
409  // fill graphs with computed displacements
410  for (Int_t i = 0; i < 20; i++) {
411  halfChShiftXVsStep[i]->SetPoint(iStep+1, iStep+1, halfChShiftNB[i]);
412  halfChShiftXVsStep[i]->SetPointError(iStep+1, 0., halfChShiftNBErr[i]);
413  halfChShiftYVsStep[i]->SetPoint(iStep+1, iStep+1, halfChShiftB[i]);
414  halfChShiftYVsStep[i]->SetPointError(iStep+1, 0., halfChShiftBErr[i]);
415  }
416 
417  // add the new measurements of DE displacements from the output file if in use
418  if (shiftDE && !AddDEShift(iStep, deShiftNB, deShiftB)) {
419  missingInfo = kTRUE;
420  break;
421  }
422 
423  }
424 
425  // check if missing important results from previous steps
426  if (missingInfo) continue;
427 
428  // keep previous steps and remove the others
429  gSystem->Exec("mkdir __TMP__");
430  for (Int_t iStep = 0; iStep < firstStep; iStep++)
431  if (!gSystem->AccessPathName(Form("chamberResolution_step%d.root", iStep)))
432  gSystem->Exec(Form("mv chamberResolution_step%d.root __TMP__", iStep));
433  gSystem->Exec("rm -f chamberResolution_step*[0-9].root");
434  gSystem->Exec("mv __TMP__/chamberResolution_step*[0-9].root .");
435  gSystem->Exec("rm -rf __TMP__");
436 
437  return kTRUE;
438  }
439 
440 }
441 
442 //______________________________________________________________________________
443 void LoadAlirootOnProof(TString& aaf, TString rootVersion, TString aliphysicsVersion, Int_t iStep)
444 {
446 
447  // set general environment and close previous session
448  if (iStep == 0) gEnv->SetValue("XSec.GSI.DelegProxy","2");
449  else gProof->Close("s");
450 
451  // connect
452  if (aaf == "saf3") TProof::Open("pod://");
453  else {
454  TString location = (aaf == "caf") ? "alice-caf.cern.ch" : "nansafmaster2.in2p3.fr"; //"localhost:1093"
455  TString nWorkers = (aaf == "caf") ? "workers=80" : ""; //"workers=3x"
456  TString user = (gSystem->Getenv("alien_API_USER") == NULL) ? "" : Form("%s@",gSystem->Getenv("alien_API_USER"));
457  TProof::Mgr(Form("%s%s",user.Data(), location.Data()))->SetROOTVersion(Form("VO_ALICE@ROOT::%s",rootVersion.Data()));
458  TProof::Open(Form("%s%s/?N",user.Data(), location.Data()), nWorkers.Data());
459  }
460  if (!gProof) return;
461 
462  // set environment and load libraries on workers
463  TList* list = new TList();
464  list->Add(new TNamed("ALIROOT_MODE", "base"));
465  list->Add(new TNamed("ALIROOT_ENABLE_ALIEN", "1"));
466  if (aaf == "saf3") {
467  TString home = gSystem->Getenv("HOME");
468  gProof->UploadPackage(Form("%s/AliceVaf.par", home.Data()));
469  gProof->EnablePackage(Form("%s/AliceVaf.par", home.Data()), list, (iStep!=0));
470  } else gProof->EnablePackage(Form("VO_ALICE@AliPhysics::%s",aliphysicsVersion.Data()), list, kTRUE);
471 // gProof->UploadPackage("$ALICE_PHYSICS/PARfiles/PWGPPMUONdep.par");
472 // gProof->EnablePackage("$ALICE_PHYSICS/PARfiles/PWGPPMUONdep.par", kTRUE);
473 // gProof->UploadPackage("PWGPPMUONdep.par");
474 // gProof->EnablePackage("PWGPPMUONdep.par", kTRUE);
475 
476 }
477 
478 //______________________________________________________________________________
480  Bool_t matchTrig, Bool_t applyAccCut, Bool_t applyPDCACut,
481  Double_t minMomentum, Double_t minPt, Bool_t isMC, Bool_t correctForSystematics,
482  Int_t extrapMode, Double_t clusterResNB[10], Double_t clusterResB[10],
483  Bool_t shiftHalfCh, Double_t halfChShiftNB[20], Double_t halfChShiftB[20],
484  Bool_t shiftDE, Double_t deShiftNB[200], Double_t deShiftB[200])
485 {
487 
488  // Create the analysis manager
489  AliAnalysisManager *mgr = new AliAnalysisManager("MuonResolutionAnalysis");
490  //mgr->SetNSysInfo(100);
491  //mgr->SetDebugLevel(3);
492 
493  // ESD input handler
494  AliESDInputHandler* esdH = new AliESDInputHandler();
495  esdH->SetReadFriends(kFALSE);
496  esdH->SetInactiveBranches("*");
497  esdH->SetActiveBranches("MuonTracks MuonClusters MuonPads AliESDRun. AliESDHeader. AliMultiplicity. AliESDFMD. AliESDVZERO. SPDVertex. PrimaryVertex. AliESDZDC.");
498  mgr->SetInputEventHandler(esdH);
499 
500  // event selection
501  UInt_t eventSelectionMask = 0;
502  if (selectPhysics) {
503  gROOT->LoadMacro("$ALICE_PHYSICS/OADB/macros/AddTaskPhysicsSelection.C");
504  AliPhysicsSelectionTask* physicsSelection = reinterpret_cast<AliPhysicsSelectionTask*>(gROOT->ProcessLineSync(TString::Format("AddTaskPhysicsSelection(%d)", isMC)));
505  if (!physicsSelection) {
506  Error("CreateAnalysisTrain","AliPhysicsSelectionTask not created!");
507  return 0x0;
508  }
509  //if (!isMC) physicsSelection->GetPhysicsSelection()->SetUseBXNumbers(kFALSE); // Needed to merge runs with different running scheme
510  eventSelectionMask |= AliMuonEventCuts::kPhysicsSelected;
511  }
512  if (selectTrigger) eventSelectionMask |= AliMuonEventCuts::kSelectedTrig;
513 
514  // multiplicity/centrality selection
515  gROOT->LoadMacro("$ALICE_PHYSICS/OADB/COMMON/MULTIPLICITY/macros/AddTaskMultSelection.C");
516  if (!gROOT->ProcessLineFast("AddTaskMultSelection()")) {
517  Error("CreateAnalysisTrain","AliCentralitySelectionTask not created!");
518  return 0x0;
519  }
520 
521  // Muon Resolution analysis
522  TString outputFileName = Form("chamberResolution_step%d.root", iStep);
523  AliAnalysisManager::SetCommonFileName(outputFileName.Data());
524  AliAnalysisTaskMuonResolution *muonResolution = AddTaskMuonResolution(minMomentum, minPt, correctForSystematics, extrapMode);
525  if (!muonResolution) {
526  Error("CreateAnalysisTrain","AliAnalysisTaskMuonResolution not created!");
527  return 0x0;
528  }
529  /*if (mode == kLocal) muonResolution->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
530  else muonResolution->SetDefaultStorage("raw://");*/
531  muonResolution->SetDefaultStorage("local:///cvmfs/alice-ocdb.cern.ch/calibration/data/2016/OCDB");
532  if (mode != kProof) muonResolution->ShowProgressBar();
533  muonResolution->PrintClusterRes(kTRUE, kTRUE);
534  muonResolution->SetStartingResolution(clusterResNB, clusterResB);
535  muonResolution->RemoveMonoCathodClusters(kTRUE, kFALSE);
536 // muonResolution->FitResiduals(kFALSE);
537 // muonResolution->ImproveTracks(kTRUE);
538 // muonResolution->ReAlign("", -1, -1, "alien://folder=/alice/cern.ch/user/h/hupereir/CDB/LHC16_realign");
539 // muonResolution->ReAlign("", 1, 0, "");
540 // muonResolution->SetAlignStorage("", 5);
541 // muonResolution->SetMuonSign(-1);
542 
543  if (shiftHalfCh) {
544  muonResolution->SetHalfChShift(halfChShiftNB, halfChShiftB);
545  muonResolution->ShiftHalfCh();
546  muonResolution->PrintHalfChShift();
547  }
548  if (shiftDE) {
549  muonResolution->SetDEShift(deShiftNB, deShiftB);
550  muonResolution->ShiftDE();
551  muonResolution->PrintDEShift();
552  }
553 
554  if (eventSelectionMask != 0) {
555  AliMuonEventCuts eventCuts("muEventCuts", "muEventCuts");
556  eventCuts.SetFilterMask(eventSelectionMask);
557  if (selectPhysics) eventCuts.SetPhysicsSelectionMask(AliVEvent::kAny);
558 // if (selectPhysics) eventCuts.SetPhysicsSelectionMask(AliVEvent::kMuonUnlikeLowPt7);
559  if (selectTrigger) eventCuts.SetTrigClassPatterns(eventCuts.GetDefaultTrigClassPatterns());
560  muonResolution->SetMuonEventCuts(eventCuts);
561  }
562 
563  UInt_t trackSelectionMask = 0;
564  if (matchTrig) trackSelectionMask |= AliMuonTrackCuts::kMuMatchLpt;
565  if (applyAccCut) trackSelectionMask |= AliMuonTrackCuts::kMuEta | AliMuonTrackCuts::kMuThetaAbs;
566  if (applyPDCACut) trackSelectionMask |= AliMuonTrackCuts::kMuPdca;
567  if (trackSelectionMask != 0) {
568  AliMuonTrackCuts trackCuts("muTrackCuts", "muTrackCuts");
569  trackCuts.SetAllowDefaultParams();
570  if (isMC) trackCuts.SetIsMC();
571  trackCuts.SetFilterMask(trackSelectionMask);
572  muonResolution->SetMuonTrackCuts(trackCuts);
573  }
574 
575  return muonResolution;
576 
577 }
578 
579 //______________________________________________________________________________
580 Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10])
581 {
583 
584  TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
585 
586  if (!outFile || !outFile->IsOpen()) {
587  Error("GetChamberResolution","output file does not exist!");
588  return kFALSE;
589  }
590 
591  TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
592  TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0;
593  TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerChSigma")) : 0x0;
594 
595  if (!gCombinedResidualXPerChSigma || !gCombinedResidualYPerChSigma) {
596  Error("GetChamberResolution","resolution graphs do not exist!");
597  return kFALSE;
598  }
599 
600  Double_t dummy;
601  for (Int_t i = 0; i < 10; i++) {
602  gCombinedResidualXPerChSigma->GetPoint(i, dummy, clusterResNB[i]);
603  gCombinedResidualYPerChSigma->GetPoint(i, dummy, clusterResB[i]);
604  clusterResNBErr[i] = gCombinedResidualXPerChSigma->GetErrorY(i);
605  clusterResBErr[i] = gCombinedResidualYPerChSigma->GetErrorY(i);
606  }
607 
608  outFile->Close();
609  //delete outFile;
610 
611  return kTRUE;
612 }
613 
614 //______________________________________________________________________________
615 Bool_t AddHalfChShift(Int_t iStep, Double_t halfChShiftNB[20], Double_t halfChShiftB[20], Double_t halfChShiftNBErr[20], Double_t halfChShiftBErr[20])
616 {
618 
619  TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
620 
621  if (!outFile || !outFile->IsOpen()) {
622  Error("AddHalfChShift","output file does not exist!");
623  return kFALSE;
624  }
625 
626  TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
627  TGraphErrors* gResidualXPerHalfChMean = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerHalfChMean_ClusterIn")) : 0x0;
628  TGraphErrors* gResidualYPerHalfChMean = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerHalfChMean_ClusterIn")) : 0x0;
629 
630  if (!gResidualXPerHalfChMean || !gResidualYPerHalfChMean) {
631  Error("AddHalfChShift","half-chamber shift graphs do not exist!");
632  return kFALSE;
633  }
634 
635  Double_t dummy, dx, dy;
636  for (Int_t i = 0; i < 20; i++) {
637  gResidualXPerHalfChMean->GetPoint(i, dummy, dx);
638  halfChShiftNB[i] += dx;
639  halfChShiftNBErr[i] = gResidualXPerHalfChMean->GetErrorY(i);
640  gResidualYPerHalfChMean->GetPoint(i, dummy, dy);
641  halfChShiftB[i] += dy;
642  halfChShiftBErr[i] = gResidualYPerHalfChMean->GetErrorY(i);
643  }
644 
645  outFile->Close();
646  //delete outFile;
647 
648  return kTRUE;
649 }
650 
651 //______________________________________________________________________________
652 Bool_t AddDEShift(Int_t iStep, Double_t deShiftNB[200], Double_t deShiftB[200])
653 {
655 
656  TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
657 
658  if (!outFile || !outFile->IsOpen()) {
659  Error("AddDEShift","output file does not exist!");
660  return kFALSE;
661  }
662 
663  TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
664  TGraphErrors* gResidualXPerDEMean = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerDEMean_ClusterIn")) : 0x0;
665  TGraphErrors* gResidualYPerDEMean = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerDEMean_ClusterIn")) : 0x0;
666 
667  if (!gResidualXPerDEMean || !gResidualYPerDEMean) {
668  Error("AddDEShift","DE shift graphs do not exist!");
669  return kFALSE;
670  }
671 
672  Double_t dummy, dx, dy;
673  nDE = gResidualXPerDEMean->GetN();
674  for (Int_t i = 0; i < nDE; i++) {
675  gResidualXPerDEMean->GetPoint(i, dummy, dx);
676  deShiftNB[i] += dx;
677  gResidualYPerDEMean->GetPoint(i, dummy, dy);
678  deShiftB[i] += dy;
679  }
680 
681  outFile->Close();
682  //delete outFile;
683 
684  return kTRUE;
685 }
686 
687 //______________________________________________________________________________
688 void AddMCHViews(TString smode, TFile* file)
689 {
691 
692  if ( ! AliMpDDLStore::Instance(false) )
693  {
694  Warning("AddMCHViews","mapping was not loaded. Loading it from $ALICE_ROOT/OCDB");
695  if (smode == "saf3") AliCDBManager::Instance()->SetDefaultStorage("raw://");
696  else AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
697  AliCDBManager::Instance()->SetRun(999999999);
698  }
699 
700  AliMpCDB::LoadAll();
701 
702  TObjArray* summary = static_cast<TObjArray*>(file->FindObjectAny("ChamberRes"));
703  if (!summary) {
704  Error("AddMCHViews","resolution graphs do not exist!");
705  return;
706  }
707 
708  TGraphErrors* g = 0x0;
709  g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerDESigma"));
710  if (g) {
711  file->cd();
712  AliMUONTrackerData* data = ConvertGraph(*g, "resoX");
713  data->Write();
714  delete data;
715  }
716 
717  g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerDESigma"));
718  if (g) {
719  file->cd();
720  AliMUONTrackerData* data = ConvertGraph(*g, "resoY");
721  data->Write();
722  delete data;
723  }
724 
725  g = static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerDEMean_ClusterOut"));
726  if (g) {
727  file->cd();
728  AliMUONTrackerData* data = ConvertGraph(*g, "shiftX");
729  data->Write();
730  delete data;
731  }
732 
733  g = static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerDEMean_ClusterOut"));
734  if (g) {
735  file->cd();
736  AliMUONTrackerData* data = ConvertGraph(*g, "shiftY");
737  data->Write();
738  delete data;
739  }
740 }
741 
742 //______________________________________________________________________________
743 AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name)
744 {
746 
747  AliMUON2DMap deValues(kFALSE);
748 
749  for ( Int_t i = 0 ; i < g.GetN(); ++i )
750  {
751  double y = g.GetY()[i];
752  double ey = g.GetEY()[i];
753  int detElemId;
754  sscanf(g.GetXaxis()->GetBinLabel(i+1),"%d",&detElemId);
755 
756  AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
757 
758  AliMUONVCalibParam* param = new AliMUONCalibParamND(5, 1, detElemId, 0);
759 
760  Double_t sumn = 1000.0;
761  Double_t sumw = sumn*y;
762  Double_t sumw2 = (sumn-1)*ey*ey+sumw*sumw/sumn;
763 
764  param->SetValueAsDouble(0,0,sumw);
765  param->SetValueAsDouble(0,1,sumw2);
766  param->SetValueAsDouble(0,2,sumn);
767  param->SetValueAsDouble(0,3,de->NofChannels());
768  param->SetValueAsDouble(0,4,1);
769 
770  deValues.Add(param);
771  }
772 
773  AliMUONTrackerData* data = new AliMUONTrackerData(name,name,deValues,1);
774  data->SetDimensionName(0,name);
775 
776  return data;
777 }
778 
779 //______________________________________________________________________________
781 {
782  if (smode == "local") {
783  if ( input.EndsWith(".xml") ) return kInteractif_xml;
784  else if ( input.EndsWith(".txt") ) return kInteractif_ESDList;
785  else if ( input.EndsWith(".root") ) return kLocal;
786  } else if (smode == "caf" || smode == "saf" || smode == "saf3") return kProof;
787  return -1;
788 }
789 
790 //______________________________________________________________________________
791 TChain* CreateChainFromCollection(const char *xmlfile)
792 {
793  // Create a chain from the collection of tags.
794  if (!TGrid::Connect("alien://")) return NULL;
795 
796  TGridCollection* coll = TAlienCollection::Open(xmlfile);
797  if (!coll) {
798  Error("CreateChainFromCollection", "Cannot create the AliEn collection");
799  return NULL;
800  }
801 
802  TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE);
803  AliTagAnalysis *tagAna = new AliTagAnalysis("ESD");
804  tagAna->ChainGridTags(tagResult);
805 
806  AliRunTagCuts *runCuts = new AliRunTagCuts();
807  AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();
808  AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
809  AliEventTagCuts *evCuts = new AliEventTagCuts();
810 
811  // Check if the cuts configuration file was provided
812  if (!gSystem->AccessPathName("ConfigureCuts.C"))
813  gROOT->ProcessLine(Form(".x ConfigureCuts.C((AliRunTagCuts*)%p, (AliLHCTagCuts*)%p, (AliDetectorTagCuts*)%p,"
814  " (AliEventTagCuts*)%p)", runCuts, lhcCuts, detCuts, evCuts));
815 
816  TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts);
817  if (!chain || !chain->GetNtrees()) return NULL;
818  chain->ls();
819  return chain;
820 }
821 
822 //______________________________________________________________________________
823 TChain* CreateChainFromFile(const char *rootfile)
824 {
825  // Create a chain using the root file.
826  TChain* chain = new TChain("esdTree");
827  chain->Add(rootfile);
828  if (!chain->GetNtrees()) return NULL;
829  chain->ls();
830  return chain;
831 }
832 
833 //______________________________________________________________________________
834 TChain* CreateChainFromESDList(const char *esdList)
835 {
836  // Create a chain using tags from the run list.
837  TChain* chain = new TChain("esdTree");
838  ifstream inFile(esdList);
839  TString inFileName;
840  if (inFile.is_open()) {
841  while (! inFile.eof() ) {
842  inFileName.ReadLine(inFile,kFALSE);
843  if(!inFileName.EndsWith(".root")) continue;
844  chain->Add(inFileName.Data());
845  }
846  }
847  inFile.close();
848  if (!chain->GetNtrees()) return NULL;
849  chain->ls();
850  return chain;
851 }
852 
853 //______________________________________________________________________________
855 {
856  printf("*******************************\n");
857  printf("*** Getting the Chain ***\n");
858  printf("*******************************\n");
859  if(mode == kInteractif_xml) return CreateChainFromCollection(input.Data());
860  else if (mode == kInteractif_ESDList) return CreateChainFromESDList(input.Data());
861  else if (mode == kLocal) return CreateChainFromFile(input.Data());
862  else return NULL;
863 }
864 
Bool_t AddDEShift(Int_t iStep, Double_t deShiftNB[200], Double_t deShiftB[200])
double Double_t
Definition: External.C:58
Int_t GetMode(TString smode, TString input)
TChain * CreateChain(Int_t mode, TString input)
Bool_t Resume(Int_t &firstStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10], Bool_t shiftHalfCh, Double_t halfChShiftNB[20], Double_t halfChShiftB[20], Double_t halfChShiftNBErr[20], Double_t halfChShiftBErr[20], Bool_t shiftDE, Double_t deShiftNB[200], Double_t deShiftB[200], TGraphErrors *clusterResXVsStep[10], TGraphErrors *clusterResYVsStep[10], TGraphErrors *halfChShiftXVsStep[20], TGraphErrors *halfChShiftYVsStep[20])
TSystem * gSystem
TList * list
void LoadAlirootOnProof(TString &aaf, TString rootVersion, TString aliphysicsVersion, Int_t iStep)
TChain * CreateChainFromESDList(const char *esdList)
TObjArray * GetCanvases()
return the list of summary canvases
Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10])
Int_t nDE
TChain * CreateChainFromCollection(const char *xmlfile)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
AliAnalysisTaskMuonResolution * AddTaskMuonResolution(Double_t minMomentum=0., Double_t minPt=0., Bool_t correctForSystematics=kTRUE, Int_t extrapMode=1)
Int_t mode
Definition: anaM.C:41
void GetHalfChShift(Double_t valNB[20], Double_t valB[20]) const
void MuonResolution(TString smode, TString inputFileName, TString rootVersion, TString aliphysicsVersion, Int_t nSteps, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig, Bool_t applyAccCut, Bool_t applyPDCACut, Double_t minMomentum, Double_t minPt, Bool_t isMC, Bool_t correctForSystematics, Int_t extrapMode, Bool_t shiftHalfCh, Bool_t shiftDE, Int_t nevents)
void AddMCHViews(TString smode, TFile *file)
void GetStartingResolution(Double_t valNB[10], Double_t valB[10]) const
Bool_t isMC
TChain * CreateChainFromFile(const char *rootfile)
Muon spectrometer resolution.
TFile * file
Int_t nevents[nsamples]
Bool_t AddHalfChShift(Int_t iStep, Double_t halfChShiftNB[20], Double_t halfChShiftB[20], Double_t halfChShiftNBErr[20], Double_t halfChShiftBErr[20])
bool Bool_t
Definition: External.C:53
void GetDEShift(Double_t valNB[200], Double_t valB[200]) const
AliAnalysisTaskMuonResolution * CreateAnalysisTrain(Int_t mode, Int_t iStep, Bool_t selectPhysics, Bool_t selectTrigger, Bool_t matchTrig, Bool_t applyAccCut, Bool_t applyPDCACut, Double_t minMomentum, Double_t minPt, Bool_t isMC, Bool_t correctForSystematics, Int_t extrapMode, Double_t clusterResNB[10], Double_t clusterResB[10], Bool_t shiftHalfCh, Double_t halfChShiftNB[20], Double_t halfChShiftB[20], Bool_t shiftDE, Double_t deShiftNB[200], Double_t deShiftB[200])
AliMUONTrackerData * ConvertGraph(TGraphErrors &g, const char *name)