major dielectron update (included also the data and plotting macros for paper)
[u/mrichter/AliRoot.git] / PWG3 / dielectron / macros / FitCDFLocal.C
1 TH1F *psproperMCsecJPSI();
2 void LoadLib();
3
4 void FitCDFLocal(){
5
6   ///////////////////////////////////////////////////////////////////
7   //
8   // Example macro to read local N-Tuples of JPSI 
9   // and bkg candidates and perform log-likelihood 
10   // minimization using the minimization handler class
11   //
12   // Origin: C. Di Giglio
13   //
14   ///////////////////////////////////////////////////////////////////
15   LoadLib();
16   AliDielectronBtoJPSItoEleCDFfitFCN *likely_obj = 0x0;
17   
18   // Those initial parameters are evaluated from the fits done separately on
19   // invariant mass and pseudoproper distributions. For the moment are setted by
20   // hand 
21
22   Double_t *resParam = new Double_t[6];
23   resParam[0] = 2586.32; resParam[1]= -0.163616; resParam[2] = 31.026770; resParam[3] = 433.85; 
24   resParam[4] = 4.67; resParam[5] = 73.31;
25
26   Double_t *bkgParam  = new Double_t[7];
27   bkgParam[0] = 7063.8; bkgParam[1] = 23214.45; bkgParam[2] = 24742.65; bkgParam[3] = 15188.57; 
28   bkgParam[4] = 0.009579; bkgParam[5] = 0.010422; bkgParam[6] = 0.000423;
29
30   Double_t *invMassParam = new Double_t[9];
31   invMassParam[0] = 3.092; invMassParam[1] =  0.0306; invMassParam[2] =  0.4939; 
32   invMassParam[3] = 1.36; invMassParam[4] = 41.97; invMassParam[5] = 80.43; invMassParam[6] = 0.975;
33   invMassParam[7] = 1.61; invMassParam[8] = 2.129;
34   
35   Double_t Fsig = 0.454;
36
37   TH1F* hCsiMCPithya = new TH1F();
38
39   hCsiMCPithya = psproperMCsecJPSI();  
40   Double_t integral = 0;
41   for(int i=1;i<hCsiMCPithya->GetNbinsX()+1; i++) integral += (hCsiMCPithya->GetBinContent(i)*hCsiMCPithya->GetBinWidth(i));
42   hCsiMCPithya->Scale(1./integral);
43
44   Double_t* x=0x0; Double_t* m=0x0; Int_t n=0;
45   AliDielectronBtoJPSItoEle* aBtoJPSItoEle =new AliDielectronBtoJPSItoEle();
46   Double_t paramInputValues [20] =   
47                                    /*paramInputValues[0]   ----> fWeightRes;
48                                      paramInputValues[1]   ----> fPos;
49                                      paramInputValues[2]   ----> fNeg;
50                                      paramInputValues[3]   ----> fSym;
51                                      paramInputValues[4]   ----> fOneOvLamPlus;
52                                      paramInputValues[5]   ----> fOneOvLamMinus;
53                                      paramInputValues[6]   ----> fOneOvLamSym;
54                                      paramInputValues[7]   ----> fB;
55                                      paramInputValues[8]   ----> fFsig;
56                                      paramInputValues[9]   ----> fMmean;
57                                      paramInputValues[10]  ----> fNexp;
58                                      paramInputValues[11]  ----> fSigma;
59                                      paramInputValues[12]  ----> fAlpha;
60                                      paramInputValues[13]  ----> fNorm;
61                                      paramInputValues[14]  ----> fBkgNorm;
62                                      paramInputValues[15]  ----> fBkgMean;
63                                      paramInputValues[16]  ----> fBkgSlope;
64                                      paramInputValues[17]  ----> fBkgConst;
65                                      paramInputValues[18]  ----> fGaus1Norm; 
66                                      paramInputValues[19]  ----> fGaus2Norm;*/                                     
67                                    { bkgParam[9],  // fWeightRes [0]
68                                      bkgParam[10], // Fpos [1]
69                                      bkgParam[11], // FNeg [2]
70                                      bkgParam[12], // FSym [3]
71                                      bkgParam[6], // LamdaPos [4]
72                                      bkgParam[7], // LambdaNeg [5]
73                                      bkgParam[8], // LamdaSym [6]
74                                      0.10, // fB [7]
75                                      Fsig, // FSig [8]
76                                      invMassParam[0], // FMean [9]
77                                      invMassParam[4], // fNexp [10]
78                                      invMassParam[1], // fNsigma [11]
79                                      invMassParam[2], // fAlpha [12]
80                                      invMassParam[3], //fNorm [13]
81                                      invMassParam[5], // fBkgNorm [14]
82                                      invMassParam[6], // fBkgMean [15]
83                                      invMassParam[7], // fBkgSlope [16]
84                                      invMassParam[8], // fBkgConst [17]
85                                      resParam[0], // Gaus1Norm [18]
86                                      resParam[3] // Gaus2Norm [19]
87                                       }; // Starting values for parameters 
88
89   // retrieve the TNtupla from the file result.root
90   TFile f("result.root");
91   TList * list = (TList*)f.Get("resultAOD");
92   TNtuple *nt=(TNtuple*)list->FindObject("fNtupleJPSI");
93   nt->ls();
94   aBtoJPSItoEle->ReadCandidates(nt,x,m,n); // read N-Tuples
95   printf("+++\n+++ Number of total candidates (prim J/psi + secondary J/psi + bkg) ---> %d candidates \n+++\n",n);
96
97   aBtoJPSItoEle->SetFitHandler(x,m,n); // Set the fit handler with given values of x, m, # of candidates 
98
99   aBtoJPSItoEle->CloneMCtemplate(hCsiMCPithya);    // clone MC template and copy internally
100                                                    // in this way any model can be setted from outside
101
102   aBtoJPSItoEle->SetCsiMC(); // Pass the MC template to the CDF fit function
103
104   AliDielectronBtoJPSItoEleCDFfitHandler* fithandler = aBtoJPSItoEle->GetCDFFitHandler(); // Get the fit handler
105
106   //
107   // Set some fit options through the handler class
108   //
109   fithandler->SetErrorDef(0.5); // tells Minuit that the error interval is the one in which
110                                 // the function differs from the minimum for less than setted value
111
112   fithandler->SetPrintStatus(kTRUE);
113
114   fithandler->SetParamStartValues(paramInputValues);
115
116   Double_t *resConst = new Double_t[4];
117   
118   //resolution constants: fixed from MC
119   resConst[0]= resParam[1]; // mean1
120   resConst[1]= resParam[2]; // sigma1
121   resConst[2]= resParam[4]; // mean2
122   resConst[3]= resParam[5]; // sigma2
123
124   
125   fithandler->SetResolutionConstants(resConst); 
126   fithandler->SetCrystalBallFunction(kTRUE);
127   
128   // fix parameter
129   fithandler->FixParam(0,kTRUE);
130   fithandler->FixParam(1,kTRUE); // Fix Bkg  weights  
131   fithandler->FixParam(2,kTRUE);
132   fithandler->FixParam(3,kTRUE);
133
134   fithandler->FixParam(4,kTRUE);
135   fithandler->FixParam(5,kTRUE); // Fix bkg exponential 
136   fithandler->FixParam(6,kTRUE);
137   //fithandler->FixParam(8,kTRUE); // Fsig fix
138   
139   fithandler->FixParam(9,kTRUE);
140   fithandler->FixParam(10,kTRUE); // Fix CristalBall param
141   fithandler->FixParam(11,kTRUE);
142   fithandler->FixParam(12,kTRUE);  
143   fithandler->FixParam(13,kTRUE); 
144
145   fithandler->FixParam(14,kTRUE);
146   fithandler->FixParam(15,kTRUE); // Fix Bkg invMass param
147   fithandler->FixParam(16,kTRUE);
148   fithandler->FixParam(17,kTRUE);
149   
150   fithandler->FixParam(18,kTRUE); // Fix resolution weights
151   fithandler->FixParam(19,kTRUE);
152   
153   //invariant mass function normalized between invMassMin, invMassMax  -->
154   Double_t invMassMin = 2.0; // GeV
155   Double_t invMassMax = 4.0; // GeV 
156   Double_t massLow = (TDatabasePDG::Instance()->GetParticle(443)->Mass()) - invMassMin;
157   Double_t massHigh = invMassMax - (TDatabasePDG::Instance()->GetParticle(443)->Mass());
158   fithandler->SetMassWndLow(massLow);
159   fithandler->SetMassWndHigh(massHigh);
160
161   likely_obj = fithandler->LikelihoodPointer();
162   likely_obj->SetAllParameters(paramInputValues); 
163   likely_obj->ComputeMassIntegral();
164   likely_obj->PrintStatus(); 
165  
166   aBtoJPSItoEle->DoMinimization();
167   f.Close();
168   return;
169 }
170
171
172 TH1F *psproperMCsecJPSI(){
173    TH1F *hDecayTimeMCjpsifromB = new TH1F("hDecayTimeMCjpsifromB","B --> J/#Psi MC pseudo proper decay length",150,-3000,3000
174 ); 
175    hDecayTimeMCjpsifromB->SetBinContent(67,2);
176    hDecayTimeMCjpsifromB->SetBinContent(68,2);
177    hDecayTimeMCjpsifromB->SetBinContent(69,2);
178    hDecayTimeMCjpsifromB->SetBinContent(70,4);
179    hDecayTimeMCjpsifromB->SetBinContent(71,5);
180    hDecayTimeMCjpsifromB->SetBinContent(72,17);
181    hDecayTimeMCjpsifromB->SetBinContent(73,35);
182    hDecayTimeMCjpsifromB->SetBinContent(74,96);
183    hDecayTimeMCjpsifromB->SetBinContent(75,439);
184    hDecayTimeMCjpsifromB->SetBinContent(76,11428);
185    hDecayTimeMCjpsifromB->SetBinContent(77,9094);
186    hDecayTimeMCjpsifromB->SetBinContent(78,7769);
187    hDecayTimeMCjpsifromB->SetBinContent(79,6475);
188    hDecayTimeMCjpsifromB->SetBinContent(80,5695);
189    hDecayTimeMCjpsifromB->SetBinContent(81,4884);
190    hDecayTimeMCjpsifromB->SetBinContent(82,4294);
191    hDecayTimeMCjpsifromB->SetBinContent(83,3780);
192    hDecayTimeMCjpsifromB->SetBinContent(84,3321);
193    hDecayTimeMCjpsifromB->SetBinContent(85,2827);
194    hDecayTimeMCjpsifromB->SetBinContent(86,2531);
195    hDecayTimeMCjpsifromB->SetBinContent(87,2256);
196    hDecayTimeMCjpsifromB->SetBinContent(88,2099);
197    hDecayTimeMCjpsifromB->SetBinContent(89,1782);
198    hDecayTimeMCjpsifromB->SetBinContent(90,1592);
199    hDecayTimeMCjpsifromB->SetBinContent(91,1478);
200    hDecayTimeMCjpsifromB->SetBinContent(92,1286);
201    hDecayTimeMCjpsifromB->SetBinContent(93,1145);
202    hDecayTimeMCjpsifromB->SetBinContent(94,980);
203    hDecayTimeMCjpsifromB->SetBinContent(95,933);
204    hDecayTimeMCjpsifromB->SetBinContent(96,865);
205    hDecayTimeMCjpsifromB->SetBinContent(97,761);
206    hDecayTimeMCjpsifromB->SetBinContent(98,654);
207    hDecayTimeMCjpsifromB->SetBinContent(99,622);
208    hDecayTimeMCjpsifromB->SetBinContent(100,587);
209    hDecayTimeMCjpsifromB->SetBinContent(101,572);
210    hDecayTimeMCjpsifromB->SetBinContent(102,449);
211    hDecayTimeMCjpsifromB->SetBinContent(103,416);
212    hDecayTimeMCjpsifromB->SetBinContent(104,346);
213    hDecayTimeMCjpsifromB->SetBinContent(105,361);
214    hDecayTimeMCjpsifromB->SetBinContent(105,361);
215    hDecayTimeMCjpsifromB->SetBinContent(106,315);
216    hDecayTimeMCjpsifromB->SetBinContent(107,265);
217    hDecayTimeMCjpsifromB->SetBinContent(108,239);
218    hDecayTimeMCjpsifromB->SetBinContent(109,247);
219    hDecayTimeMCjpsifromB->SetBinContent(110,189);
220    hDecayTimeMCjpsifromB->SetBinContent(111,223);
221    hDecayTimeMCjpsifromB->SetBinContent(112,174);
222    hDecayTimeMCjpsifromB->SetBinContent(113,184);
223    hDecayTimeMCjpsifromB->SetBinContent(114,171);
224    hDecayTimeMCjpsifromB->SetBinContent(115,153);
225    hDecayTimeMCjpsifromB->SetBinContent(116,151);
226    hDecayTimeMCjpsifromB->SetBinContent(117,126);
227    hDecayTimeMCjpsifromB->SetBinContent(118,105);
228    hDecayTimeMCjpsifromB->SetBinContent(119,98);
229    hDecayTimeMCjpsifromB->SetBinContent(120,88);
230    hDecayTimeMCjpsifromB->SetBinContent(121,92);
231    hDecayTimeMCjpsifromB->SetBinContent(122,86);
232    hDecayTimeMCjpsifromB->SetBinContent(123,73);
233    hDecayTimeMCjpsifromB->SetBinContent(124,78);
234    hDecayTimeMCjpsifromB->SetBinContent(125,63);
235    hDecayTimeMCjpsifromB->SetBinContent(126,69);
236    hDecayTimeMCjpsifromB->SetBinContent(127,55);
237    hDecayTimeMCjpsifromB->SetBinContent(128,43);
238    hDecayTimeMCjpsifromB->SetBinContent(129,46);
239    hDecayTimeMCjpsifromB->SetBinContent(130,38);
240    hDecayTimeMCjpsifromB->SetBinContent(131,34);
241    hDecayTimeMCjpsifromB->SetBinContent(132,33);
242    hDecayTimeMCjpsifromB->SetBinContent(133,51);
243    hDecayTimeMCjpsifromB->SetBinContent(134,25);
244    hDecayTimeMCjpsifromB->SetBinContent(135,44);
245    hDecayTimeMCjpsifromB->SetBinContent(136,24);
246    hDecayTimeMCjpsifromB->SetBinContent(137,25);
247    hDecayTimeMCjpsifromB->SetBinContent(138,21);
248    hDecayTimeMCjpsifromB->SetBinContent(139,23);
249    hDecayTimeMCjpsifromB->SetBinContent(140,16);
250    hDecayTimeMCjpsifromB->SetBinContent(141,24);
251    hDecayTimeMCjpsifromB->SetBinContent(142,14);
252    hDecayTimeMCjpsifromB->SetBinContent(143,21);
253    hDecayTimeMCjpsifromB->SetBinContent(144,22);
254    hDecayTimeMCjpsifromB->SetBinContent(145,22);
255    hDecayTimeMCjpsifromB->SetBinContent(146,13);
256    hDecayTimeMCjpsifromB->SetBinContent(147,10);
257    hDecayTimeMCjpsifromB->SetBinContent(148,15);
258    hDecayTimeMCjpsifromB->SetBinContent(149,12);
259    hDecayTimeMCjpsifromB->SetBinContent(150,12);
260    hDecayTimeMCjpsifromB->SetBinContent(151,231);
261    return hDecayTimeMCjpsifromB;
262 }
263
264 void LoadLib(){
265   gSystem->Load("libTree.so");
266   gSystem->Load("libGeom.so");
267   gSystem->Load("libPhysics.so");
268   gSystem->Load("libVMC.so");
269   gSystem->Load("libMinuit.so");
270   gSystem->Load("libSTEERBase.so");
271   gSystem->Load("libESD.so");
272   gSystem->Load("libAOD.so");
273   gSystem->Load("libANALYSIS.so");
274   gSystem->Load("libANALYSISalice.so");
275   gSystem->Load("libCORRFW.so");
276   gSystem->Load("libPWG3dielectron.so");
277  }