small fix to remove TF1Helper which sneaked in again, inadvertently...
[u/mrichter/AliRoot.git] / PWG3 / muondep / RunMuonResolution.C
1 //--------------------------------------------------------------------------
2 // Base macro for submitting muon Resolution analysis.
3 //
4 // In case it is not run with full aliroot, it needs the following libraries:
5 //  - libSTEERBase.so
6 //  - libESD.so
7 //  - libAOD.so
8 //  - libANALYSIS.so
9 //  - libANALYSISalice.so
10 //  - libGui.so
11 //  - libMinuit.so
12 //  - libProofPlayer.so
13 //  - libXMLParser.so
14 //  - libRAWDatabase.so
15 //  - libCDB.so
16 //  - libSTEER.so
17 //  - libMUONcore.so
18 //  - libMUONmapping.so
19 //  - libMUONcalib.so
20 //  - libMUONgeometry.so
21 //  - libMUONtrigger.so
22 //  - libMUONraw.so
23 //  - libMUONbase.so
24 //  - libMUONrec.so
25 //  - libCORRFW.so
26 //  - libPWG3muondep.so
27 //
28 // The macro reads ESDs and store outputs in standard output file (AnalysisResults.root)
29 // It also needs to load magnetic field, mapping, geometry (+alignment), and recoonstruction parameters from the OCDB
30 //
31 // Author: Philippe Pillot - SUBATECH Nantes
32 //--------------------------------------------------------------------------
33
34 TString aliroot="VO_ALICE@AliRoot::v4-19-15-AN";
35
36 enum {kLocal, kInteractif_xml, kInteractif_ESDList, kProof};
37
38 // global declaration needed in Proof mode, certainly because of the interpretor
39 TStopwatch* localTimer;
40 TMultiGraph* mgClusterResXVsStep;
41 TMultiGraph* mgClusterResYVsStep;
42 Double_t clusterResNB[10];
43 Double_t clusterResB[10];
44 Double_t clusterResNBErr[10];
45 Double_t clusterResBErr[10];
46
47 void RunMuonResolution(TString smode = "local", TString inputFileName = "AliESDs.root", Int_t nSteps = 10,
48                        Bool_t selectPhysics = kFALSE, Bool_t matchTrig = kFALSE, Double_t minMomentum = 0.,
49                        Bool_t correctForSystematics = kTRUE, Int_t extrapMode = 1)
50 {
51   /// Compute the cluster resolution by studying cluster-track residual, deconvoluting from track resolution
52   /// if extrapMode == 0: extrapolate from the closest cluster
53   /// if extrapMode == 1: extrapolate from the previous cluster except between stations 2-3-4
54   /// if correctForSystematics == kTRUE: the systematic shifts of the residuals is included in the resolution
55   
56   // timer start...
57   localTimer = new TStopwatch;
58   
59   // check parameters
60   nSteps = TMath::Max(nSteps,1);
61   if (extrapMode != 0 && extrapMode != 1) {
62     Error("RunMuonResolution","incorrect extrapolation mode!");
63     return;
64   }
65   
66   // Check runing mode
67   Int_t mode = GetMode(smode, inputFileName);
68   if(mode < 0){
69     Error("RunMuonResolution","Please provide either an ESD root file a collection of ESDs or a dataset.");
70     return;
71   }
72   
73   gSystem->Load("libVMC");
74   gSystem->Load("libSTEERBase");
75   gSystem->Load("libESD");
76   gSystem->Load("libAOD");
77   gSystem->Load("libANALYSIS");
78   gSystem->Load("libANALYSISalice");
79   
80   // Load additional libraries
81   gSystem->Load("libProofPlayer");
82   TString extraLibs="Physics:Minuit:XMLParser:Gui:RAWDatabase:CDB:STEER:MUONcore:MUONmapping:MUONcalib:MUONgeometry:MUONtrigger:MUONraw:MUONbase:MUONrec:CORRFW:PWG3muondep";
83   TObjArray* libs = extraLibs.Tokenize(":");
84   for (Int_t i = 0; i < libs->GetEntriesFast(); i++)
85     gSystem->Load(Form("lib%s",static_cast<TObjString*>(libs->UncheckedAt(i))->GetName()));
86   
87   // check for old output file to removed
88   char remove = '';
89   if (!gSystem->Exec("ls chamberResolution_step*[0-9].root")) {
90     cout<<"above files must be removed from the current directory. Delete? [y=yes, n=no] "<<flush;
91     while (remove != 'y' && remove != 'n') cin>>remove;
92     if (remove == 'y') gSystem->Exec("rm -f chamberResolution_step*[0-9].root");
93     else {
94       Error("RunMuonResolution","cannot proceed with these files there otherwise results will be mixed up!");
95       return;
96     }
97   }
98   
99   // OCDB access
100   if(mode == kLocal) AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
101   else if (mode != kProof) {
102     if (!TGrid::Connect("alien://")) return;
103     if(mode == kInteractif_ESDList || !gSystem->AccessPathName("ConfigureCuts.C"))
104       AliCDBManager::Instance()->SetDefaultStorage("raw://");
105   }
106   
107   // Create input object
108   TObject* inputObj = 0x0;
109   if (mode == kProof) inputObj = new TObjString(inputFileName);
110   else inputObj = CreateChain(mode, inputFileName);
111   if (!inputObj) return;
112   
113   // set starting chamber resolution (if -1 they will be loaded from recoParam in the task)
114   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
115     clusterResNB[i] = -1.;
116     clusterResB[i] = -1.;
117     clusterResNBErr[i] = 0.;
118     clusterResBErr[i] = 0.;
119   }
120   
121   // output graphs
122   mgClusterResXVsStep = new TMultiGraph("mgClusterResXVsStep","cluster X-resolution versus step;step;#sigma_{X} (cm)");
123   mgClusterResYVsStep = new TMultiGraph("mgClusterResYVsStep","cluster Y-resolution versus step;step;#sigma_{Y} (cm)");
124   TGraphErrors* gClusterResXVsStep[10];
125   TGraphErrors* gClusterResYVsStep[10];
126   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
127     gClusterResXVsStep[i] = new TGraphErrors(nSteps+1);
128     gClusterResXVsStep[i]->SetName(Form("gResX_ch%d",i+1));
129     gClusterResXVsStep[i]->SetMarkerStyle(kFullDotMedium);
130     gClusterResXVsStep[i]->SetMarkerColor(i+1+i/9);
131     mgClusterResXVsStep->Add(gClusterResXVsStep[i],"lp");
132     
133     gClusterResYVsStep[i] = new TGraphErrors(nSteps+1);
134     gClusterResYVsStep[i]->SetName(Form("gResY_ch%d",i+1));
135     gClusterResYVsStep[i]->SetMarkerStyle(kFullDotMedium);
136     gClusterResYVsStep[i]->SetMarkerColor(i+1+i/9);
137     mgClusterResYVsStep->Add(gClusterResYVsStep[i],"lp");
138   }
139   
140   // loop over step
141   for (Int_t iStep = 0; iStep < nSteps; iStep++) {
142     cout<<"step "<<iStep+1<<"/"<<nSteps<<endl;
143     
144     // Connect to proof if needed and prepare environment
145     if (mode == kProof) {
146       
147       // set general environment and close previous session
148       if (iStep == 0) {
149         gEnv->SetValue("XSec.GSI.DelegProxy","2");
150         TProof::AddEnvVar("ALIROOT_EXTRA_LIBS", extraLibs.Data());
151       } else gProof->Close("s");
152       
153       // connect
154       if (gSystem->Getenv("alien_API_USER") == NULL) TProof::Open("alice-caf.cern.ch","workers=20");
155       else TProof::Open(Form("%s@alice-caf.cern.ch",gSystem->Getenv("alien_API_USER")),"workers=20");
156       if (!gProof) return;
157       
158       // set environment and compile task on workers
159       gProof->EnablePackage(aliroot.Data());
160       
161       // prepare OCDB access on workers
162       gProof->Exec("AliCDBManager::Instance()->SetDefaultStorage(\"raw://\")");
163       
164     }
165     
166     // run one step
167     run(mode, iStep, inputObj, extrapMode, correctForSystematics, minMomentum, matchTrig, selectPhysics, clusterResNB, clusterResB);
168     
169     // fill graph with starting resolutions from the task at first step
170     if (iStep == 0) for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
171       gClusterResXVsStep[i]->SetPoint(0, 0, clusterResNB[i]);
172       gClusterResXVsStep[i]->SetPointError(0, 0., clusterResNBErr[i]);
173       gClusterResYVsStep[i]->SetPoint(0, 0, clusterResB[i]);
174       gClusterResYVsStep[i]->SetPointError(0, 0., clusterResBErr[i]);
175     }
176     
177     // read the chamber resolution from the output file
178     if (!GetChamberResolution(iStep, clusterResNB, clusterResB, clusterResNBErr, clusterResBErr)) return;
179     
180     // fill graphs with computed resolutions
181     for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
182       gClusterResXVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResNB[i]);
183       gClusterResXVsStep[i]->SetPointError(iStep+1, 0., clusterResNBErr[i]);
184       gClusterResYVsStep[i]->SetPoint(iStep+1, iStep+1, clusterResB[i]);
185       gClusterResYVsStep[i]->SetPointError(iStep+1, 0., clusterResBErr[i]);
186     }
187     
188   }
189   
190   // copy final results in results.root file
191   gSystem->Exec(Form("cp chamberResolution_step%d.root results.root", nSteps-1));
192   
193   // display convergence
194   TCanvas* convergence = new TCanvas("convergence","convergence");
195   convergence->Divide(1,2);
196   convergence->cd(1);
197   mgClusterResXVsStep->Draw("ap");
198   convergence->cd(2);
199   mgClusterResYVsStep->Draw("ap");
200   
201   // save convergence plots
202   TFile* outFile = TFile::Open("results.root","UPDATE");
203   if (!outFile || !outFile->IsOpen()) return;
204   outFile->cd();
205   mgClusterResXVsStep->Write();
206   mgClusterResYVsStep->Write();
207   convergence->Write();
208   outFile->Close();
209   
210   // print results
211   printf("\nchamber resolution:\n");
212   printf(" - non-bending:");
213   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",clusterResNB[i]);
214   printf("\n -     bending:");
215   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",clusterResB[i]);
216   printf("\n\n");
217   
218   // ...timer stop
219   localTimer->Stop();
220   localTimer->Print();
221   
222 }
223
224 //______________________________________________________________________________
225 void run(Int_t mode, Int_t iStep, TObject* input, Int_t extrapMode, Bool_t correctForSystematics,
226          Double_t minMomentum, Bool_t matchTrig, Bool_t selectPhysics, Double_t clusterResNB[10], Double_t clusterResB[10])
227 {
228   /// launch the analysis with these parameters
229   
230   // Create the analysis manager
231   AliAnalysisManager *mgr = new AliAnalysisManager("MuonResolutionAnalysis");
232   
233   // ESD input handler
234   AliESDInputHandler* esdH = new AliESDInputHandler();
235   esdH->SetReadFriends(kFALSE);
236   mgr->SetInputEventHandler(esdH);
237   
238   // event selection
239   if (selectPhysics) {
240     gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
241     AliPhysicsSelectionTask* physicsSelection = AddTaskPhysicsSelection();
242     if(!physicsSelection) {
243       Error("RunMuonResolution","AliPhysicsSelectionTask not created!");
244       return;
245     }
246     physicsSelection->GetPhysicsSelection()->SetUseMuonTriggers();
247   }
248   
249   // Muon Resolution analysis
250   gROOT->LoadMacro("$ALICE_ROOT/PWG3/muondep/AddTaskMuonResolution.C");
251   AliAnalysisManager::SetCommonFileName(Form("chamberResolution_step%d.root", iStep));
252   AliAnalysisTaskMuonResolution* muonResolution = AddTaskMuonResolution(selectPhysics, matchTrig, minMomentum, correctForSystematics, extrapMode);
253   if(!muonResolution) {
254     Error("RunMuonResolution","AliAnalysisTaskMuonResolution not created!");
255     return;
256   }
257   muonResolution->SetStartingResolution(clusterResNB, clusterResB);
258   
259   // Enable debug printouts
260   //mgr->SetDebugLevel(2);
261   
262   // start local analysis
263   if (mgr->InitAnalysis()) {
264     mgr->PrintStatus();
265     if (mode == kProof) mgr->StartAnalysis("proof", Form("%s",static_cast<TObjString*>(input)->GetName()));
266     else mgr->StartAnalysis("local", static_cast<TChain*>(input));
267   }
268   
269   // save the summary canvases
270   if (muonResolution->GetCanvases()) {
271     TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"UPDATE");
272     if (outFile && outFile->IsOpen()) {
273       muonResolution->GetCanvases()->Write();
274       outFile->Close();
275     }
276   }
277   
278   // return starting chamber resolution from the task
279   muonResolution->GetStartingResolution(clusterResNB, clusterResB);
280   
281   // clean memory
282   delete mgr;
283   TObject::SetObjectStat(kFALSE);
284 }
285
286 //______________________________________________________________________________
287 Bool_t GetChamberResolution(Int_t iStep, Double_t clusterResNB[10], Double_t clusterResB[10], Double_t clusterResNBErr[10], Double_t clusterResBErr[10])
288 {
289   /// read the chamber resolution from the output file
290   
291   TFile* outFile = TFile::Open(Form("chamberResolution_step%d.root", iStep),"READ");
292   
293   if (!outFile || !outFile->IsOpen()) {
294     Error("GetChamberResolution","output file does not exist!");
295     return kFALSE;
296   }
297   
298   TObjArray* summary = static_cast<TObjArray*>(outFile->FindObjectAny("ChamberRes"));
299   TGraphErrors* gCombinedResidualXPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualXPerChSigma")) : 0x0;
300   TGraphErrors* gCombinedResidualYPerChSigma = (summary) ? static_cast<TGraphErrors*>(summary->FindObject("gCombinedResidualYPerChSigma")) : 0x0;
301   
302   if (!gCombinedResidualXPerChSigma || !gCombinedResidualYPerChSigma) {
303     Error("GetChamberResolution","resolution graphs do not exist!");
304     return kFALSE;
305   }
306   
307   Double_t dummy;
308   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
309     gCombinedResidualXPerChSigma->GetPoint(i, dummy, clusterResNB[i]);
310     gCombinedResidualYPerChSigma->GetPoint(i, dummy, clusterResB[i]);
311     clusterResNBErr[i] = gCombinedResidualXPerChSigma->GetErrorY(i);
312     clusterResBErr[i] = gCombinedResidualYPerChSigma->GetErrorY(i);
313   }
314   
315   outFile->Close();
316   
317   return kTRUE;
318 }
319
320 //______________________________________________________________________________
321 Int_t GetMode(TString smode, TString input)
322 {
323   if (smode == "local") {
324     if ( input.EndsWith(".xml") ) return kInteractif_xml;
325     else if ( input.EndsWith(".txt") ) return kInteractif_ESDList;
326     else if ( input.EndsWith(".root") ) return kLocal;    
327   } else if (smode == "proof") return kProof;
328   return -1;
329 }
330
331 //______________________________________________________________________________
332 TChain* CreateChainFromCollection(const char *xmlfile)
333 {
334   // Create a chain from the collection of tags.
335   TAlienCollection* coll = TAlienCollection::Open(xmlfile);
336   if (!coll) {
337     Error("CreateChainFromCollection", Form("Cannot create an AliEn collection from %s!", xmlfile));
338     return NULL;
339   }
340   
341   TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE);
342   AliTagAnalysis *tagAna = new AliTagAnalysis("ESD");
343   tagAna->ChainGridTags(tagResult);
344   
345   AliRunTagCuts      *runCuts = new AliRunTagCuts();
346   AliLHCTagCuts      *lhcCuts = new AliLHCTagCuts();
347   AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
348   AliEventTagCuts    *evCuts  = new AliEventTagCuts();
349   
350   // Check if the cuts configuration file was provided
351   if (!gSystem->AccessPathName("ConfigureCuts.C")) {
352     gROOT->LoadMacro("ConfigureCuts.C");
353     ConfigureCuts(runCuts, lhcCuts, detCuts, evCuts);
354   }
355   
356   TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts);
357   if (!chain || !chain->GetNtrees()) return NULL;
358   chain->ls();
359   return chain;
360 }
361
362 //______________________________________________________________________________
363 TChain* CreateChainFromFile(const char *rootfile)
364 {
365   // Create a chain using the root file.
366   TChain* chain = new TChain("esdTree");
367   chain->Add(rootfile);
368   if (!chain->GetNtrees()) return NULL;
369   chain->ls();
370   return chain;
371 }
372
373 //______________________________________________________________________________
374 TChain* CreateChainFromESDList(const char *esdList)
375 {
376   // Create a chain using tags from the run list.
377   TChain* chain = new TChain("esdTree");
378   ifstream inFile(esdList);
379   TString inFileName;
380   if (inFile.is_open()) {
381     while (! inFile.eof() ) {
382       inFileName.ReadLine(inFile,kFALSE);
383       if(!inFileName.EndsWith(".root")) continue;
384       chain->Add(inFileName.Data());
385     }
386   }
387   inFile.close();
388   if (!chain->GetNtrees()) return NULL;
389   chain->ls();
390   return chain;
391 }
392
393 //______________________________________________________________________________
394 TChain* CreateChain(Int_t mode, TString input)
395 {
396   printf("*******************************\n");
397   printf("*** Getting the Chain       ***\n");
398   printf("*******************************\n");
399   if(mode == kInteractif_xml) return CreateChainFromCollection(input.Data());
400   else if (mode == kInteractif_ESDList) return CreateChainFromESDList(input.Data());
401   else if (mode == kLocal) return CreateChainFromFile(input.Data());
402   else return NULL;
403 }
404