Added some more scripts
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / MakeELossFit.C
1 #ifndef __CINT__
2 #include "AliForwardCorrectionManager.h"
3 #include "AliFMDCorrELossFit.h"
4 #include <TAxis.h>
5 #include <TFile.h>
6 #include <TList.h>
7 #include <TError.h>
8 #include <TH1.h>
9 #include <TF1.h>
10 #include <TClonesArray.h>
11 #else
12 class TAxis;
13 class AliFMDCirrELossFit;
14 class TH1;
15
16 #endif
17
18 /**
19  * Class to make correction object and save to file 
20  *
21  * Run like                                     
22  *
23  * @verbatim 
24  * gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
25  * gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/Compile.C");
26  * Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/MakeELossFit.C"); 
27  * MakeELossFit mef(sys, cms, field, mc, "AnalysisResults.root"); 
28  * mef.Run();
29  * @endverbatim 
30  * where @e sys, the collision system, is one of 
31  * - 1: pp
32  * - 2: PbPb
33  * 
34  * @e cms is the center of mass energy per nuclean in GeV (e.g., 2760
35  * for first PbPb data), @e field is (signed) L3 magnetic in kG, and 
36  * @e mc is wether this correction applies to MC data or not. 
37  * 
38  * The class generates a file like 
39  * @verbatim
40  * elossfits<sys>_<cms>GeV_<field>kG_<realmc>.root
41  * @endverbatim 
42  * in the working directory. This file can be moved to 
43  * @verbatim
44  * $(ALICE_ROOT)/PWG2/FORWARD/corrections/ELossFits
45  * @endverbatim 
46  * and stored in SVN for later use. 
47  *
48  * @ingroup pwg2_forward_analysis_scripts
49  */
50
51 class MakeELossFit 
52 {
53 protected: 
54 public:
55   TList* fFitter;
56   TAxis* fAxis;
57   TClonesArray fFits;
58   UShort_t     fSys;
59   UShort_t     fCMS;
60   Short_t      fField;
61   Bool_t       fMC;
62
63   //__________________________________________________________________
64   MakeELossFit(UShort_t    sys, 
65                UShort_t    cms, 
66                Short_t     field, 
67                Bool_t      mc, 
68                const char* filename="AnalysisResults.root") 
69     : fFitter(0), 
70       fAxis(0),
71       fFits("AliFMDCorrELossFit::ELossFit"),
72       fSys(sys), 
73       fCMS(cms), 
74       fField(field), 
75       fMC(mc)
76   {
77     TFile* file = TFile::Open(filename, "READ");
78     if (!file) { 
79       Error("MakeELossFit", "Failed to open file %s", filename);
80       return;
81     }
82     TList* forward = static_cast<TList*>(file->Get("Forward"));
83     // static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
84     if (!forward) { 
85       Error("MakeELossFit", "Couldn't get forward list from %s", filename);
86       return;
87     }
88
89     fFitter = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
90     if (!fFitter) { 
91       Error("MakeELossFit", "Couldn't get fitter folder");
92       return;
93     }
94
95     fAxis = static_cast<TAxis*>(fFitter->FindObject("etaAxis"));
96     if (!fAxis) { 
97       Error("MakeELossFit", "Couldn't get eta axis");
98       return;
99     }
100     file->Close();
101
102 #if 0
103     AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
104     mgr.Init(sys, cms, field, mc, 0, true);
105 #endif
106   }
107   //__________________________________________________________________
108   TList* GetRing(UShort_t d, Char_t r) const
109   {
110     TList* rL = static_cast<TList*>(fFitter->FindObject(Form("FMD%d%c",d,r)));
111     if (!rL) { 
112       Warning("DrawFits", "List FMD%d%c not found", d, r);
113       return 0;
114     }
115     return rL;
116   }
117   //__________________________________________________________________
118   TList* GetEDists(UShort_t d, Char_t r) const
119   {
120     TList* rList = GetRing(d, r);
121     if (!rList) 
122       return 0;
123     // rList->ls();
124
125     TList* edists = static_cast<TList*>(rList->FindObject("EDists"));
126     if (!edists) { 
127       Error("PrintFits", "Couldn't get EDists list for FMD%d%c", d,r);
128       return 0; 
129     }
130     return edists;
131   }
132   //__________________________________________________________________
133   TH1* GetEDist(UShort_t d, Char_t r, UShort_t etabin)
134   {
135     TList* eList = GetEDists(d, r);
136     if (!eList) { 
137       Warning("GetEDist", "No list for FMD%d%c", d, r);
138       return 0;
139     }
140     // eList->ls();
141
142     TString cmp(Form("FMD%d%c_etabin%03d", d, r, etabin));
143
144     TIter next(eList);
145     TObject* o = 0;
146     while ((o = next())) { 
147       if (!cmp.CompareTo(o->GetName())) {
148         return static_cast<TH1*>(o);
149       }
150     }
151     return 0;
152   }
153 #ifndef __CINT__
154   //__________________________________________________________________
155   AliFMDCorrELossFit::ELossFit* FindBestFit(TH1* dist) 
156   {
157     TList* funcs = dist->GetListOfFunctions();
158     TIter  next(funcs);
159     TF1*   func = 0;
160     fFits.Clear();
161     Int_t  i = 0;
162     Info("FindBestFit", "%s", dist->GetName());
163     while ((func = static_cast<TF1*>(next()))) { 
164       AliFMDCorrELossFit::ELossFit* fit = 
165         new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
166       fit->CalculateQuality(10, .20, 1e-7);
167     }
168     fFits.Sort(false);
169     // fFits.Print();
170     return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
171   }
172 #endif
173   //__________________________________________________________________
174   Bool_t Run()
175   {
176     if (!fFitter || !fAxis) { 
177       Error("Run", "Missing objects");
178       return kFALSE;
179     }
180     AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
181     obj->SetEtaAxis(*fAxis);
182
183     for (UShort_t d = 1; d <= 3; d++) { 
184       Info("Run", "detector is FMD%d", d);
185       UShort_t nQ = (d == 1 ? 1 : 2);
186       for (UShort_t q = 0; q < nQ; q++) { 
187         Char_t r = (q == 0 ? 'I' : 'O');
188         Int_t nBin = fAxis->GetNbins();
189         Info("Run", " ring is FMD%d%c - %d bins", d, r, nBin);
190         
191         Bool_t oneSeen = kFALSE;
192         for (UShort_t b = 1; b <= nBin; b++) { 
193           TH1* h = GetEDist(d, r, b);
194           if (oneSeen && !h) break;
195           if (!h)            continue;
196           if (!oneSeen)      oneSeen = true;
197
198           AliFMDCorrELossFit::ELossFit* best = FindBestFit(h);
199           best->Print("");
200           best->fDet  = d; 
201           best->fRing = r;
202           best->fBin  = b;
203           // Double_t eta = fAxis->GetBinCenter(b);
204           Info("Run", "    bin=%3d->%6.4f", b, eta);
205           obj->SetFit(d, r, b, new AliFMDCorrELossFit::ELossFit(*best));
206         }
207       }
208     }
209     AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
210     TString fname(mgr.GetFileName(AliForwardCorrectionManager::kELossFits,
211                                   fSys, fCMS, fField, fMC));
212     TFile* output = TFile::Open(fname.Data(), "RECREATE");
213     if (!output) { 
214       Warning("Run", "Failed to open output file %s", fname.Data());
215       return kFALSE;
216     }
217     obj->Write(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
218     output->Write();
219     output->Close();
220     Info("Run", "File %s created.  It should be copied to %s and stored in SVN",
221          fname.Data(),mgr.GetFileDir(AliForwardCorrectionManager::kELossFits));
222     
223     return kTRUE;
224   }
225 };
226 //
227 // EOF
228 //