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