AliPhysics  32b88a8 (32b88a8)
 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/2015/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/LHC15_realign_all_4");
539 // muonResolution->ReAlign("", 6, -1, "");
540 // muonResolution->SetAlignStorage("", 5);
541 
542  if (shiftHalfCh) {
543  muonResolution->SetHalfChShift(halfChShiftNB, halfChShiftB);
544  muonResolution->ShiftHalfCh();
545  muonResolution->PrintHalfChShift();
546  }
547  if (shiftDE) {
548  muonResolution->SetDEShift(deShiftNB, deShiftB);
549  muonResolution->ShiftDE();
550  muonResolution->PrintDEShift();
551  }
552 
553  if (eventSelectionMask != 0) {
554  AliMuonEventCuts eventCuts("muEventCuts", "muEventCuts");
555  eventCuts.SetFilterMask(eventSelectionMask);
556  if (selectPhysics) eventCuts.SetPhysicsSelectionMask(AliVEvent::kAny);
557 // if (selectPhysics) eventCuts.SetPhysicsSelectionMask(AliVEvent::kMuonUnlikeLowPt7);
558  if (selectTrigger) eventCuts.SetTrigClassPatterns(eventCuts.GetDefaultTrigClassPatterns());
559  muonResolution->SetMuonEventCuts(eventCuts);
560  }
561 
562  UInt_t trackSelectionMask = 0;
563  if (matchTrig) trackSelectionMask |= AliMuonTrackCuts::kMuMatchLpt;
564  if (applyAccCut) trackSelectionMask |= AliMuonTrackCuts::kMuEta | AliMuonTrackCuts::kMuThetaAbs;
565  if (applyPDCACut) trackSelectionMask |= AliMuonTrackCuts::kMuPdca;
566  if (trackSelectionMask != 0) {
567  AliMuonTrackCuts trackCuts("muTrackCuts", "muTrackCuts");
568  trackCuts.SetAllowDefaultParams();
569  if (isMC) trackCuts.SetIsMC();
570  trackCuts.SetFilterMask(trackSelectionMask);
571  muonResolution->SetMuonTrackCuts(trackCuts);
572  }
573 
574  return muonResolution;
575 
576 }
577 
578 //______________________________________________________________________________
579 Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10])
580 {
582 
583  TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
584 
585  if (!outFile || !outFile->IsOpen()) {
586  Error("GetChamberResolution","output file does not exist!");
587  return kFALSE;
588  }
589 
590  TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
591  TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0;
592  TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerChSigma")) : 0x0;
593 
594  if (!gCombinedResidualXPerChSigma || !gCombinedResidualYPerChSigma) {
595  Error("GetChamberResolution","resolution graphs do not exist!");
596  return kFALSE;
597  }
598 
599  Double_t dummy;
600  for (Int_t i = 0; i < 10; i++) {
601  gCombinedResidualXPerChSigma->GetPoint(i, dummy, clusterResNB[i]);
602  gCombinedResidualYPerChSigma->GetPoint(i, dummy, clusterResB[i]);
603  clusterResNBErr[i] = gCombinedResidualXPerChSigma->GetErrorY(i);
604  clusterResBErr[i] = gCombinedResidualYPerChSigma->GetErrorY(i);
605  }
606 
607  outFile->Close();
608  //delete outFile;
609 
610  return kTRUE;
611 }
612 
613 //______________________________________________________________________________
614 Bool_t AddHalfChShift(Int_t iStep, Double_t halfChShiftNB[20], Double_t halfChShiftB[20], Double_t halfChShiftNBErr[20], Double_t halfChShiftBErr[20])
615 {
617 
618  TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
619 
620  if (!outFile || !outFile->IsOpen()) {
621  Error("AddHalfChShift","output file does not exist!");
622  return kFALSE;
623  }
624 
625  TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
626  TGraphErrors* gResidualXPerHalfChMean = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerHalfChMean_ClusterIn")) : 0x0;
627  TGraphErrors* gResidualYPerHalfChMean = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerHalfChMean_ClusterIn")) : 0x0;
628 
629  if (!gResidualXPerHalfChMean || !gResidualYPerHalfChMean) {
630  Error("AddHalfChShift","half-chamber shift graphs do not exist!");
631  return kFALSE;
632  }
633 
634  Double_t dummy, dx, dy;
635  for (Int_t i = 0; i < 20; i++) {
636  gResidualXPerHalfChMean->GetPoint(i, dummy, dx);
637  halfChShiftNB[i] += dx;
638  halfChShiftNBErr[i] = gResidualXPerHalfChMean->GetErrorY(i);
639  gResidualYPerHalfChMean->GetPoint(i, dummy, dy);
640  halfChShiftB[i] += dy;
641  halfChShiftBErr[i] = gResidualYPerHalfChMean->GetErrorY(i);
642  }
643 
644  outFile->Close();
645  //delete outFile;
646 
647  return kTRUE;
648 }
649 
650 //______________________________________________________________________________
651 Bool_t AddDEShift(Int_t iStep, Double_t deShiftNB[200], Double_t deShiftB[200])
652 {
654 
655  TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
656 
657  if (!outFile || !outFile->IsOpen()) {
658  Error("AddDEShift","output file does not exist!");
659  return kFALSE;
660  }
661 
662  TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
663  TGraphErrors* gResidualXPerDEMean = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerDEMean_ClusterIn")) : 0x0;
664  TGraphErrors* gResidualYPerDEMean = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerDEMean_ClusterIn")) : 0x0;
665 
666  if (!gResidualXPerDEMean || !gResidualYPerDEMean) {
667  Error("AddDEShift","DE shift graphs do not exist!");
668  return kFALSE;
669  }
670 
671  Double_t dummy, dx, dy;
672  nDE = gResidualXPerDEMean->GetN();
673  for (Int_t i = 0; i < nDE; i++) {
674  gResidualXPerDEMean->GetPoint(i, dummy, dx);
675  deShiftNB[i] += dx;
676  gResidualYPerDEMean->GetPoint(i, dummy, dy);
677  deShiftB[i] += dy;
678  }
679 
680  outFile->Close();
681  //delete outFile;
682 
683  return kTRUE;
684 }
685 
686 //______________________________________________________________________________
687 void AddMCHViews(TString smode, TFile* file)
688 {
690 
691  if ( ! AliMpDDLStore::Instance(false) )
692  {
693  Warning("AddMCHViews","mapping was not loaded. Loading it from $ALICE_ROOT/OCDB");
694  if (smode == "saf3") AliCDBManager::Instance()->SetDefaultStorage("raw://");
695  else AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
696  AliCDBManager::Instance()->SetRun(999999999);
697  }
698 
699  AliMpCDB::LoadAll();
700 
701  TObjArray* summary = static_cast<TObjArray*>(file->FindObjectAny("ChamberRes"));
702  if (!summary) {
703  Error("AddMCHViews","resolution graphs do not exist!");
704  return;
705  }
706 
707  TGraphErrors* g = 0x0;
708  g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerDESigma"));
709  if (g) {
710  file->cd();
711  AliMUONTrackerData* data = ConvertGraph(*g, "resoX");
712  data->Write();
713  delete data;
714  }
715 
716  g = static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerDESigma"));
717  if (g) {
718  file->cd();
719  AliMUONTrackerData* data = ConvertGraph(*g, "resoY");
720  data->Write();
721  delete data;
722  }
723 
724  g = static_cast<TGraphErrors*>(summary->FindObject("gResidualXPerDEMean_ClusterOut"));
725  if (g) {
726  file->cd();
727  AliMUONTrackerData* data = ConvertGraph(*g, "shiftX");
728  data->Write();
729  delete data;
730  }
731 
732  g = static_cast<TGraphErrors*>(summary->FindObject("gResidualYPerDEMean_ClusterOut"));
733  if (g) {
734  file->cd();
735  AliMUONTrackerData* data = ConvertGraph(*g, "shiftY");
736  data->Write();
737  delete data;
738  }
739 }
740 
741 //______________________________________________________________________________
742 AliMUONTrackerData* ConvertGraph(TGraphErrors& g, const char* name)
743 {
745 
746  AliMUON2DMap deValues(kFALSE);
747 
748  for ( Int_t i = 0 ; i < g.GetN(); ++i )
749  {
750  double y = g.GetY()[i];
751  double ey = g.GetEY()[i];
752  int detElemId;
753  sscanf(g.GetXaxis()->GetBinLabel(i+1),"%d",&detElemId);
754 
755  AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
756 
757  AliMUONVCalibParam* param = new AliMUONCalibParamND(5, 1, detElemId, 0);
758 
759  Double_t sumn = 1000.0;
760  Double_t sumw = sumn*y;
761  Double_t sumw2 = (sumn-1)*ey*ey+sumw*sumw/sumn;
762 
763  param->SetValueAsDouble(0,0,sumw);
764  param->SetValueAsDouble(0,1,sumw2);
765  param->SetValueAsDouble(0,2,sumn);
766  param->SetValueAsDouble(0,3,de->NofChannels());
767  param->SetValueAsDouble(0,4,1);
768 
769  deValues.Add(param);
770  }
771 
772  AliMUONTrackerData* data = new AliMUONTrackerData(name,name,deValues,1);
773  data->SetDimensionName(0,name);
774 
775  return data;
776 }
777 
778 //______________________________________________________________________________
780 {
781  if (smode == "local") {
782  if ( input.EndsWith(".xml") ) return kInteractif_xml;
783  else if ( input.EndsWith(".txt") ) return kInteractif_ESDList;
784  else if ( input.EndsWith(".root") ) return kLocal;
785  } else if (smode == "caf" || smode == "saf" || smode == "saf3") return kProof;
786  return -1;
787 }
788 
789 //______________________________________________________________________________
790 TChain* CreateChainFromCollection(const char *xmlfile)
791 {
792  // Create a chain from the collection of tags.
793  if (!TGrid::Connect("alien://")) return NULL;
794 
795  TGridCollection* coll = TAlienCollection::Open(xmlfile);
796  if (!coll) {
797  Error("CreateChainFromCollection", "Cannot create the AliEn collection");
798  return NULL;
799  }
800 
801  TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE);
802  AliTagAnalysis *tagAna = new AliTagAnalysis("ESD");
803  tagAna->ChainGridTags(tagResult);
804 
805  AliRunTagCuts *runCuts = new AliRunTagCuts();
806  AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();
807  AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
808  AliEventTagCuts *evCuts = new AliEventTagCuts();
809 
810  // Check if the cuts configuration file was provided
811  if (!gSystem->AccessPathName("ConfigureCuts.C"))
812  gROOT->ProcessLine(Form(".x ConfigureCuts.C((AliRunTagCuts*)%p, (AliLHCTagCuts*)%p, (AliDetectorTagCuts*)%p,"
813  " (AliEventTagCuts*)%p)", runCuts, lhcCuts, detCuts, evCuts));
814 
815  TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts);
816  if (!chain || !chain->GetNtrees()) return NULL;
817  chain->ls();
818  return chain;
819 }
820 
821 //______________________________________________________________________________
822 TChain* CreateChainFromFile(const char *rootfile)
823 {
824  // Create a chain using the root file.
825  TChain* chain = new TChain("esdTree");
826  chain->Add(rootfile);
827  if (!chain->GetNtrees()) return NULL;
828  chain->ls();
829  return chain;
830 }
831 
832 //______________________________________________________________________________
833 TChain* CreateChainFromESDList(const char *esdList)
834 {
835  // Create a chain using tags from the run list.
836  TChain* chain = new TChain("esdTree");
837  ifstream inFile(esdList);
838  TString inFileName;
839  if (inFile.is_open()) {
840  while (! inFile.eof() ) {
841  inFileName.ReadLine(inFile,kFALSE);
842  if(!inFileName.EndsWith(".root")) continue;
843  chain->Add(inFileName.Data());
844  }
845  }
846  inFile.close();
847  if (!chain->GetNtrees()) return NULL;
848  chain->ls();
849  return chain;
850 }
851 
852 //______________________________________________________________________________
854 {
855  printf("*******************************\n");
856  printf("*** Getting the Chain ***\n");
857  printf("*******************************\n");
858  if(mode == kInteractif_xml) return CreateChainFromCollection(input.Data());
859  else if (mode == kInteractif_ESDList) return CreateChainFromESDList(input.Data());
860  else if (mode == kLocal) return CreateChainFromFile(input.Data());
861  else return NULL;
862 }
863 
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:40
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)