]>
Commit | Line | Data |
---|---|---|
57a6d361 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /*********************************************** | |
17 | ||
18 | Lambda Analysis Module | |
19 | ---------------------- | |
20 | ||
21 | This Class enables the analysis of output files | |
22 | created with AliAnalysisTaskExtractV0 and | |
23 | AliAnalysisTaskExtractPerformanceV0 grid tasks. | |
24 | ||
25 | It constructs corrected Lambda, AntiLambda and | |
26 | K0Short spectra. | |
27 | ||
28 | This version: 27th April 2012 | |
29 | ||
30 | To be compiled as a library | |
31 | (.L AliV0Module.cxx++ or something to that effect) | |
32 | ||
33 | --- David Dobrigkeit Chinellato | |
34 | daviddc@ifi.unicamp.br | |
35 | ||
36 | ***********************************************/ | |
37 | ||
38 | //--- For C++ ---- | |
39 | #include <TApplication.h> | |
d1d5d7e6 | 40 | #include <TMatrix.h> |
41 | #include <TMatrixD.h> | |
57a6d361 | 42 | #include <TROOT.h> |
43 | #include <TMath.h> | |
44 | #include <TF1.h> | |
45 | #include <TH1F.h> | |
46 | #include <TH2F.h> | |
47 | #include <TH3F.h> | |
48 | #include <TList.h> | |
49 | #include <TLine.h> | |
50 | #include <TFile.h> | |
51 | #include <TTree.h> | |
52 | #include <TCanvas.h> | |
53 | #include <TString.h> | |
54 | #include <THashList.h> | |
55 | #include <TDirectoryFile.h> | |
56 | #include <cstdlib> | |
57 | using namespace std; | |
58 | ||
59 | //--- For ROOT --- | |
60 | #include "AliV0Module.h" | |
61 | ||
62 | AliV0Module::AliV0Module() | |
63 | { | |
64 | //Dummy Constructor that sets some default values | |
65 | //so that nothing too horrible will happen (hopefully) | |
66 | fWhichParticle = "Lambda"; //Default | |
67 | fCINT1BoverINELratio = 1; | |
68 | fRapidityBoundary = 0.5; | |
69 | fRealDataFile = ""; | |
70 | fMCDataFile = ""; | |
71 | fFeedDownDataFile = ""; | |
72 | fOutputDataFile = ""; | |
73 | ||
74 | //Selections | |
75 | fCutV0Radius = -1; | |
76 | fCutDCANegToPV = -1; | |
77 | fCutDCAPosToPV = -1; | |
78 | fCutDCAV0Daughters = 1000; | |
79 | fCutV0CosPA = -2; | |
80 | fCutProperLifetime = 1e+6; | |
81 | fCutTPCPIDNSigmas = 1e+6; | |
82 | fCutNSigmasForSignalExtraction = 5; | |
83 | fCutLeastNumberOfCrossedRows = 70; | |
84 | fCutDaughterEta = 0.8; | |
85 | fCutCompetingV0Rejection = -1; | |
86 | ||
87 | //Set Feeddown Treatment | |
88 | fFDSwitch = "UseMCRatio"; | |
89 | ||
90 | //Default: Use bin counting | |
91 | fFitBackgroundSwitch = kFALSE; | |
92 | ||
d1d5d7e6 | 93 | //Default: Min-Bias |
94 | fPerformMultiplicityStudy = kFALSE; | |
95 | fLoMultBound = -1; | |
96 | fHiMultBound = 10000; | |
97 | ||
98 | //Default: no cut in Armenteros | |
99 | fSpecialArmenterosCutK0s = kFALSE; | |
100 | ||
57a6d361 | 101 | //Pt Bins: undefined |
102 | fptbinnumb = -1; | |
103 | } | |
104 | ||
105 | AliV0Module::AliV0Module(TString fParticleType) | |
106 | { | |
107 | // Allows definition of Particle Type in analysis. | |
108 | // Possible Options are "Lambda", "AntiLambda" and "K0Short". | |
109 | // If some other string is given, this constructor will | |
110 | // default to "Lambda". | |
111 | fWhichParticle = fParticleType; | |
112 | if(fWhichParticle!="Lambda"&&fWhichParticle!="AntiLambda"&&fWhichParticle!="K0Short"){ | |
113 | cout<<"Particle Type "<<fParticleType<<" unknown. Set to lambda."<<endl; | |
114 | fWhichParticle = "Lambda"; | |
115 | } | |
116 | fCINT1BoverINELratio = 1; | |
117 | fRapidityBoundary = 0.5; | |
118 | fRealDataFile = ""; | |
119 | fMCDataFile = ""; | |
120 | fFeedDownDataFile = ""; | |
121 | fOutputDataFile = ""; | |
122 | ||
123 | //Selections | |
124 | fCutV0Radius = -1; | |
125 | fCutDCANegToPV = -1; | |
126 | fCutDCAPosToPV = -1; | |
127 | fCutDCAV0Daughters = 1000; | |
128 | fCutV0CosPA = -2; | |
129 | fCutProperLifetime = 1e+6; | |
130 | fCutTPCPIDNSigmas = 1e+6; | |
131 | fCutNSigmasForSignalExtraction = 5; | |
132 | fCutLeastNumberOfCrossedRows = 70; | |
133 | fCutDaughterEta = 0.8; | |
134 | fCutCompetingV0Rejection = -1; | |
135 | ||
136 | //Set Feeddown Treatment | |
137 | fFDSwitch = "UseMCRatio"; | |
138 | ||
139 | //Default: Use bin counting | |
140 | fFitBackgroundSwitch = kFALSE; | |
141 | ||
d1d5d7e6 | 142 | //Default: Min-Bias |
143 | fPerformMultiplicityStudy = kFALSE; | |
144 | fLoMultBound = -1; | |
145 | fHiMultBound = 10000; | |
146 | ||
147 | //Default: no cut in Armenteros | |
148 | fSpecialArmenterosCutK0s = kFALSE; | |
149 | ||
57a6d361 | 150 | //Pt Bins: undefined |
151 | fptbinnumb = -1; | |
152 | } | |
153 | ||
154 | /*********************************************** | |
155 | --- Setters For Configuration --- | |
156 | ***********************************************/ | |
157 | ||
158 | // Filename Setters | |
159 | void AliV0Module::SetRealDataFile ( TString RealDataFilename ){ | |
160 | //Set root file containing real data candidates. | |
161 | fRealDataFile = RealDataFilename; } | |
162 | void AliV0Module::SetMCDataFile ( TString MCDataFilename ){ | |
163 | //Set root file containing Monte Carlo data (for efficiency computation). | |
164 | fMCDataFile = MCDataFilename; } | |
165 | void AliV0Module::SetFeedDownDataFile( TString FeedDownDataFilename ){ | |
166 | //Set root file containing Monte Carlo data (for feeddown computation). | |
167 | fFeedDownDataFile = FeedDownDataFilename; } | |
168 | void AliV0Module::SetOutputFile ( TString OutputFilename ){ | |
169 | //Set root filename for the analysis output. | |
170 | fOutputDataFile = OutputFilename; | |
171 | cout<<"[AliV0Module] Set output file to \'"<<OutputFilename<<"\'."<<endl; | |
172 | } | |
173 | ||
174 | // Bin Limit Setter | |
175 | void AliV0Module::SetPtBinLimits(Long_t got_fptbinnumb,const Double_t *got_fptbinlimits){ | |
176 | //Function to set pt binning. First argument is the number of pt bins, second is | |
177 | //an array with bin limits. | |
178 | fptbinnumb = got_fptbinnumb; | |
179 | for(int ix = 0;ix<fptbinnumb+1;ix++){ | |
180 | fptbinlimits[ix] = got_fptbinlimits[ix]; | |
181 | } | |
182 | for(int ix = 0;ix<fptbinnumb;ix++){ | |
183 | fptX[ix] = (fptbinlimits[ix+1] + fptbinlimits[ix])/2.; | |
184 | } | |
185 | cout<<"[AliV0Module] Received "<<fptbinnumb<<" pt bins, set accordingly."<<endl; | |
186 | } | |
187 | ||
188 | // Rapidity Window Setter | |
189 | void AliV0Module::SetRapidityWindow(Double_t got_fRapidityBoundary){ | |
190 | //Set Rapidity Boundary used in analysis. | |
191 | //Value provided will be used as upper limit for the modulus of y (|y|< given value) | |
192 | fRapidityBoundary = got_fRapidityBoundary; | |
193 | cout<<"[AliV0Module] Received "<<got_fRapidityBoundary<<" as rapidity limits, set accordingly."<<endl; | |
194 | } | |
195 | ||
196 | // CINT1B/INEL Setter for normalization to yields | |
197 | void AliV0Module::SetCINT1BoverINEL(Double_t got_fCINT1BoverINEL){ | |
198 | //Set CINT1B/INEL ratio (determined by collaboration). | |
199 | fCINT1BoverINELratio = got_fCINT1BoverINEL; | |
200 | cout<<"[AliV0Module] Received CINT1B/INEL = "<<got_fCINT1BoverINEL<<" to normalize with, set accordingly."<<endl; | |
201 | } | |
202 | ||
203 | // Topological Selection Setters | |
204 | void AliV0Module::SetCutV0Radius(Double_t cut){ | |
205 | //Set minimum decay radius for the V0 in centimeters. | |
206 | //Note: this is (R_{2D}, cylindrical, centered around ALICE detector) | |
207 | fCutV0Radius = cut; | |
208 | cout<<"[AliV0Module] Received V0 Radius (min value) = "<<cut<<endl; | |
209 | } | |
210 | void AliV0Module::SetCutDCANegToPV(Double_t cut){ | |
211 | //Set minimum distance of closest approach between V0 negative daughter | |
212 | //track and primary vertex (in centimeters). | |
213 | fCutDCANegToPV = cut; | |
214 | cout<<"[AliV0Module] Received DCA Negative track to PV (min value) = "<<cut<<endl; | |
215 | } | |
216 | void AliV0Module::SetCutDCAPosToPV(Double_t cut){ | |
217 | //Set minimum distance of closest approach between V0 positive daughter | |
218 | //track and primary vertex (in centimeters). | |
219 | fCutDCAPosToPV = cut; | |
220 | cout<<"[AliV0Module] Received DCA Positive track to PV (min value) = "<<cut<<endl; | |
221 | } | |
222 | void AliV0Module::SetCutDCAV0Daughters(Double_t cut){ | |
223 | //Set minimum distance of closest approach between V0 daughter | |
224 | //tracks, in sigmas. This is not in centimeters because if the | |
225 | //tracks have been determined with ITS refit the resolution is | |
226 | //greatly improved; thus, the cut may be tighter in that case. | |
227 | //Using sigmas will therefore take the tracking method in | |
228 | //consideration. | |
229 | fCutDCAV0Daughters = cut; | |
230 | cout<<"[AliV0Module] Received DCA V0 Daughters (max value) = "<<cut<<endl; | |
231 | } | |
232 | void AliV0Module::SetCutV0CosPA(Double_t cut){ | |
233 | //Set minimum value for the cosine of pointing angle of the V0. | |
234 | fCutV0CosPA = cut; | |
235 | cout<<"[AliV0Module] Received V0 Cosine of Pointing Angle (min value) = "<<cut<<endl; | |
236 | } | |
237 | ||
238 | // Other Selection Setters | |
239 | void AliV0Module::SetCutProperLifetime(Double_t cut){ | |
240 | //Set maximum value for m*L/p variable for the V0. | |
241 | //This is the "proper lifetime" selection and is usually called a | |
242 | //"c*tau cut". Should be set to a value larger than the c*tau for | |
243 | //the V0 considered. | |
244 | fCutProperLifetime = cut; | |
245 | cout<<"[AliV0Module] Received proper lifetime cut (max value) = "<<cut<<endl; | |
246 | } | |
247 | void AliV0Module::SetCutTPCPIDNSigmas(Double_t cut){ | |
248 | //Set maximum deviation from the expected energy loss in | |
249 | //the TPC, in multiples of sigmas as computed from the AliPIDResponse | |
250 | //object. Selection is only used in real data and should thus | |
251 | //be very loose to ensure negligible signal loss. | |
252 | fCutTPCPIDNSigmas = cut; | |
253 | cout<<"[AliV0Module] Received TPC N-sigmas selection (max dist from BB curve) = "<<cut<<endl; | |
254 | } | |
255 | void AliV0Module::SetCutSigmaForSignalExtraction(Double_t cut){ | |
256 | //Set number of sigmas for the signal extraction method. The value | |
257 | //provided is the number of sigmas from the peak position used for | |
258 | //the peak sampling, i.e. the peak will be from [<m>-cut*sigma,<m>+cut*sigma] | |
259 | //while the background sampling regions will be from | |
260 | //[<m>-2*cut*sigma,<m>+2*cut*sigma. | |
261 | fCutNSigmasForSignalExtraction = cut; | |
262 | cout<<"[AliV0Module] Received N-sigmas for sig. ext.: peak is (-"<<cut<<",+"<<cut<<") in sigmas"<<endl; | |
263 | } | |
264 | void AliV0Module::SetCutLeastNumberOfCrossedRows(Double_t cut){ | |
265 | //Set smallest allowed number of TPC clusters for the V0 daughter tracks. | |
266 | fCutLeastNumberOfCrossedRows = cut; | |
267 | cout<<"[AliV0Module] Received Least Nbr of crossed rows (min value) = "<<cut<<endl; | |
268 | } | |
269 | void AliV0Module::SetCutLeastNumberOfCrossedRowsOverFindable(Double_t cut){ | |
270 | //Set smallest allowed number of TPC clusters for the V0 daughter tracks. | |
271 | fCutLeastNumberOfCrossedRowsOverFindable = cut; | |
272 | cout<<"[AliV0Module] Received Least Nbr of crossed rows over findable clusters (min value) = "<<cut<<endl; | |
273 | } | |
274 | ||
275 | void AliV0Module::SetCutDaughterEta(Double_t cut){ | |
276 | //Set maximum eta value allowed for the V0 daughter tracks (|eta|<cut). | |
277 | fCutDaughterEta = cut; | |
278 | cout<<"[AliV0Module] Received Daughter |eta| cut (max value) = "<<cut<<endl; | |
279 | } | |
280 | ||
281 | void AliV0Module::SetCutCompetingV0Rejection(Double_t cut){ | |
282 | //Set rejection window around invariant mass of competing V0 species. | |
283 | //If negative, no rejection will occur. If analysis is for Lambdas, K0s | |
284 | //will be rejected and vice-versa. Typical values revolve around | |
285 | //0.003 - 0.010 GeV/c^2. | |
286 | fCutCompetingV0Rejection = cut; | |
287 | cout<<"[AliV0Module] Received Competing V0 Rejection window of +/- (in GeV/c^2) = "<<cut<<endl; | |
288 | } | |
289 | ||
290 | void AliV0Module::SetFeeddownTreatment ( TString fFDMethod ){ | |
291 | //Set method used to compute charged and neutral Xi baryon feeddown to Lambda. | |
292 | //Methods allowed are: | |
293 | //"NoFD" - No Feeddown subtraction performed. | |
294 | //"DoubleChargedXi" - Feeddown is computed for the charged Xi, and is then | |
295 | //subtracted twice from real data. Assumes charged and neutral Xi enter | |
296 | //the analysis in the same way and are produced at a ratio 1:1. | |
297 | //"UseMCRatio" - The Feeddown matrix F_{ij} is filled not only with Lambdas | |
298 | //coming from charged Xis but also from neutral ones. The scaling is performed | |
299 | //according to the measured charged Xis, and since the matrix is more populated | |
300 | //a larger subtraction will occurr. This method assumes that MC correctly | |
301 | //reproduces the ratio between charged and neutral Xi baryons. Warning: this | |
302 | //should not be used with charged Xi triggered Monte Carlo data. | |
303 | fFDSwitch = fFDMethod; | |
304 | cout<<"[AliV0Module] Received Feeddown treatment method: "<<fFDMethod<<endl; | |
305 | if( fFDMethod == "NoFD" ) | |
306 | cout<<"[AliV0Module] ---> No Feeddown correction will be performed."<<endl; | |
307 | if( fFDMethod == "DoubleChargedXi" ) | |
308 | cout<<"[AliV0Module] ---> Feeddown performed by doubling charged Xi correction."<<endl; | |
309 | if( fFDMethod == "UseMCRatio" ) | |
310 | cout<<"[AliV0Module] ---> Feeddown performed by using MC neutral/charged Xi."<<endl; | |
311 | } | |
312 | ||
313 | void AliV0Module::SetFitBackground ( Bool_t fitBgSwitch ){ | |
314 | //Turns on background fitting for signal extraction instead of pure | |
315 | //bin counting. Useful, among other things, for systematics. | |
316 | fFitBackgroundSwitch = fitBgSwitch; | |
317 | } | |
318 | ||
d1d5d7e6 | 319 | void AliV0Module::SetPerformMultiplicityStudy ( Bool_t lPerformMultStudy ){ |
320 | //Turns on selection according to multiplicity of the event. | |
321 | //Boundaries are set in charged track multiplicity (pp) or in | |
322 | //centrality percentiles (PbPb). This is a requirement for studying | |
323 | //PbPb. | |
324 | fPerformMultiplicityStudy = lPerformMultStudy; | |
325 | } | |
326 | ||
327 | void AliV0Module::SetLowMultValue ( Int_t lLoMultBound ){ | |
328 | //Lower boundary (inclusive) in integer number for mult selection. | |
329 | //Note: If in PbPb and you want, say, 10-20% centrality, set this | |
330 | //to 10. | |
331 | fLoMultBound = lLoMultBound; | |
332 | } | |
333 | ||
334 | void AliV0Module::SetHighMultValue ( Int_t lHiMultBound ){ | |
335 | //Lower boundary (inclusive) in integer number for mult selection. | |
336 | //Note: If in PbPb and you want, say, 10-20% centrality, set this | |
337 | //to 10. | |
338 | fHiMultBound = lHiMultBound; | |
339 | } | |
340 | ||
341 | void AliV0Module::SetSpecialArmenterosCutK0s ( Bool_t lSpecialArmenterosCutK0s ){ | |
342 | //Special armenteros cut: |alpha|<5*pt_{arm} | |
343 | fSpecialArmenterosCutK0s = lSpecialArmenterosCutK0s; | |
344 | } | |
345 | ||
57a6d361 | 346 | void AliV0Module::SetDefaultCuts(){ |
347 | //Sets Default cuts for analysis. (adjusted for adequate pp analysis) | |
348 | cout<<"[AliV0Module] Setting default cuts for particle species: "<<fWhichParticle<<endl; | |
349 | //Set Cuts - topological | |
350 | SetCutV0Radius (0.500); | |
351 | SetCutDCANegToPV (0.060); | |
352 | SetCutDCAPosToPV (0.060); | |
353 | SetCutDCAV0Daughters (1.000); | |
354 | if ( fWhichParticle == "K0Short") | |
355 | SetCutV0CosPA (0.970); | |
356 | if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ) | |
357 | SetCutV0CosPA (0.995); | |
358 | ||
359 | ||
360 | //Set Cuts - other | |
d1d5d7e6 | 361 | if ( fWhichParticle == "K0Short") |
362 | SetCutProperLifetime ( 20); | |
363 | if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ) | |
364 | SetCutProperLifetime ( 30); | |
365 | ||
57a6d361 | 366 | SetCutTPCPIDNSigmas ( 5); |
d1d5d7e6 | 367 | SetCutSigmaForSignalExtraction ( 6); |
57a6d361 | 368 | SetCutLeastNumberOfCrossedRows ( 70); |
369 | SetCutLeastNumberOfCrossedRowsOverFindable (0.8); | |
370 | SetCutDaughterEta (0.8); | |
371 | if ( fWhichParticle == "K0Short") | |
372 | SetCutCompetingV0Rejection (0.005); | |
373 | if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ) | |
374 | SetCutCompetingV0Rejection (0.010); | |
375 | } | |
376 | ||
377 | TString AliV0Module::IntToString(int input){ | |
378 | //Integer to TString Converter | |
379 | char dummyChar[50]; | |
380 | sprintf(dummyChar, "%d", input); | |
381 | TString outputstring = dummyChar; | |
382 | return outputstring; | |
383 | } | |
384 | ||
385 | TString AliV0Module::DoubleToString(double input){ | |
386 | //Double to TString Converter | |
387 | char dummyChar[50]; | |
388 | sprintf(dummyChar, "%.3f", input); | |
389 | TString outputstring = dummyChar; | |
390 | return outputstring; | |
391 | } | |
392 | ||
393 | Double_t AliV0Module::ErrorInRatio ( Double_t A, Double_t Aerr, Double_t B, Double_t Berr ){ | |
394 | //Error in a Ratio | |
395 | if(B!=0){ | |
396 | Double_t errorfromtop = Aerr*Aerr / (B*B) ; | |
397 | Double_t errorfrombottom = ((A*A)/(B*B*B*B)) * Berr * Berr; | |
398 | return TMath::Sqrt( errorfromtop + errorfrombottom ); | |
399 | } | |
400 | return 1; | |
401 | } | |
402 | ||
403 | Double_t AliV0Module::MyGeant3FlukaCorrectionForProtons(const Double_t *x, const Double_t *par){ | |
404 | //Parametrization Used for Geant3/Fluka Correction for protons | |
405 | //Credit: Antonin Maire | |
406 | return 1 - par[0]*TMath::Exp(par[1]*x[0]) + par[2]; | |
407 | } | |
408 | ||
409 | ||
410 | Double_t AliV0Module::MyGeant3FlukaCorrectionForAntiProtons(const Double_t *x, const Double_t *par){ | |
411 | // Parametrization Used for Geant3/Fluka Correction for antiprotons | |
412 | // Credit: Antonin Maire | |
413 | // | |
414 | // La fonction A*TMath::Erf( B.x ) ne marche pas à bas pt. | |
415 | // Différentes fonctions ont été testées. | |
416 | // La fonction suivante semble donner de bons résultats | |
417 | // On peut jouer sur la puissance n du terme "1/x^n" pour repousser le pt auquel la fonction prend la valeur 1.0 | |
418 | // Ici, pour n = 0.2, la fonction prend la valeur 0.9990 en pt = 10 GeV/c | |
419 | ||
420 | return 1 - par[0]*TMath::Exp(par[1]*x[0]) + par[2] + par[3]*1/TMath::Power(x[0], 0.2)*TMath::Log(x[0]); | |
421 | } | |
422 | ||
14d51191 | 423 | Double_t AliV0Module::MyLevyPtXi(const Double_t *pt, const Double_t *par) |
57a6d361 | 424 | { |
425 | //Levy Fit Function | |
426 | Double_t lMass = 1.32171; //Xi Mass | |
427 | Double_t ldNdy = par[0]; | |
428 | Double_t lTemp = par[1]; | |
429 | Double_t lPower = par[2]; | |
430 | ||
431 | Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (lPower*lTemp*(lPower*lTemp+lMass*(lPower-2))); | |
432 | Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp); | |
433 | ||
434 | return ldNdy * pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower); | |
435 | //return ldNdy * lBigCoef * TMath::Power(lInPower,(-1)*lPower); | |
436 | } | |
437 | ||
438 | Double_t AliV0Module::MyBgPol1(const Double_t *x, const Double_t *par) | |
439 | { | |
440 | //Function for background fitting, rejects peak region | |
441 | if ( x[0] > par[2] && x[0] < par[3]) { | |
442 | TF1::RejectPoint(); | |
443 | return 0; | |
444 | } | |
445 | return par[0] + par[1]*x[0]; | |
446 | } | |
447 | ||
448 | Double_t AliV0Module::MyBgPolToEval1(const Double_t *x, const Double_t *par) | |
449 | { | |
450 | //Just a plain linear function. | |
451 | return par[0] + par[1]*x[0]; | |
452 | } | |
453 | ||
454 | Double_t AliV0Module::RoundToThousandth( const Double_t lToRound ){ | |
455 | //Round any number to a hundredth... | |
456 | //Well, within machine precision... | |
457 | return TMath::Nint( 1000 * lToRound ) * 0.001; | |
458 | } | |
459 | ||
460 | void AliV0Module::DoAnalysis(){ | |
461 | //---------------------- | |
462 | // V0 Analysis Function | |
463 | //---------------------- | |
464 | // | |
465 | //Consists of the following steps: | |
466 | // | |
467 | // (1) Loop over Real data candidates, acquire peak position and widths. | |
468 | // (2) Loop over Real data candidates, extract signal with variable extraction | |
469 | // areas. Two loops: totally avoids binning data, allowing for sub-MeV/c^2 | |
470 | // granularity in signal extraction areas. | |
471 | // (3) Loop over MC data reconstructed candidates, find associated-to-MC primary | |
472 | // candidates for efficiency numerator. | |
473 | // (4) (if Lambda or AntiLambda) Perform Geant3/Fluka correction. | |
474 | // (5) Get generated primary V0 histograms from MC file for efficiency denominator. | |
475 | // (6) (if Lambda or AntiLambda + FD correction enabled) Open MC file for feeddown | |
476 | // subtraction. Loop over candidates, find Lambda associated with primary Xi | |
477 | // to fill feeddown matrix. Scale the feeddown contribution according to real-life | |
478 | // measured Xi- production under specified feeddown subtraction scheme (see | |
479 | // function SetFeeddownTreatment). Perform subtraction of raw Lambda or AntiLambda | |
480 | // count estimated to be coming from charged or neutral Xi. | |
481 | // (7) Perform detection efficiency correction and compute final corrected spectra. | |
482 | // | |
483 | // | |
484 | // Normalization: | |
485 | // --- Number of Inelastic Events estimated to be: | |
486 | // <Number of CINT1B triggers> / <CINT1B/INEL ratio> | |
487 | // | |
488 | // --- Efficiency denominator: filled at pre-physics selection stage. All signal | |
489 | // losses are therefore estimated from Monte Carlo. | |
490 | // | |
491 | // --- Pileup: not included in MC. Fraction of events removed by IsPileupFromSPD | |
492 | // is artificially removed from Number of Inelastic Events as well. | |
493 | // | |
494 | // Output: written to specified file (with SetOutputFile(...)). | |
495 | // --- includes spectra and efficiencies in base directory. | |
496 | // --- includes a number of different subdirectories storing plots such as | |
497 | // such as invariant mass histos, signal extraction, pt resolution, G3/F | |
498 | // related distributions, feeddown related distributions, and so on. | |
499 | ||
500 | //Set Batch Mode: Ignore all canvases in this section | |
501 | gROOT->SetBatch (kTRUE); | |
502 | ||
503 | cout<<"======================================"<<endl; | |
504 | cout<<" -------- Spectra Extraction -------- "<<endl; | |
505 | cout<<"======================================"<<endl; | |
506 | cout<<endl; | |
507 | if(fptbinnumb == -1){ | |
508 | cout<<"[AliV0Module] It's unclear what kind of pt binning you want."<<endl; | |
509 | cout<<"[AliV0Module] Most likely you forgot to set it using SetPtBinLimits..."<<endl; | |
510 | cout<<"[AliV0Module] Analysis will NOT be done. Returning."<<endl; | |
511 | return; | |
512 | } | |
513 | cout<<"--------------- Configuration --------------------------"<<endl; | |
514 | cout<<" Analysed Particle.............: "<<fWhichParticle<<endl; | |
515 | cout<<" Rapidity Window...............: "<<fRapidityBoundary<<endl; | |
516 | cout<<" CINT1B/INEL Ratio used........: "<<fCINT1BoverINELratio<<endl; | |
517 | cout<<" V0 Decay Radius...............: "<<fCutV0Radius<<endl; | |
518 | cout<<" DCA Negative track to PV......: "<<fCutDCANegToPV<<endl; | |
519 | cout<<" DCA Positive track to PV......: "<<fCutDCAPosToPV<<endl; | |
520 | cout<<" DCA V0 Daughters..............: "<<fCutDCAV0Daughters<<endl; | |
521 | cout<<" Cosine of Pointing Angle V0...: "<<fCutV0CosPA<<endl; | |
522 | cout<<" Proper Lifetime cut (cm)......: "<<fCutProperLifetime<<endl; | |
523 | cout<<" TPC dE/dx sigmas cut (Real)...: "<<fCutTPCPIDNSigmas<<endl; | |
524 | cout<<" CutNSigmasForSignalExtraction.: "<<fCutNSigmasForSignalExtraction<<endl; | |
525 | cout<<" Least # of Crossed Rows.......: "<<fCutLeastNumberOfCrossedRows<<endl; | |
526 | cout<<" Least # of C. R. over findable: "<<fCutLeastNumberOfCrossedRowsOverFindable<<endl; | |
527 | cout<<" Daughter Track |eta| < .......: "<<fCutDaughterEta<<endl; | |
528 | cout<<" Competing V0 Reject. (GeV/c^2): "<<fCutCompetingV0Rejection<<endl; | |
529 | cout<<"--------------- File Names -----------------------------"<<endl; | |
530 | cout<<" Real Data File................: "<<fRealDataFile<<endl; | |
531 | cout<<" MC File.......................: "<<fMCDataFile<<endl; | |
532 | cout<<" MC File for feeddown..........: "<<fFeedDownDataFile<<endl; | |
533 | cout<<" Analysis output file..........: "<<fOutputDataFile<<endl; | |
534 | cout<<"--------------------------------------------------------"<<endl; | |
535 | cout<<endl; | |
536 | ||
537 | //=== Real data loop 1: Acquire peak positions, widths === | |
538 | //--- Preparing... --- | |
539 | //defining helping histogram - only used to FindBin Index!======== | |
540 | TH1F* fHistPt = new TH1F("fHistPt","Dummy;p_{T} (GeV/c);Counts",fptbinnumb,fptbinlimits); | |
541 | //Peak Position Histograms | |
542 | TH1F* fHistPeakPosition = new TH1F("fHistPeakPosition","Peak Position (Real);p_{T} (GeV/c);Peak Position",fptbinnumb,fptbinlimits); | |
543 | TH1F* fHistPeakPositionMC = new TH1F("fHistPeakPositionMC","Peak Position (MC);p_{T} (GeV/c);Peak Position",fptbinnumb,fptbinlimits); | |
544 | //Signal to Noise Histograms | |
545 | TH1F* fHistSigToNoise = new TH1F("fHistSigToNoise","Signal to Noise Ratio (real);p_{T} (GeV/c);Sig / Bg",fptbinnumb,fptbinlimits); | |
546 | TH1F* fHistSigToNoiseMC = new TH1F("fHistSigToNoiseMC","Signal to Noise Ratio (MC);p_{T} (GeV/c);Sig / Bg",fptbinnumb,fptbinlimits); | |
547 | ||
548 | //Signal Extraction Range Histogram | |
549 | TH1F* fHistSignalExtractionRange = new TH1F("fHistSignalExtractionRange","Sig. Ext. Range;p_{T} (GeV/c);Range",fptbinnumb,fptbinlimits); | |
550 | ||
551 | //Resolution Histogram (filled with MC) | |
552 | TH2F* f2dHistPtResolution = new TH2F("f2dHistPtResolution","p_{t} Resolution;p_{t} (reco);p_{t} (mc)",fptbinnumb,fptbinlimits, fptbinnumb, fptbinlimits); | |
553 | ||
d1d5d7e6 | 554 | TH2F* f2dHistPtBlur = new TH2F("f2dHistPtBlur","p_{t} blurring matrix;p_{t} (reco);p_{t} (mc)",fptbinnumb,fptbinlimits, fptbinnumb, fptbinlimits); |
555 | TH2F* f2dHistPtSharpen = new TH2F("f2dHistPtSharpen","p_{t} sharpening matrix;p_{t} (reco);p_{t} (mc)",fptbinnumb,fptbinlimits, fptbinnumb, fptbinlimits); | |
556 | ||
57a6d361 | 557 | /******************************************************** |
558 | ||
559 | ---> Let's Remember the limits of the data we're analyzing! | |
560 | ---> Important so that we don't try signal extraction outside | |
561 | ---> the bundaries of available data in the Tree object. | |
562 | ||
563 | From AliAnalysisTaskExtractV0.cxx | |
564 | ||
565 | //Second Selection: rough 20-sigma band, parametric. | |
566 | //K0Short: Enough to parametrize peak broadening with linear function. | |
567 | Double_t UpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*tree_lPt; | |
568 | Double_t LowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*tree_lPt; | |
569 | ||
570 | //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening) | |
571 | //[0]+[1]*x+[2]*TMath::Exp(-[3]*x) | |
572 | Double_t UpperLimitLambda = (1.13688e+00) + (5.27838e-03)*tree_lPt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*tree_lPt); | |
573 | Double_t LowerLimitLambda = (1.09501e+00) - (5.23272e-03)*tree_lPt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*tree_lPt); | |
574 | ||
575 | ********************************************************/ | |
576 | ||
577 | TF1 *fKDataUpper = new TF1("fKDataUpper","[0]+[1]*x",0,20); | |
578 | TF1 *fKDataLower = new TF1("fKDataLower","[0]+[1]*x",0,20); | |
579 | TF1 *fLDataUpper = new TF1("fLDataUpper","[0]+[1]*x+[2]*TMath::Exp([3]*x)",0,20); | |
580 | TF1 *fLDataLower = new TF1("fLDataLower","[0]+[1]*x+[2]*TMath::Exp([3]*x)",0,20); | |
581 | ||
582 | fKDataUpper->SetParameter(0, 5.63707e-01); | |
583 | fKDataUpper->SetParameter(1, 1.14979e-02); | |
584 | fKDataLower->SetParameter(0, 4.30006e-01); | |
585 | fKDataLower->SetParameter(1,-1.10029e-02); | |
586 | ||
587 | fLDataUpper->SetParameter(0, 1.13688e+00); | |
588 | fLDataUpper->SetParameter(1, 5.27838e-03); | |
589 | fLDataUpper->SetParameter(2, 8.42220e-02); | |
590 | fLDataUpper->SetParameter(3,-3.80595e+00); | |
591 | ||
592 | fLDataLower->SetParameter(0, 1.09501e+00); | |
593 | fLDataLower->SetParameter(1,-5.23272e-03); | |
594 | fLDataLower->SetParameter(2,-7.52690e-02); | |
595 | fLDataLower->SetParameter(3,-3.46339e+00); | |
596 | ||
597 | Int_t lWeAreAtBin = 0; | |
598 | Double_t lParticleMass = -1; //needed for proper lifetime selection | |
599 | Double_t lCompetingParticleMass = -1; //needed for Competing V0 Rejection | |
600 | Double_t lHistoLowerBoundary = -1; | |
601 | Double_t lHistoUpperBoundary = -1; | |
602 | Long_t lHistoNBins = -1; | |
603 | if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
604 | lParticleMass = 1.115683; | |
605 | lCompetingParticleMass = 0.4976; | |
606 | lHistoLowerBoundary = 1.116 - 0.1; | |
607 | lHistoUpperBoundary = 1.116 + 0.1; | |
608 | lHistoNBins = 200; //1MeV/c^2 binning | |
609 | } | |
610 | if ( fWhichParticle == "K0Short" ){ | |
611 | lParticleMass = 0.4976; | |
612 | lCompetingParticleMass = 1.115683; | |
613 | lHistoLowerBoundary = 0.498 - 0.15; | |
614 | lHistoUpperBoundary = 0.498 + 0.15; | |
615 | lHistoNBins = 300; //1MeV/c^2 binning | |
616 | } | |
617 | //Setting Up Base Histograms to use=============================== | |
618 | TH1F* lHistoFullV0[100]; | |
619 | TH1F* lHistoFullV0MC[100]; | |
620 | TCanvas* lCanvasHistoFullV0[100]; | |
621 | TCanvas* lCanvasHistoFullV0MC[100]; | |
d1d5d7e6 | 622 | TH1F *lHistResolution[100]; |
623 | TF1 *lHistResolutionGaussian[100]; | |
624 | char lHistResolutionName[100]; | |
625 | TH1F* fHistResolutionVsPt = new TH1F("fHistResolutionVsPt","Resolution vs p_{t};p_{t} (GeV/c);#Delta(p_{t}) = p_{t}-p_{t}^{true} (GeV/c)",fptbinnumb,fptbinlimits); | |
626 | TH1F* fHistResolutionVsPtDivByBinWidth = new TH1F("fHistResolutionVsPtDivByBinWidth","Resolution vs p_{t};p_{t} (GeV/c);#Delta(p_{t})/#Delta^{bin}(p_{t}) = (p_{t}-p_{t}^{true})/#Delta^{bin}(p_{t}) (GeV/c)",fptbinnumb,fptbinlimits); | |
627 | TH1F* fHistResolutionVsPtWithGaussians = new TH1F("fHistResolutionVsPtWithGaussians","Resolution vs p_{t};p_{t} (GeV/c);#Delta(p_{t}) = p_{t}-p_{t}^{true} (GeV/c)",fptbinnumb,fptbinlimits); | |
628 | TH1F* fHistResolutionVsPtDivByBinWidthWithGaussians = new TH1F("fHistResolutionVsPtDivByBinWidthWithGaussians","Resolution vs p_{t};p_{t} (GeV/c);#Delta(p_{t})/#Delta^{bin}(p_{t}) = (p_{t}-p_{t}^{true})/#Delta^{bin}(p_{t}) (GeV/c)",fptbinnumb,fptbinlimits); | |
629 | for(Int_t ibin=0;ibin<fptbinnumb;ibin++){ | |
630 | } | |
57a6d361 | 631 | char histname[80]; TString bindescription = ""; |
632 | for(Int_t ihist=0;ihist<100;ihist++){ | |
633 | //Histo For Real Data | |
634 | sprintf(histname,"lHistoFullV0%i",ihist); | |
635 | if(fWhichParticle == "Lambda") bindescription="#Lambda, bin #"; | |
636 | if(fWhichParticle == "AntiLambda") bindescription="#bar{#Lambda}, bin #"; | |
637 | if(fWhichParticle == "K0Short") bindescription="K^{0}_{S}, bin #"; | |
638 | bindescription.Append(IntToString( ihist )); | |
639 | if ( ihist < fptbinnumb ){ | |
640 | bindescription.Append(", "); | |
641 | bindescription.Append(DoubleToString(fptbinlimits[ ihist ])); | |
642 | bindescription.Append("-"); | |
643 | bindescription.Append(DoubleToString(fptbinlimits[ ihist+1 ])); | |
644 | bindescription.Append("GeV/c"); | |
645 | } | |
646 | lHistoFullV0[ihist] = new TH1F(histname,"Candidates;M(appropriate) (GeV/c^{2});Counts", lHistoNBins,lHistoLowerBoundary,lHistoUpperBoundary); | |
647 | lHistoFullV0[ihist]->SetTitle(bindescription); | |
648 | sprintf(histname,"lCanvasHistoFullV0%i",ihist); | |
649 | lCanvasHistoFullV0[ihist] = new TCanvas(histname, bindescription, 800, 600); | |
650 | lCanvasHistoFullV0[ihist] -> SetFillColor(kWhite); | |
651 | lCanvasHistoFullV0[ihist] -> SetLeftMargin( 0.175 ); | |
652 | lCanvasHistoFullV0[ihist] -> SetBottomMargin( 0.175 ); | |
57a6d361 | 653 | |
654 | //Histo for MC | |
655 | sprintf(histname,"lHistoFullV0MC%i",ihist); | |
656 | if(fWhichParticle == "Lambda") bindescription="#Lambda, MC, bin #"; | |
657 | if(fWhichParticle == "AntiLambda") bindescription="#bar{#Lambda}, MC, bin #"; | |
658 | if(fWhichParticle == "K0Short") bindescription="K^{0}_{S}, MC, bin #"; | |
659 | bindescription.Append(IntToString( ihist )); | |
660 | if ( ihist < fptbinnumb ){ | |
661 | bindescription.Append(", "); | |
662 | bindescription.Append(DoubleToString(fptbinlimits[ ihist ])); | |
663 | bindescription.Append("-"); | |
664 | bindescription.Append(DoubleToString(fptbinlimits[ ihist+1 ])); | |
665 | bindescription.Append("GeV/c"); | |
666 | } | |
667 | lHistoFullV0MC[ihist] = new TH1F(histname,"Candidates;M(appropriate) (GeV/c^{2});Counts", lHistoNBins,lHistoLowerBoundary,lHistoUpperBoundary); | |
668 | lHistoFullV0MC[ihist]->SetTitle(bindescription); | |
669 | sprintf(histname,"lCanvasHistoFullV0MC%i",ihist); | |
670 | lCanvasHistoFullV0MC[ihist] = new TCanvas(histname, bindescription, 800, 600); | |
671 | lCanvasHistoFullV0MC[ihist] -> SetFillColor(kWhite); | |
672 | lCanvasHistoFullV0MC[ihist] -> SetLeftMargin( 0.175 ); | |
673 | lCanvasHistoFullV0MC[ihist] -> SetBottomMargin( 0.175 ); | |
d1d5d7e6 | 674 | |
675 | //Histo for resolution tests | |
676 | sprintf(lHistResolutionName,"lHistResolution%i",ihist); | |
677 | Long_t lNumberOfBinsResolution = 5000; | |
678 | if ( (fptbinlimits[ihist+1]+fptbinlimits[ihist])*0.5 > 5 ) | |
679 | lNumberOfBinsResolution = 500; | |
680 | if ( (fptbinlimits[ihist+1]+fptbinlimits[ihist])*0.5 > 10 ) | |
681 | lNumberOfBinsResolution = 50; | |
682 | ||
683 | lHistResolution[ihist] = new TH1F ( lHistResolutionName,bindescription,lNumberOfBinsResolution, -5, 5); //histo with 5MeV/c precision! | |
684 | sprintf(lHistResolutionName,"lHistResolutionGaussian%i",ihist); | |
685 | lHistResolutionGaussian[ihist] = new TF1(lHistResolutionName, "[0]*TMath::Gaus(x,[1],[2])",-5,5); | |
686 | ||
57a6d361 | 687 | } |
688 | //================================================================ | |
689 | cout<<endl; | |
690 | ||
691 | ||
692 | cout<<"--------------- Open Real Data File --------------------"<<endl; | |
693 | TFile* file = TFile::Open(fRealDataFile, "READ"); | |
694 | file->cd("PWG2CheckLambda_PP"); | |
695 | TList* v0list = (TList*)file->FindObjectAny("clistV0"); | |
696 | TTree* lTree; | |
697 | TH1F* fHistV0MultiplicityForTrigEvt; | |
698 | TH1F* fHistV0MultiplicityForSelEvtNoTPCOnly; | |
699 | TH1F* fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup; | |
700 | lTree = (TTree*)file->FindObjectAny("fTree"); | |
701 | fHistV0MultiplicityForTrigEvt = (TH1F*)v0list->FindObject("fHistV0MultiplicityForTrigEvt"); | |
702 | fHistV0MultiplicityForSelEvtNoTPCOnly = (TH1F*)v0list->FindObject("fHistV0MultiplicityForSelEvtNoTPCOnly"); | |
703 | fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = (TH1F*)v0list->FindObject("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup"); | |
704 | ||
705 | Long_t lNCandidates = lTree->GetEntries(); | |
706 | Long_t lNEvents = fHistV0MultiplicityForTrigEvt->GetEntries(); | |
707 | Double_t lEvtFracAfterPileup = | |
708 | ((double)(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup->GetEntries()))/ | |
709 | ((double)(fHistV0MultiplicityForSelEvtNoTPCOnly->GetEntries())); | |
710 | ||
711 | Double_t lNInelasticEvents = lEvtFracAfterPileup * (lNEvents / fCINT1BoverINELratio); | |
712 | ||
713 | cout<<" Number of CINT1B triggers.....: "<<lNEvents<<endl; | |
714 | cout<<" Pileup Fraction...............: "<<1 - lEvtFracAfterPileup<<endl; | |
715 | cout<<" Estimated N_{INEL}............: "<<lNInelasticEvents<<endl; | |
716 | cout<<" V0 Candidates.................: "<<lNCandidates<<endl; | |
717 | cout<<"--------------------------------------------------------"<<endl; | |
718 | ||
d1d5d7e6 | 719 | Double_t lNEventsBeforePileup = 0; |
720 | Double_t lNEventsAfterPileup = 0; | |
721 | if( fPerformMultiplicityStudy == kTRUE ){ | |
722 | lNInelasticEvents = 0; | |
723 | ||
724 | TH1F* fHistMultiplicity = (TH1F*)v0list->FindObject("fHistMultiplicityForTrigEvt"); | |
725 | TH1F* fHistMultiplicityBeforePileup = (TH1F*)v0list->FindObject("fHistMultiplicityNoTPCOnly"); | |
726 | TH1F* fHistMultiplicityAfterPileup = (TH1F*)v0list->FindObject("fHistMultiplicityNoTPCOnlyNoPileup"); | |
727 | ||
728 | ||
729 | ||
730 | for( Int_t ibin = fHistMultiplicity->FindBin(fLoMultBound); ibin<fHistMultiplicity->FindBin(fHiMultBound)+1; ibin++){ | |
731 | lNInelasticEvents += fHistMultiplicity->GetBinContent(ibin); | |
732 | lNEventsBeforePileup += fHistMultiplicityBeforePileup->GetBinContent(ibin); | |
733 | lNEventsAfterPileup += fHistMultiplicityAfterPileup ->GetBinContent(ibin); | |
734 | } | |
735 | cout<<" Number of events, no pileup rejection..: "<<lNInelasticEvents<<endl; | |
736 | cout<<" Number of events, this multiplicity....: "<<lNInelasticEvents * ( lNEventsAfterPileup / lNEventsBeforePileup ) <<endl; | |
737 | lNInelasticEvents = lNInelasticEvents * ( lNEventsAfterPileup / lNEventsBeforePileup ); | |
738 | cout<<"--------------------------------------------------------"<<endl; | |
739 | ||
740 | } | |
741 | ||
57a6d361 | 742 | //Variable Definition============================================= |
743 | //Kinematic | |
744 | Float_t lPt, lRap, lPtMC, lNegEta, lPosEta; | |
745 | //Invariant Masses | |
746 | Float_t lInvariantMass, lInvariantMassCompetingOne, lInvariantMassCompetingTwo; //wildcard | |
747 | //DCA Variables | |
748 | Float_t lDcaV0Daughters; | |
749 | Float_t lDcaPosToPrimVertex, lDcaNegToPrimVertex; | |
750 | //Cosine of Pointing Angle variable | |
751 | Float_t lV0CosinePointingAngle; | |
752 | //Decay Radius and distance over total momentum | |
753 | Float_t lV0Radius, lDistOverTotMom; | |
754 | //Least Number of TPC Clusters | |
755 | Int_t lLeastNbrCrossedRows; | |
756 | Float_t lLeastNbrCrossedRowsOverFindable; | |
757 | //TPC dE/dx acquired with AliPIDResponse class | |
758 | Float_t lNSigmasPosProton,lNSigmasNegProton,lNSigmasPosPion,lNSigmasNegPion; | |
d1d5d7e6 | 759 | Float_t lArmPt,lArmAlpha; |
760 | //Multiplicity Variable | |
761 | Int_t lMultiplicity = -1; | |
57a6d361 | 762 | //================================================================ |
763 | ||
764 | //Linking to Tree================================================= | |
765 | ||
766 | //--- Base Variables ---------------------------------------------- | |
767 | lTree->SetBranchAddress("fTreeVariablePosEta",&lPosEta); | |
768 | lTree->SetBranchAddress("fTreeVariableNegEta",&lNegEta); | |
769 | lTree->SetBranchAddress("fTreeVariablePt",&lPt); | |
770 | if ( fWhichParticle == "Lambda" ) lTree->SetBranchAddress("fTreeVariableInvMassLambda",&lInvariantMass); | |
771 | if ( fWhichParticle == "AntiLambda" ) lTree->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMass); | |
772 | if ( fWhichParticle == "K0Short" ) lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMass); | |
773 | if ( fWhichParticle == "Lambda" ){ //For symmetry of computation... | |
774 | lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne); | |
775 | lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo); | |
776 | } | |
777 | if ( fWhichParticle == "AntiLambda" ){ //For symmetry of computation... | |
778 | lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne); | |
779 | lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo); | |
780 | } | |
781 | if ( fWhichParticle == "K0Short" ){ | |
782 | lTree->SetBranchAddress("fTreeVariableInvMassLambda" ,&lInvariantMassCompetingOne); | |
783 | lTree->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMassCompetingTwo); | |
784 | } | |
785 | if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ) | |
786 | lTree->SetBranchAddress("fTreeVariableRapLambda",&lRap); | |
787 | if ( fWhichParticle == "K0Short" ) | |
788 | lTree->SetBranchAddress("fTreeVariableRapK0Short",&lRap); | |
789 | lTree->SetBranchAddress("fTreeVariableDistOverTotMom",&lDistOverTotMom); | |
790 | lTree->SetBranchAddress("fTreeVariableLeastNbrCrossedRows",&lLeastNbrCrossedRows); | |
791 | lTree->SetBranchAddress("fTreeVariableLeastRatioCrossedRowsOverFindable",&lLeastNbrCrossedRowsOverFindable); | |
792 | //--- TPC dEdx Variables ------------------------------------------ | |
793 | lTree->SetBranchAddress("fTreeVariableNSigmasPosProton",&lNSigmasPosProton); | |
794 | lTree->SetBranchAddress("fTreeVariableNSigmasNegProton",&lNSigmasNegProton); | |
795 | lTree->SetBranchAddress("fTreeVariableNSigmasPosPion",&lNSigmasPosPion); | |
796 | lTree->SetBranchAddress("fTreeVariableNSigmasNegPion",&lNSigmasNegPion); | |
797 | //--- Topological selection variables ----------------------------- | |
798 | lTree->SetBranchAddress("fTreeVariableV0Radius",&lV0Radius); | |
799 | lTree->SetBranchAddress("fTreeVariableDcaNegToPrimVertex",&lDcaNegToPrimVertex); | |
800 | lTree->SetBranchAddress("fTreeVariableDcaPosToPrimVertex",&lDcaPosToPrimVertex); | |
801 | lTree->SetBranchAddress("fTreeVariableDcaV0Daughters",&lDcaV0Daughters); | |
802 | lTree->SetBranchAddress("fTreeVariableV0CosineOfPointingAngle",&lV0CosinePointingAngle); | |
d1d5d7e6 | 803 | //--- Multiplicity Variable ---------------------------------------- |
804 | lTree->SetBranchAddress("fTreeVariableMultiplicity",&lMultiplicity); | |
805 | //--- Armenteros-Podolansky ---------------------------------------- | |
806 | lTree->SetBranchAddress("fTreeVariableAlphaV0",&lArmAlpha); | |
807 | lTree->SetBranchAddress("fTreeVariablePtArmV0",&lArmPt); | |
57a6d361 | 808 | //================================================================ |
809 | ||
810 | ||
811 | cout<<"--------------- Real Data File Loop 1 ------------------"<<endl; | |
812 | Long_t lOneTenthOfNCandidates = ((double)(lNCandidates) / 10. ); | |
813 | for(Long_t icand = 0;icand<lNCandidates;icand++){ | |
814 | lTree->GetEntry(icand); | |
815 | if( icand % lOneTenthOfNCandidates == 0 ) | |
816 | cout<<" Currently at candidate........: "<<icand<<" / "<<lNCandidates<<" ( "<<(long)(((double)(icand)/(double)(lNCandidates))*(100.+1e-3))<<"% )"<<endl; | |
817 | if(TMath::Abs(lRap)<fRapidityBoundary && | |
818 | TMath::Abs(lNegEta) <= fCutDaughterEta && | |
819 | TMath::Abs(lPosEta) <= fCutDaughterEta && | |
820 | lV0Radius >= fCutV0Radius && | |
821 | lDcaNegToPrimVertex >= fCutDCANegToPV && | |
822 | lDcaPosToPrimVertex >= fCutDCAPosToPV && | |
823 | lDcaV0Daughters <= fCutDCAV0Daughters && | |
824 | lV0CosinePointingAngle >= fCutV0CosPA && | |
825 | lParticleMass*lDistOverTotMom <= fCutProperLifetime && | |
826 | lLeastNbrCrossedRows >= fCutLeastNumberOfCrossedRows && | |
827 | lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable && | |
828 | TMath::Abs(lInvariantMassCompetingOne - lCompetingParticleMass) > fCutCompetingV0Rejection && | |
829 | TMath::Abs(lInvariantMassCompetingTwo - lCompetingParticleMass) > fCutCompetingV0Rejection && | |
25231937 | 830 | ( fSpecialArmenterosCutK0s == kFALSE || ( fSpecialArmenterosCutK0s == kTRUE && lArmPt*5>TMath::Abs(lArmAlpha) ) ) && |
57a6d361 | 831 | ( //official response code |
832 | ( fWhichParticle == "Lambda" | |
833 | && TMath::Abs(lNSigmasNegPion) <= fCutTPCPIDNSigmas | |
834 | && TMath::Abs(lNSigmasPosProton) <= fCutTPCPIDNSigmas) || | |
835 | ( fWhichParticle == "AntiLambda" | |
836 | && TMath::Abs(lNSigmasPosPion) <= fCutTPCPIDNSigmas | |
837 | && TMath::Abs(lNSigmasNegProton) <= fCutTPCPIDNSigmas) || | |
838 | ( fWhichParticle == "K0Short" | |
839 | && TMath::Abs(lNSigmasPosPion) <= fCutTPCPIDNSigmas | |
840 | && TMath::Abs(lNSigmasNegPion) <= fCutTPCPIDNSigmas) | |
d1d5d7e6 | 841 | ) && |
842 | ( //Multiplicity Switch | |
843 | fPerformMultiplicityStudy == kFALSE || //No mult selection, or... | |
844 | (fPerformMultiplicityStudy == kTRUE && //inside mult bin | |
845 | lMultiplicity>=fLoMultBound && | |
846 | lMultiplicity<=fHiMultBound) | |
57a6d361 | 847 | ) |
848 | ){ | |
849 | lWeAreAtBin = fHistPt->FindBin(lPt)-1; | |
850 | if(lWeAreAtBin == -1) lWeAreAtBin = 99; //UnderFlow, special treatment | |
851 | lHistoFullV0[lWeAreAtBin]->Fill(lInvariantMass); //fill with specific inv mass | |
852 | } | |
853 | } | |
854 | cout<<"--------------- Loop Completed -------------------------"<<endl; | |
855 | cout<<endl; | |
856 | ||
857 | cout<<"--------------- Peak Finding (gauss+linear) ------------"<<endl; | |
858 | //== Variables for holding peak position, width ============== | |
859 | Double_t lPeakPosition[100]; | |
860 | Double_t lPeakWidth[100]; | |
861 | Double_t lLeftBgLeftLimit[100]; | |
862 | Double_t lLeftBgRightLimit[100]; | |
863 | Double_t lPeakLeftLimit[100]; | |
864 | Double_t lPeakRightLimit[100]; | |
865 | Double_t lRightBgLeftLimit[100]; | |
866 | Double_t lRightBgRightLimit[100]; | |
867 | //May be needed if bg areas are different in the future | |
868 | //Double_t lScaleFactor[100]; | |
869 | ||
870 | TLine *lLineLeftMost[100]; | |
871 | TLine *lLineLeft[100]; | |
872 | TLine *lLineRight[100]; | |
873 | TLine *lLineRightMost[100]; | |
874 | ||
875 | TLine *lLineLeftMostMC[100]; | |
876 | TLine *lLineLeftMC[100]; | |
877 | TLine *lLineRightMC[100]; | |
878 | TLine *lLineRightMostMC[100]; | |
879 | ||
d1d5d7e6 | 880 | Double_t lMiddle = 0; |
881 | Double_t lUpperFit = 0; | |
882 | Double_t lLowerFit = 0; | |
883 | ||
57a6d361 | 884 | char fgausname[100]; |
885 | TF1 *fgausPt[100]; | |
886 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++){ | |
887 | cout<<"---> Peak Finding, bin #"<<ibin<<"..."<<endl; | |
888 | sprintf(fgausname,"fGausPt%i",ibin); | |
889 | if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
d1d5d7e6 | 890 | lMiddle = 0.5 * ( fLDataLower->Eval( fHistPt->GetBinCenter(ibin+1) ) + fLDataUpper->Eval( fHistPt->GetBinCenter(ibin+1) ) ); |
891 | lUpperFit = lMiddle + 0.7 * ( fLDataUpper->Eval( fHistPt->GetBinCenter(ibin+1) ) - lMiddle ); | |
892 | lLowerFit = lMiddle - 0.7 * ( lMiddle - fLDataLower->Eval( fHistPt->GetBinCenter(ibin+1) ) ); | |
893 | fgausPt[ibin]= new TF1(fgausname,"[0]*TMath::Gaus(x,[1],[2])+[3]*x+[4]", lLowerFit, lUpperFit ); | |
57a6d361 | 894 | fgausPt[ibin]->SetParameter(1,1.116); |
895 | fgausPt[ibin]->SetParameter(2,0.0025); | |
896 | fgausPt[ibin]->SetParLimits(1,1,1.2); | |
d1d5d7e6 | 897 | fgausPt[ibin]->SetParLimits(2,0.001,0.01); |
57a6d361 | 898 | } |
899 | if ( fWhichParticle == "K0Short"){ | |
900 | fgausPt[ibin]= new TF1(fgausname,"[0]*TMath::Gaus(x,[1],[2])+[3]*x+[4]",0.498-0.06,0.498+0.060); | |
901 | fgausPt[ibin]->SetParameter(1,0.498); | |
902 | fgausPt[ibin]->SetParameter(2,0.004); | |
903 | } | |
904 | fgausPt[ibin]->SetParameter(0,lHistoFullV0[ibin]->GetMaximum() * 0.9); | |
905 | ||
906 | ||
907 | fgausPt[ibin]->SetParameter(3,0); | |
908 | fgausPt[ibin]->SetParameter(4,lHistoFullV0[ibin]->GetMaximum() * 0.1); | |
909 | lHistoFullV0[ibin]->Fit(fgausname,"QREM0"); | |
910 | lPeakPosition[ibin] = fgausPt[ibin]->GetParameter(1); | |
d1d5d7e6 | 911 | lPeakWidth[ibin] = TMath::Abs( fgausPt[ibin]->GetParameter(2) ); |
57a6d361 | 912 | cout<<"---> ["<<fptbinlimits[ibin]<<" - "<<fptbinlimits[ibin+1]<<" GeV/c]\tPeak at: "<<lPeakPosition[ibin]<<", sigma = "<<lPeakWidth[ibin]<<endl; |
913 | //Find Corresponding Limits In this bin | |
914 | lLeftBgLeftLimit[ibin] = lPeakPosition[ibin] - 2.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin]; | |
915 | lLeftBgRightLimit[ibin] = lPeakPosition[ibin] - 1.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin]; | |
916 | lPeakLeftLimit[ibin] = lPeakPosition[ibin] - 1.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin]; | |
917 | lPeakRightLimit[ibin] = lPeakPosition[ibin] + 1.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin]; | |
918 | lRightBgLeftLimit[ibin] = lPeakPosition[ibin] + 1.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin]; | |
919 | lRightBgRightLimit[ibin]= lPeakPosition[ibin] + 2.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin]; | |
920 | //lScaleFactor[ibin] = (lPeakRightLimit[ibin] - lPeakLeftLimit[ibin]) / ( (lLeftBgRightLimit[ibin]-lLeftBgLeftLimit[ibin]) + (lRightBgRightLimit[ibin] - lRightBgLeftLimit[ibin]) ); | |
921 | fHistPeakPosition->SetBinContent(ibin+1, lPeakPosition[ibin]); | |
922 | fHistPeakPosition->SetBinError(ibin+1, lPeakWidth[ibin]); | |
923 | //Create Signal Extraction Range Histogram | |
924 | fHistSignalExtractionRange->SetBinContent(ibin+1, lPeakPosition[ibin]); | |
925 | fHistSignalExtractionRange->SetBinError(ibin+1, 2.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin] ); | |
926 | //Create appropriate TLine Objects for Canvases | |
927 | lLineLeftMost[ibin] = new TLine( lLeftBgLeftLimit[ibin], 0, lLeftBgLeftLimit[ibin], lHistoFullV0[ibin]->GetMaximum() * 0.95 ); | |
928 | lLineLeft[ibin] = new TLine( lLeftBgRightLimit[ibin], 0, lLeftBgRightLimit[ibin], lHistoFullV0[ibin]->GetMaximum() * 0.95 ); | |
929 | lLineRight[ibin] = new TLine( lRightBgLeftLimit[ibin], 0, lRightBgLeftLimit[ibin], lHistoFullV0[ibin]->GetMaximum() * 0.95 ); | |
930 | lLineRightMost[ibin] = new TLine( lRightBgRightLimit[ibin], 0, lRightBgRightLimit[ibin], lHistoFullV0[ibin]->GetMaximum() * 0.95 ); | |
931 | ||
932 | //Preparing Canvas for storing... | |
933 | lCanvasHistoFullV0[ibin]->cd(); | |
934 | lHistoFullV0[ibin]->Draw(); | |
935 | lLineLeftMost[ibin]->Draw(); | |
936 | lLineLeft[ibin]->Draw(); | |
937 | lLineRight[ibin]->Draw(); | |
938 | lLineRightMost[ibin]->Draw(); | |
939 | fgausPt[ibin]->Draw("same"); | |
940 | ||
941 | } | |
942 | cout<<"--------------- Peak Finding Finished ------------------"<<endl; | |
943 | ||
944 | //Defining Signal holding variables=============================== | |
945 | Double_t lSigRealV0[100]; | |
946 | Double_t lSigErrRealV0[100]; | |
947 | Double_t lSigMCV0[100]; | |
948 | Double_t lSigErrMCV0[100]; | |
949 | //================================================================ | |
950 | Double_t lLeftPlusRightBgV0[100]; | |
951 | Double_t lSigPlusCenterBgV0[100]; | |
952 | Double_t lLeftPlusRightBgV0MC[100]; | |
953 | Double_t lSigPlusCenterBgV0MC[100]; | |
954 | for(Long_t n1 = 0; n1<100; n1++){ | |
955 | lLeftPlusRightBgV0[n1]=0; | |
956 | lSigPlusCenterBgV0[n1]=0; | |
957 | lLeftPlusRightBgV0MC[n1]=0; | |
958 | lSigPlusCenterBgV0MC[n1]=0; | |
959 | } | |
960 | //================================================================ | |
961 | cout<<endl; | |
962 | cout<<"--------------- Real Data File Loop 2 ------------------"<<endl; | |
963 | for(Long_t icand = 0;icand<lNCandidates;icand++){ | |
964 | lTree->GetEntry(icand); | |
965 | if( icand % lOneTenthOfNCandidates == 0 ) | |
966 | cout<<" Currently at candidate........: "<<icand<<" / "<<lNCandidates<<" ( "<<(long)(((double)(icand)/(double)(lNCandidates))*(100.+1e-3))<<"% )"<<endl; | |
967 | if(TMath::Abs(lRap)<fRapidityBoundary && | |
968 | TMath::Abs(lNegEta) <= fCutDaughterEta && | |
969 | TMath::Abs(lPosEta) <= fCutDaughterEta && | |
970 | lV0Radius >= fCutV0Radius && | |
971 | lDcaNegToPrimVertex >= fCutDCANegToPV && | |
972 | lDcaPosToPrimVertex >= fCutDCAPosToPV && | |
973 | lDcaV0Daughters <= fCutDCAV0Daughters && | |
974 | lV0CosinePointingAngle >= fCutV0CosPA && | |
975 | lParticleMass*lDistOverTotMom <= fCutProperLifetime && | |
976 | lLeastNbrCrossedRows >= fCutLeastNumberOfCrossedRows && | |
977 | lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable && | |
978 | TMath::Abs(lInvariantMassCompetingOne - lCompetingParticleMass) > fCutCompetingV0Rejection && | |
979 | TMath::Abs(lInvariantMassCompetingTwo - lCompetingParticleMass) > fCutCompetingV0Rejection && | |
25231937 | 980 | ( fSpecialArmenterosCutK0s == kFALSE || ( fSpecialArmenterosCutK0s == kTRUE && lArmPt*5>TMath::Abs(lArmAlpha) ) ) && |
57a6d361 | 981 | ( //official response code |
982 | ( fWhichParticle == "Lambda" | |
983 | && TMath::Abs(lNSigmasNegPion) <= fCutTPCPIDNSigmas | |
984 | && TMath::Abs(lNSigmasPosProton) <= fCutTPCPIDNSigmas) || | |
985 | ( fWhichParticle == "AntiLambda" | |
986 | && TMath::Abs(lNSigmasPosPion) <= fCutTPCPIDNSigmas | |
987 | && TMath::Abs(lNSigmasNegProton) <= fCutTPCPIDNSigmas) || | |
988 | ( fWhichParticle == "K0Short" | |
989 | && TMath::Abs(lNSigmasPosPion) <= fCutTPCPIDNSigmas | |
990 | && TMath::Abs(lNSigmasNegPion) <= fCutTPCPIDNSigmas) | |
d1d5d7e6 | 991 | ) && |
992 | ( //Multiplicity Switch | |
993 | fPerformMultiplicityStudy == kFALSE || //No mult selection, or... | |
994 | (fPerformMultiplicityStudy == kTRUE && //inside mult bin | |
995 | lMultiplicity>=fLoMultBound && | |
996 | lMultiplicity<=fHiMultBound) | |
57a6d361 | 997 | ) |
998 | ){ // Start Entry Loop | |
999 | lWeAreAtBin = fHistPt->FindBin(lPt)-1; | |
1000 | if(lWeAreAtBin == -1) lWeAreAtBin = 99; //UnderFlow, special treatment | |
1001 | ||
1002 | //Extract left and right background areas and peak | |
1003 | //--- Left Background Sample | |
1004 | if(lInvariantMass>lLeftBgLeftLimit[lWeAreAtBin] && lInvariantMass<lLeftBgRightLimit[lWeAreAtBin] ){ lLeftPlusRightBgV0[lWeAreAtBin]++; } | |
1005 | if(lInvariantMass>lRightBgLeftLimit[lWeAreAtBin] && lInvariantMass<lRightBgRightLimit[lWeAreAtBin] ){ lLeftPlusRightBgV0[lWeAreAtBin]++; } | |
1006 | //--- Peak Region | |
1007 | if(lInvariantMass>lPeakLeftLimit[lWeAreAtBin] && lInvariantMass<lPeakRightLimit[lWeAreAtBin] ){ lSigPlusCenterBgV0[lWeAreAtBin]++; } | |
1008 | } // End Entry Loop | |
1009 | } | |
1010 | cout<<"--------------- Loop Completed -------------------------"<<endl; | |
1011 | cout<<endl; | |
1012 | cout<<"--------------- Memory Cleanup -------------------------"<<endl; | |
1013 | if ( v0list ) { | |
1014 | v0list->Delete(); | |
1015 | delete v0list; | |
1016 | } | |
1017 | ||
1018 | file->Close("R"); | |
1019 | file->Delete(); | |
1020 | delete file; | |
1021 | cout<<endl; | |
1022 | ||
1023 | TF1 *lfitNoise[50]; | |
1024 | TF1 *lSampleNoise[50]; | |
1025 | char lFitNameOne[50]; | |
1026 | ||
1027 | if( fFitBackgroundSwitch ){ | |
1028 | cout<<"--------------- Backgrounds: fitting with linear -------"<<endl; | |
1029 | for(long i=0; i<fptbinnumb; i++){ | |
1030 | //Define Function to Fit Background | |
1031 | sprintf(lFitNameOne,"lfitNoise%i",((int)(i))); | |
1032 | //cout<<"creating fitnoise, named "<<lFitNameOne<<endl; | |
1033 | lfitNoise[i] = new TF1(lFitNameOne, this, &AliV0Module::MyBgPol1, | |
1034 | RoundToThousandth ( lLeftBgLeftLimit[i] ), | |
1035 | RoundToThousandth ( lRightBgRightLimit[i] ), 4 , "AliV0Module", "MyBgPol1"); | |
1036 | lfitNoise[i] -> FixParameter(2,RoundToThousandth ( lLeftBgRightLimit[i] ) ); | |
1037 | lfitNoise[i] -> FixParameter(3,RoundToThousandth ( lRightBgLeftLimit[i] ) ); | |
1038 | lfitNoise[i] -> SetParameter(0,lLeftPlusRightBgV0[i] * lHistoFullV0[i]->GetBinWidth(5) / (lRightBgLeftLimit[i]-lLeftBgRightLimit[i] + 1e-6 ) ); | |
1039 | lfitNoise[i] -> SetParameter(1,0); | |
1040 | cout<<"Guessed Parameter 0 fot "<<lFitNameOne<<" to be "<<lLeftPlusRightBgV0[i] * lHistoFullV0[i]->GetBinWidth(5) / (lRightBgLeftLimit[i]-lLeftBgRightLimit[i] + 1e-6 )<<endl; | |
1041 | sprintf(lFitNameOne,"lSampleNoise%i",((int)(i))); | |
1042 | ||
1043 | //Define Function to Sample Background | |
1044 | //cout<<"creating sample "<<i<<endl; | |
1045 | lSampleNoise[i] = new TF1(lFitNameOne, this, &AliV0Module::MyBgPolToEval1, | |
1046 | RoundToThousandth ( lLeftBgLeftLimit[i] ), | |
1047 | RoundToThousandth ( lRightBgRightLimit[i] ), 2 , "AliV0Module", "MyBgPolToEval1"); | |
1048 | } | |
1049 | for(long i=0; i<fptbinnumb; i++){ | |
1050 | //cout<<"Fitting function for bin "<<i<<", get name = "<<lfitNoise[i]->GetName()<<endl; | |
1051 | sprintf(lFitNameOne,"lfitNoise%i",((int)(i))); | |
1052 | lHistoFullV0[i] -> Fit(lFitNameOne,"LLrie+0"); | |
1053 | lSampleNoise[i]->SetParameter(0, lfitNoise[i]->GetParameter(0) ); | |
1054 | lSampleNoise[i]->SetParameter(1, lfitNoise[i]->GetParameter(1) ); | |
1055 | } | |
1056 | for(long i=0; i<fptbinnumb; i++){ | |
1057 | cout<<"Overriding Background info: Was "<<lLeftPlusRightBgV0[i]<<", is now "<<lSampleNoise[i]->Integral( lPeakLeftLimit[i], lPeakRightLimit[i] )/lHistoFullV0[i]->GetBinWidth(5)<<endl; | |
1058 | lLeftPlusRightBgV0[i] = lSampleNoise[i]->Integral( lPeakLeftLimit[i], lPeakRightLimit[i] )/lHistoFullV0[i]->GetBinWidth(5); | |
1059 | } | |
1060 | cout<<"--------------- Fitting Finished! ----------------------"<<endl; | |
1061 | } | |
1062 | ||
1063 | //============================================================= | |
1064 | // Compute Signal + Sig to noise | |
1065 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1066 | //Signal computation | |
1067 | lSigRealV0[ipoint] = lSigPlusCenterBgV0[ipoint] - lLeftPlusRightBgV0[ipoint]; | |
1068 | lSigErrRealV0[ipoint] = TMath::Sqrt(lSigPlusCenterBgV0[ipoint]+lLeftPlusRightBgV0[ipoint]); | |
1069 | //Error Formula: Equivalent to Sqrt(S+B+B) = Sqrt(S+2B) | |
1070 | ||
1071 | //Signal-to-noise computation | |
1072 | if( lLeftPlusRightBgV0[ipoint] != 0 ){ | |
1073 | fHistSigToNoise->SetBinContent(ipoint+1, (lSigPlusCenterBgV0[ipoint] - lLeftPlusRightBgV0[ipoint])/lLeftPlusRightBgV0[ipoint] ); | |
1074 | }else{ | |
1075 | fHistSigToNoise->SetBinContent(ipoint+1, -1) ; //-1 means: no background | |
1076 | } | |
1077 | } | |
1078 | //============================================================= | |
1079 | ||
1080 | cout<<"--------------- Extracted Signal -----------------------"<<endl; | |
1081 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1082 | cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tSignal: "<<lSigRealV0[ipoint]<<" +/- "<<lSigErrRealV0[ipoint]<<endl; | |
1083 | } | |
1084 | cout<<"--------------------------------------------------------"<<endl; | |
1085 | ||
1086 | //============================================================= | |
1087 | // Preparations for MC loop | |
1088 | //Extra info available only in MC | |
1089 | Int_t lPID = 0; | |
1090 | Float_t lNegTransvMomentumMC, lPosTransvMomentumMC; | |
1091 | Int_t lPIDPositive = 0; | |
1092 | Int_t lPIDNegative = 0; | |
1093 | Int_t lPIDMother = -1; Float_t lPtMother = -1; | |
1094 | Int_t lPrimaryStatus = 0; | |
1095 | Int_t lPrimaryStatusMother = 0; | |
1096 | //============================================================= | |
1097 | cout<<endl; | |
1098 | //Prepare containers for Geant3/fluka | |
1099 | //Will only be used and saved if we're dealing with Lambdas | |
1100 | TH1F *lProtonMomentum[100]; | |
1101 | char lNameOne[100]; | |
1102 | for(Int_t ibin=0;ibin<100;ibin++){ | |
1103 | sprintf(lNameOne,"lProtonMomentumBin%i",ibin); | |
1104 | lProtonMomentum[ibin] = new TH1F(lNameOne,"",800,0,20); | |
1105 | if(fWhichParticle == "Lambda") bindescription="#Lambda, bin #"; | |
1106 | if(fWhichParticle == "AntiLambda") bindescription="#bar{#Lambda}, bin #"; | |
1107 | if(fWhichParticle == "K0Short") bindescription="K^{0}_{S}, bin #"; | |
1108 | bindescription.Append(IntToString( ibin )); | |
1109 | if ( ibin < fptbinnumb ){ | |
1110 | bindescription.Append(", "); | |
1111 | bindescription.Append(DoubleToString(fptbinlimits[ ibin ])); | |
1112 | bindescription.Append("-"); | |
1113 | bindescription.Append(DoubleToString(fptbinlimits[ ibin+1 ])); | |
1114 | bindescription.Append("GeV/c"); | |
1115 | } | |
1116 | lProtonMomentum[ibin]->SetTitle(bindescription); | |
1117 | } | |
1118 | ||
1119 | //MC File Acquisition============================================= | |
1120 | cout<<"--------------- Opening MC file ------------------------"<<endl; | |
1121 | TFile* fileMC = TFile::Open(fMCDataFile, "READ"); | |
1122 | fileMC->cd("PWG2CheckPerformanceLambda_PP_MC"); | |
1123 | TList* v0listMC = (TList*)fileMC->FindObjectAny("clistV0MC"); | |
1124 | TTree *lTreeMC; | |
1125 | TH1F* fHistMultiplicityBeforeTrigSelMC; | |
1126 | lTreeMC = (TTree*)fileMC->FindObjectAny("fTree"); | |
1127 | fHistMultiplicityBeforeTrigSelMC = (TH1F*)v0listMC->FindObject("fHistMultiplicityBeforeTrigSel"); | |
1128 | ||
1129 | Long_t lNCandidatesMC = lTreeMC->GetEntries(); | |
1130 | Long_t lNEventsMC = fHistMultiplicityBeforeTrigSelMC->GetEntries(); | |
1131 | cout<<" MC Events (before trig sel)...: "<<lNEventsMC<<endl; | |
1132 | cout<<" Cascade MC Candidates.........: "<<lNCandidatesMC<<endl; | |
1133 | cout<<"--------------------------------------------------------"<<endl; | |
1134 | //================================================================ | |
1135 | ||
1136 | //Linking to Tree================================================= | |
1137 | //--- Base Variables ---------------------------------------------- | |
1138 | lTreeMC->SetBranchAddress("fTreeVariablePosEta",&lPosEta); | |
1139 | lTreeMC->SetBranchAddress("fTreeVariableNegEta",&lNegEta); | |
1140 | lTreeMC->SetBranchAddress("fTreeVariablePrimaryStatus",&lPrimaryStatus); | |
1141 | lTreeMC->SetBranchAddress("fTreeVariablePrimaryStatusMother",&lPrimaryStatusMother); | |
1142 | lTreeMC->SetBranchAddress("fTreeVariablePt",&lPt); | |
1143 | lTreeMC->SetBranchAddress("fTreeVariablePtXiMother",&lPtMother); | |
1144 | lTreeMC->SetBranchAddress("fTreeVariablePtMC",&lPtMC); | |
1145 | if ( fWhichParticle == "Lambda" ) lTreeMC->SetBranchAddress("fTreeVariableInvMassLambda",&lInvariantMass); | |
1146 | if ( fWhichParticle == "AntiLambda" ) lTreeMC->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMass); | |
1147 | if ( fWhichParticle == "K0Short" ) lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMass); | |
1148 | if ( fWhichParticle == "Lambda" ){ //For symmetry of computation... | |
1149 | lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne); | |
1150 | lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo); | |
1151 | } | |
1152 | if ( fWhichParticle == "AntiLambda" ){ //For symmetry of computation... | |
1153 | lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne); | |
1154 | lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo); | |
1155 | } | |
1156 | if ( fWhichParticle == "K0Short" ){ | |
1157 | lTreeMC->SetBranchAddress("fTreeVariableInvMassLambda" ,&lInvariantMassCompetingOne); | |
1158 | lTreeMC->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMassCompetingTwo); | |
1159 | } | |
1160 | //if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ) | |
1161 | // lTreeMC->SetBranchAddress("fTreeVariableRapLambda",&lRap); | |
1162 | //if ( fWhichParticle == "K0Short" ) | |
1163 | // lTreeMC->SetBranchAddress("fTreeVariableRapK0Short",&lRap); | |
1164 | ||
1165 | lTreeMC->SetBranchAddress("fTreeVariableRapMC",&lRap); | |
1166 | lTreeMC->SetBranchAddress("fTreeVariableLeastNbrCrossedRows",&lLeastNbrCrossedRows); | |
1167 | lTreeMC->SetBranchAddress("fTreeVariableLeastRatioCrossedRowsOverFindable",&lLeastNbrCrossedRowsOverFindable); | |
1168 | //--- Topological selection variables ----------------------------- | |
1169 | lTreeMC->SetBranchAddress("fTreeVariableV0Radius",&lV0Radius); | |
1170 | lTreeMC->SetBranchAddress("fTreeVariableDcaNegToPrimVertex",&lDcaNegToPrimVertex); | |
1171 | lTreeMC->SetBranchAddress("fTreeVariableDcaPosToPrimVertex",&lDcaPosToPrimVertex); | |
1172 | lTreeMC->SetBranchAddress("fTreeVariableDcaV0Daughters",&lDcaV0Daughters); | |
1173 | lTreeMC->SetBranchAddress("fTreeVariableV0CosineOfPointingAngle",&lV0CosinePointingAngle); | |
1174 | lTreeMC->SetBranchAddress("fTreeVariableDistOverTotMom",&lDistOverTotMom); | |
1175 | //---- Extra Information Available in MC only -------------------- | |
1176 | lTreeMC->SetBranchAddress("fTreeVariableNegTransvMomentumMC",&lNegTransvMomentumMC); | |
1177 | lTreeMC->SetBranchAddress("fTreeVariablePosTransvMomentumMC",&lPosTransvMomentumMC); | |
1178 | lTreeMC->SetBranchAddress("fTreeVariablePID",&lPID); | |
1179 | lTreeMC->SetBranchAddress("fTreeVariablePIDMother",&lPIDMother); | |
1180 | lTreeMC->SetBranchAddress("fTreeVariablePIDPositive",&lPIDPositive); | |
1181 | lTreeMC->SetBranchAddress("fTreeVariablePIDNegative",&lPIDNegative); | |
d1d5d7e6 | 1182 | //--- Multiplicity ------------------------------------------------ |
1183 | lTreeMC->SetBranchAddress("fTreeVariableMultiplicity",&lMultiplicity); | |
1184 | //--- Armenteros-Podolansky ---------------------------------------- | |
1185 | lTreeMC->SetBranchAddress("fTreeVariableAlphaV0",&lArmAlpha); | |
1186 | lTreeMC->SetBranchAddress("fTreeVariablePtArmV0",&lArmPt); | |
57a6d361 | 1187 | //================================================================ |
1188 | ||
1189 | //================================================================ | |
1190 | cout<<endl; | |
1191 | cout<<"--------------- MC Data File Loop ---------------------"<<endl; | |
1192 | Long_t lOneTenthOfNCandidatesMC = ((double)(lNCandidatesMC) / 10. ); | |
1193 | for(Long_t icand = 0;icand<lNCandidatesMC;icand++){ | |
1194 | lTreeMC->GetEntry(icand); | |
1195 | if( icand % lOneTenthOfNCandidatesMC == 0 ) | |
1196 | cout<<" Currently at candidate........: "<<icand<<" / "<<lNCandidatesMC<<" ( "<<(long)(((double)(icand)/(double)(lNCandidatesMC))*(100.+1e-3))<<"% )"<<endl; | |
1197 | if(TMath::Abs(lRap)<fRapidityBoundary && | |
1198 | TMath::Abs(lNegEta) <= fCutDaughterEta && | |
1199 | TMath::Abs(lPosEta) <= fCutDaughterEta && | |
1200 | lV0Radius >= fCutV0Radius && | |
1201 | lDcaNegToPrimVertex >= fCutDCANegToPV && | |
1202 | lDcaPosToPrimVertex >= fCutDCAPosToPV && | |
1203 | lDcaV0Daughters <= fCutDCAV0Daughters && | |
1204 | lV0CosinePointingAngle >= fCutV0CosPA && | |
1205 | lParticleMass*lDistOverTotMom <= fCutProperLifetime && | |
1206 | lLeastNbrCrossedRows >= fCutLeastNumberOfCrossedRows && | |
1207 | lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable && | |
1208 | TMath::Abs(lInvariantMassCompetingOne - lCompetingParticleMass) > fCutCompetingV0Rejection && | |
1209 | TMath::Abs(lInvariantMassCompetingTwo - lCompetingParticleMass) > fCutCompetingV0Rejection && | |
25231937 | 1210 | ( fSpecialArmenterosCutK0s == kFALSE || ( fSpecialArmenterosCutK0s == kTRUE && lArmPt*5>TMath::Abs(lArmAlpha) ) ) && |
57a6d361 | 1211 | ( //perfect PID association, IsPhysicalPrimary association |
1212 | ( fWhichParticle == "Lambda" | |
1213 | && lPID == 3122 //V0 is a Lambda | |
1214 | && lPIDPositive == 2212 //Pos Daughter is p | |
1215 | && lPIDNegative == -211 //Neg Daughter is pi- | |
1216 | && lPrimaryStatus == 1 | |
1217 | ) || | |
1218 | ( fWhichParticle == "AntiLambda" | |
1219 | && lPID == -3122 //V0 is an AntiLambda | |
1220 | && lPIDPositive == 211 //Pos Daughter is pi+ | |
1221 | && lPIDNegative == -2212 //Neg Daughter is antiproton | |
1222 | && lPrimaryStatus == 1 | |
1223 | ) || | |
1224 | ( fWhichParticle == "K0Short" | |
1225 | && lPID == 310 //V0 is an AntiLambda | |
1226 | && lPIDPositive == 211 //Pos Daughter is pi+ | |
1227 | && lPIDNegative == -211 //Neg Daughter is pi- | |
1228 | && lPrimaryStatus == 1 //K0Short is a primary | |
1229 | ) | |
d1d5d7e6 | 1230 | ) && |
1231 | ( //Multiplicity Switch | |
1232 | fPerformMultiplicityStudy == kFALSE || //No mult selection, or... | |
1233 | (fPerformMultiplicityStudy == kTRUE && //inside mult bin | |
1234 | lMultiplicity>=fLoMultBound && | |
1235 | lMultiplicity<=fHiMultBound) | |
57a6d361 | 1236 | ) |
1237 | ){ // Start Entry Loop | |
1238 | lWeAreAtBin = fHistPt->FindBin(lPtMC)-1; | |
1239 | if(lWeAreAtBin == -1) lWeAreAtBin = 99; //UnderFlow, special treatment | |
1240 | lHistoFullV0MC[lWeAreAtBin]->Fill(lInvariantMass); //fill with specific inv mass | |
d1d5d7e6 | 1241 | f2dHistPtResolution->Fill(lPt,lPtMC); |
57a6d361 | 1242 | //Extract left and right background areas and peak |
1243 | //--- Left Background Sample | |
1244 | if(lInvariantMass>lLeftBgLeftLimit[lWeAreAtBin] && lInvariantMass<lLeftBgRightLimit[lWeAreAtBin] ){ lLeftPlusRightBgV0MC[lWeAreAtBin]++; } | |
1245 | if(lInvariantMass>lRightBgLeftLimit[lWeAreAtBin] && lInvariantMass<lRightBgRightLimit[lWeAreAtBin] ){ lLeftPlusRightBgV0MC[lWeAreAtBin]++; } | |
1246 | //--- Peak Region | |
1247 | if(lInvariantMass>lPeakLeftLimit[lWeAreAtBin] && lInvariantMass<lPeakRightLimit[lWeAreAtBin] ){ lSigPlusCenterBgV0MC[lWeAreAtBin]++; } | |
1248 | //--- Info may be needed for geant3/fluka | |
1249 | if ( fWhichParticle == "Lambda" ) lProtonMomentum[lWeAreAtBin]->Fill( lPosTransvMomentumMC ); | |
1250 | if ( fWhichParticle == "AntiLambda" ) lProtonMomentum[lWeAreAtBin]->Fill( lNegTransvMomentumMC ); | |
d1d5d7e6 | 1251 | //--- Resolution tests |
1252 | lHistResolution[lWeAreAtBin]->Fill(lPt - lPtMC); | |
57a6d361 | 1253 | } // End Entry Loop |
1254 | } | |
1255 | cout<<"--------------- Loop Completed -------------------------"<<endl; | |
1256 | cout<<endl; | |
1257 | ||
1258 | cout<<"--------------- X-Check Peak Finding (gauss+linear) ----"<<endl; | |
1259 | //== Variables for holding peak position, width ============== | |
1260 | Double_t lPeakPositionMC[100]; | |
1261 | Double_t lPeakWidthMC[100]; | |
1262 | ||
1263 | char fgausnameMC[100]; | |
1264 | TF1 *fgausPtMC[100]; | |
1265 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++){ | |
1266 | cout<<"---> Peak Finding, bin #"<<ibin<<" (perfect count = "<<lHistoFullV0MC[ibin]->GetEntries()<<")..."<<endl; | |
1267 | sprintf(fgausnameMC,"fGausPtMC%i",ibin); | |
1268 | if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
1269 | fgausPtMC[ibin]= new TF1(fgausnameMC,"[0]*TMath::Gaus(x,[1],[2])",1.116-0.040,1.116+0.040); | |
1270 | fgausPtMC[ibin]->SetParameter(1,1.116); | |
1271 | fgausPtMC[ibin]->SetParameter(2,0.0025); | |
1272 | fgausPtMC[ibin]->SetParLimits(1,1.1,1.2); | |
1273 | fgausPtMC[ibin]->SetParLimits(2,0.001,0.02); | |
1274 | } | |
1275 | if ( fWhichParticle == "K0Short"){ | |
1276 | fgausPtMC[ibin]= new TF1(fgausnameMC,"[0]*TMath::Gaus(x,[1],[2])",0.498-0.06,0.498+0.060); | |
1277 | fgausPtMC[ibin]->SetParameter(1,0.498); | |
1278 | fgausPtMC[ibin]->SetParameter(2,0.004); | |
1279 | } | |
1280 | fgausPtMC[ibin]->SetParameter(0,lHistoFullV0MC[ibin]->GetMaximum() * 0.9); | |
1281 | //fgausPtMC[ibin]->SetParameter(3,0); | |
1282 | //fgausPtMC[ibin]->SetParameter(4,HistoFullV0MC[ibin]->GetMaximum() * 0.1); | |
1283 | lHistoFullV0MC[ibin]->Fit(fgausnameMC,"QREM0"); | |
d1d5d7e6 | 1284 | lPeakPositionMC[ibin] = TMath::Abs( fgausPtMC[ibin]->GetParameter(1) ); |
57a6d361 | 1285 | lPeakWidthMC[ibin] = fgausPtMC[ibin]->GetParameter(2); |
1286 | cout<<"---> ["<<fptbinlimits[ibin]<<" - "<<fptbinlimits[ibin+1]<<" GeV/c]\tPeak at: "<<lPeakPositionMC[ibin]<<", sigma = "<<lPeakWidthMC[ibin]<<endl; | |
1287 | fHistPeakPositionMC->SetBinContent(ibin+1, lPeakPositionMC[ibin]); | |
1288 | fHistPeakPositionMC->SetBinError(ibin+1, lPeakWidthMC[ibin]); | |
1289 | ||
1290 | //Create appropriate TLine Objects for Canvases | |
1291 | lLineLeftMostMC[ibin] = new TLine( lLeftBgLeftLimit[ibin], 0, lLeftBgLeftLimit[ibin], lHistoFullV0MC[ibin]->GetMaximum() * 0.95 ); | |
1292 | lLineLeftMC[ibin] = new TLine( lLeftBgRightLimit[ibin], 0, lLeftBgRightLimit[ibin], lHistoFullV0MC[ibin]->GetMaximum() * 0.95 ); | |
1293 | lLineRightMC[ibin] = new TLine( lRightBgLeftLimit[ibin], 0, lRightBgLeftLimit[ibin], lHistoFullV0MC[ibin]->GetMaximum() * 0.95 ); | |
1294 | lLineRightMostMC[ibin] = new TLine( lRightBgRightLimit[ibin], 0, lRightBgRightLimit[ibin], lHistoFullV0MC[ibin]->GetMaximum() * 0.95 ); | |
1295 | ||
1296 | //Preparing Canvas for storing... | |
1297 | lCanvasHistoFullV0MC[ibin]->cd(); | |
1298 | lHistoFullV0MC[ibin]->Draw(); | |
1299 | lLineLeftMostMC[ibin]->Draw(); | |
1300 | lLineLeftMC[ibin]->Draw(); | |
1301 | lLineRightMC[ibin]->Draw(); | |
1302 | lLineRightMostMC[ibin]->Draw(); | |
1303 | fgausPtMC[ibin]->Draw("same"); | |
1304 | ||
1305 | } | |
1306 | cout<<"--------------- Peak Finding Finished (in MC) ----------"<<endl; | |
1307 | cout<<endl; | |
1308 | ||
1309 | TF1 *lfitNoiseMC[50]; | |
1310 | TF1 *lSampleNoiseMC[50]; | |
1311 | char lFitNameOneMC[50]; | |
1312 | ||
1313 | if( fFitBackgroundSwitch ){ | |
1314 | cout<<"--------------- Backgrounds: fitting with linear -------"<<endl; | |
1315 | for(long i=0; i<fptbinnumb; i++){ | |
1316 | //Define Function to Fit Background | |
1317 | sprintf(lFitNameOneMC,"lfitNoiseMC%i",(int)i); | |
1318 | cout<<"creating fitnoise, named "<<lFitNameOneMC<<endl; | |
1319 | lfitNoiseMC[i] = new TF1(lFitNameOneMC, this, &AliV0Module::MyBgPol1, | |
1320 | RoundToThousandth ( lLeftBgLeftLimit[i] ), | |
1321 | RoundToThousandth ( lRightBgRightLimit[i] ), 4 , "AliV0Module", "MyBgPol1"); | |
1322 | lfitNoiseMC[i] -> FixParameter(2,RoundToThousandth ( lLeftBgRightLimit[i] ) ); | |
1323 | lfitNoiseMC[i] -> FixParameter(3,RoundToThousandth ( lRightBgLeftLimit[i] ) ); | |
1324 | lfitNoiseMC[i] -> SetParameter(0,lLeftPlusRightBgV0MC[i]*lHistoFullV0MC[i]->GetMaximum() / (lSigPlusCenterBgV0MC[i]+1e-6) ); | |
1325 | lfitNoiseMC[i] -> SetParameter(1,0); | |
1326 | sprintf(lFitNameOneMC,"lSampleNoiseMC%i",(int)i); | |
1327 | ||
1328 | //Define Function to Sample Background | |
1329 | cout<<"creating sample "<<i<<endl; | |
1330 | lSampleNoiseMC[i] = new TF1(lFitNameOneMC, this, &AliV0Module::MyBgPolToEval1, | |
1331 | RoundToThousandth ( lLeftBgLeftLimit[i] ), | |
1332 | RoundToThousandth ( lRightBgRightLimit[i] ), 2, "AliV0Module", "MyBgPolToEval1" ); | |
1333 | } | |
1334 | for(long i=0; i<fptbinnumb; i++){ | |
1335 | cout<<"Fitting function for bin "<<i<<", get name = "<<lfitNoiseMC[i]->GetName()<<endl; | |
1336 | sprintf(lFitNameOneMC,"lfitNoiseMC%i",(int)i); | |
1337 | lHistoFullV0MC[i] -> Fit(lFitNameOneMC,"LLrie+0"); | |
1338 | lSampleNoiseMC[i]->SetParameter(0, lfitNoiseMC[i]->GetParameter(0) ); | |
1339 | lSampleNoiseMC[i]->SetParameter(1, lfitNoiseMC[i]->GetParameter(1) ); | |
1340 | } | |
1341 | for(long i=0; i<fptbinnumb; i++){ | |
1342 | cout<<"Overriding Background info: Was "<<lLeftPlusRightBgV0MC[i]<<", is now "<<lSampleNoiseMC[i]->Integral( lPeakLeftLimit[i], lPeakRightLimit[i] )/lHistoFullV0MC[i]->GetBinWidth(5)<<endl; | |
1343 | lLeftPlusRightBgV0MC[i] = lSampleNoiseMC[i]->Integral( lPeakLeftLimit[i], lPeakRightLimit[i] )/lHistoFullV0MC[i]->GetBinWidth(5); | |
1344 | } | |
1345 | cout<<"--------------- Fitting Finished! ----------------------"<<endl; | |
1346 | } | |
1347 | ||
1348 | ||
1349 | //============================================================= | |
1350 | // Compute Signal | |
1351 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1352 | lSigMCV0[ipoint] = lSigPlusCenterBgV0MC[ipoint] - lLeftPlusRightBgV0MC[ipoint]; | |
1353 | lSigErrMCV0[ipoint] = TMath::Sqrt(lSigPlusCenterBgV0MC[ipoint]+lLeftPlusRightBgV0MC[ipoint]); | |
1354 | //Error Formula: Equivalent to Sqrt(S+B+B) = Sqrt(S+2B) | |
1355 | ||
1356 | //Signal-to-noise computation | |
1357 | if( lLeftPlusRightBgV0MC[ipoint] != 0 ){ | |
1358 | fHistSigToNoiseMC->SetBinContent(ipoint+1, (lSigPlusCenterBgV0MC[ipoint] - lLeftPlusRightBgV0MC[ipoint])/lLeftPlusRightBgV0MC[ipoint] ); | |
1359 | }else{ | |
1360 | fHistSigToNoiseMC->SetBinContent(ipoint+1, -1) ; //-1 means: no background | |
1361 | } | |
1362 | } | |
1363 | //============================================================= | |
1364 | cout<<"--------------- Extracted Signal (MC) ------------------"<<endl; | |
1365 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1366 | cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tSignal: "<<lSigMCV0[ipoint]<<" +/- "<<lSigErrMCV0[ipoint]<<endl; | |
1367 | } | |
1368 | cout<<"--------------------------------------------------------"<<endl; | |
1369 | cout<<endl; | |
1370 | ||
1371 | ||
1372 | cout<<"--------------- Process Generated V0s ------------------"<<endl; | |
1373 | //Filling histogram with original MC particles==================== | |
1374 | TH3F* f3dHistGenPtVsYVsMultV0 = 0x0; | |
1375 | if(fWhichParticle == "Lambda") | |
1376 | f3dHistGenPtVsYVsMultV0 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda"); | |
1377 | if(fWhichParticle == "AntiLambda") | |
1378 | f3dHistGenPtVsYVsMultV0 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda"); | |
1379 | if(fWhichParticle == "K0Short") | |
1380 | f3dHistGenPtVsYVsMultV0 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short"); | |
d1d5d7e6 | 1381 | |
1382 | TH1D *fHistDummyV0 = 0x0; | |
1383 | if( fPerformMultiplicityStudy == kFALSE ){ | |
1384 | fHistDummyV0 = f3dHistGenPtVsYVsMultV0->ProjectionX("fHistDummyV0", | |
1385 | f3dHistGenPtVsYVsMultV0->GetYaxis()->FindBin(-fRapidityBoundary+1e-2), //avoid taking previous bin | |
1386 | f3dHistGenPtVsYVsMultV0->GetYaxis()->FindBin(+fRapidityBoundary-1e-2), //avoid taking next bin | |
1387 | f3dHistGenPtVsYVsMultV0->GetZaxis()->FindBin(-1), //Not interested in Multiplicity so far, integrate | |
1388 | f3dHistGenPtVsYVsMultV0->GetZaxis()->FindBin(101) //Not interested in Multiplicity so far, integrate | |
1389 | ); | |
1390 | } | |
1391 | if( fPerformMultiplicityStudy == kTRUE ){ | |
1392 | fHistDummyV0 = f3dHistGenPtVsYVsMultV0->ProjectionX("fHistDummyV0", | |
1393 | f3dHistGenPtVsYVsMultV0->GetYaxis()->FindBin(-fRapidityBoundary+1e-2), //avoid taking previous bin | |
1394 | f3dHistGenPtVsYVsMultV0->GetYaxis()->FindBin(+fRapidityBoundary-1e-2), //avoid taking next bin | |
1395 | f3dHistGenPtVsYVsMultV0->GetZaxis()->FindBin(fLoMultBound), //Interested in Multiplicity! | |
1396 | f3dHistGenPtVsYVsMultV0->GetZaxis()->FindBin(fHiMultBound) //Interested in Multiplicity! | |
1397 | ); | |
1398 | } | |
57a6d361 | 1399 | |
1400 | TH1F *fHistMCCountbyptV0 = new TH1F("fHistMCCountbyptV0","V0 MC count;p_{T} (GeV/c);Counts",fptbinnumb,fptbinlimits); | |
1401 | ||
1402 | Double_t temppt; | |
1403 | for(long i = 1; i<fHistDummyV0->GetNbinsX()+1;i++){ | |
1404 | temppt = fHistDummyV0->GetXaxis()->GetBinCenter(i); | |
1405 | for(long filling = 0; filling<fHistDummyV0->GetBinContent(i);filling++){ | |
1406 | fHistMCCountbyptV0->Fill(temppt); | |
1407 | } | |
1408 | } | |
1409 | cout<<"--------------- Generated V0 Dump ----------------------"<<endl; | |
1410 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1411 | cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tSignal: "<<fHistMCCountbyptV0->GetBinContent(ipoint+1)<<" +/- "<<TMath::Sqrt(fHistMCCountbyptV0->GetBinContent(ipoint+1))<<endl; | |
1412 | } | |
1413 | cout<<"--------------------------------------------------------"<<endl; | |
1414 | ||
d1d5d7e6 | 1415 | //--------------------------------------------------------------- |
1416 | //extremely curious test!!! experimental only, FIXME | |
1417 | cout<<"Resolution tests ongoing - Experimental! "<<endl; | |
1418 | cout<<"--------------------------------------------------------"<<endl; | |
1419 | cout<<" Usual tests: (Pt - PtMC) histograms "<<endl; | |
1420 | cout<<" Fitting gaussians..."<<endl; | |
1421 | for(Int_t ibin=0;ibin<fptbinnumb;ibin++){ | |
1422 | lHistResolutionGaussian[ibin]->SetParameter(0, lHistResolution[ibin]->GetMaximum() ); | |
1423 | lHistResolutionGaussian[ibin]->SetParameter(1, lHistResolution[ibin]->GetMean() ); | |
1424 | lHistResolutionGaussian[ibin]->SetParameter(2, lHistResolution[ibin]->GetRMS() ); | |
1425 | sprintf(lHistResolutionName,"lHistResolutionGaussian%i",ibin); | |
1426 | lHistResolution[ibin]->Fit(lHistResolutionName,"IREM0"); | |
1427 | } | |
1428 | ||
1429 | Double_t lTotSigMCPtMC[100]; | |
1430 | Double_t lTotSigMCPtMCXCheck[100]; | |
1431 | Double_t lTotSigMCPtReco[100]; | |
1432 | Double_t lTotSigMCPtRecoXCheck[100]; | |
1433 | for(Int_t ibin=0;ibin<100;ibin++){ | |
1434 | lTotSigMCPtMC[ibin] = 0; | |
1435 | lTotSigMCPtReco[ibin] = 0; | |
1436 | lTotSigMCPtRecoXCheck[ibin] = 0; | |
1437 | } | |
1438 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ //reco pt | |
1439 | for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt | |
1440 | lTotSigMCPtMC[jbin] += ((double)f2dHistPtResolution->GetBinContent(ibin+1,jbin+1)); | |
1441 | lTotSigMCPtReco[ibin] += ((double)f2dHistPtResolution->GetBinContent(ibin+1,jbin+1)); | |
1442 | } | |
1443 | } | |
1444 | ||
1445 | ||
1446 | TMatrixD fResolutionMatrix(fptbinnumb,fptbinnumb); | |
1447 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ //reco pt | |
1448 | for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt | |
1449 | fResolutionMatrix[ibin][jbin] = ((double)f2dHistPtResolution->GetBinContent(ibin+1,jbin+1)) / lTotSigMCPtMC[jbin] ; | |
1450 | f2dHistPtBlur->SetBinContent(ibin+1,jbin+1,(float)fResolutionMatrix[ibin][jbin]); | |
1451 | } | |
1452 | } | |
1453 | fResolutionMatrix.Print(); | |
1454 | cout<<"---> Debug Test 1: Check if fResolutionMatrix * lTotSigMCPtMC is lTotSigMCPtReco"<<endl; | |
1455 | cout<<"---> Multiplying... "<<endl; | |
1456 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ //reco pt | |
1457 | for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt | |
1458 | //if resolution is perfect: lTempSigReal[ibin] = lSigRealV0[ibin] | |
1459 | // (as fResolution is delta[ibin][jbin]) | |
1460 | lTotSigMCPtRecoXCheck[ibin] += lTotSigMCPtMC[jbin]*fResolutionMatrix[ibin][jbin]; | |
1461 | } | |
1462 | } | |
1463 | ||
1464 | cout<<"---> Compare for Cross-check!..."<<endl; | |
1465 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ | |
1466 | cout<<"---> bin #"<<ibin<<"\t"<<lTotSigMCPtRecoXCheck[ibin]<<"\t"<<lTotSigMCPtReco[ibin]<<endl; | |
1467 | } | |
1468 | ||
1469 | ||
1470 | cout<<"Resolution matrix has determinant: "<<fResolutionMatrix.Determinant()<<endl; | |
1471 | if(TMath::Abs(fResolutionMatrix.Determinant()) > 1e-5){ | |
1472 | cout<<"Determinant is different from zero enough to do inversion."<<endl; | |
1473 | cout<<"Will, however, suppress data too far from diagonal!"<<endl; | |
1474 | cout<<"Truncation: meant to ensure better numerical stability"<<endl; | |
1475 | ||
1476 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ //reco pt | |
1477 | for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt | |
1478 | if ( TMath::Abs(ibin-jbin) > 1 ){ | |
1479 | fResolutionMatrix[ibin][jbin] = 0; | |
1480 | } | |
1481 | } | |
1482 | } | |
1483 | ||
1484 | fResolutionMatrix.Invert(); | |
1485 | fResolutionMatrix.Print(); | |
1486 | } | |
1487 | ||
1488 | fResolutionMatrix.Print(); | |
1489 | cout<<"---> Debug Test 2: Check if fResolutionMatrix^-1 * lTotSigMCPtReco is lTotSigMCPtMC"<<endl; | |
1490 | cout<<"---> AKA \"the big test\"!"<<endl; | |
1491 | cout<<"---> Multiplying... "<<endl; | |
1492 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ //reco pt | |
1493 | for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt | |
1494 | //if resolution is perfect: lTempSigReal[ibin] = lSigRealV0[ibin] | |
1495 | // (as fResolution is delta[ibin][jbin]) | |
1496 | lTotSigMCPtMCXCheck[ibin] += lTotSigMCPtReco[jbin]*fResolutionMatrix[ibin][jbin]; | |
1497 | } | |
1498 | } | |
1499 | cout<<"---> Compare for Cross-check!..."<<endl; | |
1500 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ | |
1501 | cout<<"---> bin #"<<ibin<<"\t"<<lTotSigMCPtMCXCheck[ibin]<<"\t"<<lTotSigMCPtMC[ibin]<<endl; | |
1502 | } | |
1503 | ||
1504 | ||
1505 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ //reco pt | |
1506 | for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt | |
1507 | f2dHistPtSharpen->SetBinContent(ibin+1,jbin+1,(float)fResolutionMatrix[ibin][jbin]); | |
1508 | } | |
1509 | } | |
1510 | ||
1511 | Double_t lTempSigRealV0[100]; | |
1512 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ | |
1513 | lTempSigRealV0[ibin] = 0; | |
1514 | } | |
1515 | //To use: de-convolute by doing matrix multiplication with lSigRealV0[]... | |
1516 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ //reco pt | |
1517 | for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt | |
1518 | //if resolution is perfect: lTempSigReal[ibin] = lSigRealV0[ibin] | |
1519 | // (as fResolution is delta[ibin][jbin]) | |
1520 | lTempSigRealV0[ibin] += lSigRealV0[jbin]*fResolutionMatrix[ibin][jbin]; | |
1521 | } | |
1522 | } | |
1523 | ||
1524 | cout<<"Compare considering resolutions"<<endl; | |
1525 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ | |
1526 | cout<<"bin #"<<ibin<<"\t"<<lTempSigRealV0[ibin]<<"\t"<<lSigRealV0[ibin]<<endl; | |
1527 | } | |
1528 | cout<<endl; | |
1529 | cout<<"inspect ratio for trends..."<<endl; | |
1530 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ | |
1531 | cout<<"bin #"<<ibin<<"\t"<<lTempSigRealV0[ibin]/lSigRealV0[ibin]<<endl; | |
1532 | } | |
1533 | //--------------------------------------------------------------- | |
1534 | ||
1535 | //============================================================= | |
57a6d361 | 1536 | // Compute Efficiency |
1537 | Double_t lEfficiency[100]; | |
1538 | Double_t lEfficiencyError[100]; | |
1539 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1540 | lEfficiency[ipoint] = lSigMCV0[ipoint] / fHistMCCountbyptV0->GetBinContent(ipoint+1); | |
1541 | lEfficiencyError[ipoint] = ErrorInRatio( | |
1542 | lSigMCV0[ipoint], | |
1543 | lSigErrMCV0[ipoint], | |
1544 | fHistMCCountbyptV0->GetBinContent(ipoint+1), | |
1545 | TMath::Sqrt(fHistMCCountbyptV0->GetBinContent(ipoint+1)) | |
1546 | ); | |
1547 | } | |
1548 | //============================================================= | |
1549 | cout<<endl; | |
1550 | cout<<"--------------- Memory Cleanup -------------------------"<<endl; | |
1551 | fHistDummyV0->Delete(); | |
1552 | fHistMCCountbyptV0->Delete(); | |
1553 | v0listMC->Delete(); | |
1554 | delete v0listMC; | |
1555 | fileMC->Close("R"); | |
1556 | fileMC->Delete(); | |
1557 | delete fileMC; | |
1558 | cout<<endl; | |
1559 | cout<<"--------------- Pure Efficiency Numbers ----------------"<<endl; | |
1560 | TH1F* fHistPureEfficiency = new TH1F("fHistPureEfficiency","Pure Efficiency;p_{T} (GeV/c);Pure Efficiency",fptbinnumb,fptbinlimits); | |
1561 | TH1F* fHistEfficiency = new TH1F("fHistEfficiency","Efficiency;p_{T} (GeV/c);Efficiency",fptbinnumb,fptbinlimits); | |
1562 | if( fWhichParticle == "K0Short" ) fHistPureEfficiency->SetTitle("K^{0}_{S} Efficiency"); | |
1563 | if( fWhichParticle == "Lambda" ) fHistPureEfficiency->SetTitle("#Lambda Efficiency (no G3/F corr.)"); | |
1564 | if( fWhichParticle == "AntiLambda" ) fHistPureEfficiency->SetTitle("#bar{#Lambda} Efficiency (no G3/F corr.)"); | |
1565 | ||
1566 | if( fWhichParticle == "Lambda" ) fHistEfficiency->SetTitle("#Lambda Efficiency (with G3/F corr.)"); | |
1567 | if( fWhichParticle == "AntiLambda" ) fHistEfficiency->SetTitle("#bar{#Lambda} Efficiency (with G3/F corr.)"); | |
1568 | ||
1569 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1570 | cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tEff.: "<<lEfficiency[ipoint]<<" +/- "<<lEfficiencyError[ipoint]<<endl; | |
1571 | fHistPureEfficiency->SetBinContent(ipoint+1, lEfficiency[ipoint]); | |
1572 | fHistPureEfficiency->SetBinError(ipoint+1, lEfficiencyError[ipoint]); | |
1573 | } | |
1574 | cout<<"--------------------------------------------------------"<<endl; | |
1575 | ||
1576 | ||
1577 | cout<<endl; | |
1578 | //If it's Lambda or AntiLambda, do Geant3/fluka correction | |
1579 | //Keep fit function outside 'if' scope | |
1580 | TF1 *fitGeant3FlukaCorr = 0x0; | |
1581 | TH1D* h1dPanosCorrections =0x0; | |
1582 | if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
1583 | cout<<"--------------- Geant3/Fluka Correction ----------------"<<endl; | |
1584 | ||
1585 | //Credit for the original version of this code: Antonin Maire (thanks!) | |
1586 | ||
1587 | TFile *fileCorrGeanT3FlukaPanos = new TFile( "aCorrectionForPPbarCrossSections_FromPanos.root", "READ"); | |
1588 | TH2D* h2dCorrForCrossSectionProtons = (TH2D*) ( fileCorrGeanT3FlukaPanos->Get("gHistCorrectionForCrossSectionProtons") ); | |
1589 | TH2D* h2dCorrCrossSectionAntiProtons = (TH2D*) ( fileCorrGeanT3FlukaPanos->Get("gHistCorrectionForCrossSectionAntiProtons") ); | |
1590 | //Find bins... | |
1591 | ||
1592 | //TH2D* h2dAsMCPtBaryonVsPtCasc = 0x0; | |
1593 | TH2D* h2dCorrCrossSection = 0x0; | |
1594 | ||
1595 | if( fWhichParticle == "Lambda" ) h2dCorrCrossSection = h2dCorrForCrossSectionProtons; | |
1596 | if( fWhichParticle == "AntiLambda" ) h2dCorrCrossSection = h2dCorrCrossSectionAntiProtons; | |
1597 | ||
1598 | Int_t myFirstYBinPanos = h2dCorrCrossSection ->GetXaxis()->FindBin(-fRapidityBoundary); | |
1599 | Int_t myLastYBinPanos = h2dCorrCrossSection ->GetXaxis()->FindBin(+fRapidityBoundary-1e-3); | |
1600 | ||
1601 | cout<<"Check limiting bins: "<<endl; | |
1602 | Printf(" - Bin %d = %f to %f", myFirstYBinPanos, h2dCorrCrossSection ->GetXaxis()->GetBinLowEdge(myFirstYBinPanos), h2dCorrCrossSection ->GetXaxis()->GetBinUpEdge(myFirstYBinPanos) ); | |
1603 | Printf(" - Bin %d = %f to %f", myLastYBinPanos, h2dCorrCrossSection ->GetXaxis()->GetBinLowEdge(myLastYBinPanos ), h2dCorrCrossSection ->GetXaxis()->GetBinUpEdge(myLastYBinPanos) ); | |
1604 | ||
1605 | // B.1.a. | |
1606 | h1dPanosCorrections = h2dCorrCrossSection->ProjectionY("h1dPanosCorrections", myFirstYBinPanos, myLastYBinPanos, "e" ); // ~mid-rapidity | |
1607 | h1dPanosCorrections->SetDirectory(0); //means manual cleanup later | |
1608 | h1dPanosCorrections->Scale( (Double_t) (1.0/(myLastYBinPanos - myFirstYBinPanos +1.0)) ); | |
1609 | h1dPanosCorrections->SetYTitle("#epsilon_{GEANT3} / #epsilon_{FLUKA}"); | |
1610 | ||
1611 | //TF1 *fitGeant3FlukaCorr = 0x0; declared outside this scope | |
1612 | if(fWhichParticle == "Lambda") fitGeant3FlukaCorr = new TF1("fitGeant3FlukaCorr", this, &AliV0Module::MyGeant3FlukaCorrectionForProtons, 0.25, 1.1, 3, "AliV0Module","MyGeant3FlukaCorrectionForProtons"); | |
1613 | if(fWhichParticle == "AntiLambda") fitGeant3FlukaCorr = new TF1("fitGeant3FlukaCorr", this, &AliV0Module::MyGeant3FlukaCorrectionForAntiProtons, 0.25, 1.1, 4, "AliV0Module", "MyGeant3FlukaCorrectionForAntiProtons"); | |
1614 | ||
1615 | h1dPanosCorrections->Fit("fitGeant3FlukaCorr","rime+0"); // Chi2 | |
1616 | ||
1617 | Printf("Test fit function..."); | |
1618 | Printf(" - fitGeant3FlukaCorr for pt = .25 GeV/c : %f", fitGeant3FlukaCorr->Eval( .25) ); | |
1619 | Printf(" - fitGeant3FlukaCorr for pt = .5 GeV/c : %f", fitGeant3FlukaCorr->Eval( .5) ); | |
1620 | Printf(" - fitGeant3FlukaCorr for pt = 1 GeV/c : %f", fitGeant3FlukaCorr->Eval( 1.0) ); | |
1621 | Printf(" - fitGeant3FlukaCorr for pt = 2 GeV/c : %f", fitGeant3FlukaCorr->Eval( 2.0) ); | |
1622 | Printf(" - fitGeant3FlukaCorr for pt = 5 GeV/c : %f", fitGeant3FlukaCorr->Eval( 5.0) ); | |
1623 | Printf(" - fitGeant3FlukaCorr for pt = 7 GeV/c : %f", fitGeant3FlukaCorr->Eval( 7.0) ); | |
1624 | Printf(" - fitGeant3FlukaCorr for pt = 8 GeV/c : %f", fitGeant3FlukaCorr->Eval( 8.0) ); | |
1625 | Printf(" - fitGeant3FlukaCorr for pt = 10 GeV/c : %f", fitGeant3FlukaCorr->Eval(10.0) ); | |
1626 | cout<<"--------------- Embedding G3/F in Efficiencies ---------"<<endl; | |
1627 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1628 | lEfficiency[ipoint] /= fitGeant3FlukaCorr->Eval( lProtonMomentum[ipoint]->GetMean() ) ; | |
1629 | lEfficiencyError[ipoint] /= fitGeant3FlukaCorr->Eval( lProtonMomentum[ipoint]->GetMean() ) ; | |
1630 | } | |
1631 | cout<<"--------------- G3/F Corrected Efficiencies ------------"<<endl; | |
1632 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1633 | cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tEff.: "<<lEfficiency[ipoint]<<" +/- "<<lEfficiencyError[ipoint]<<endl; | |
1634 | fHistEfficiency->SetBinContent(ipoint+1, lEfficiency[ipoint]); | |
1635 | fHistEfficiency->SetBinError(ipoint+1, lEfficiencyError[ipoint]); | |
1636 | } | |
1637 | cout<<"--------------------------------------------------------"<<endl; | |
1638 | cout<<endl; | |
1639 | ||
1640 | fileCorrGeanT3FlukaPanos->Close(); | |
1641 | fileCorrGeanT3FlukaPanos->Delete(); | |
1642 | delete fileCorrGeanT3FlukaPanos; | |
1643 | ||
1644 | } //end Geant3/fluka correction if | |
1645 | ||
1646 | //Declaring outside Feeddown scope to keep in memory | |
1647 | TH2F *f2dFeedDownMatrix = 0x0; | |
1648 | TH2F *f2dFeedDownEfficiency = 0x0; | |
1649 | TH2F *f2dFeedDownEfficiencyGFCorrected = 0x0; | |
1650 | TH1F *fHistFeeddownSubtraction = 0x0; | |
1651 | ||
1652 | //========================================================================= | |
1653 | // --- Feeddown correction section | |
1654 | if( fFDSwitch != "NoFD" && fWhichParticle != "K0Short"){ | |
1655 | cout<<"--------------- Feeddown Correction --------------------"<<endl; | |
1656 | //MC File Acquisition============================================= | |
1657 | cout<<"--------------- Opening MC file for Feeddown -----------"<<endl; | |
1658 | TFile* fileMCFD = TFile::Open(fFeedDownDataFile, "READ"); | |
1659 | fileMCFD->cd("PWG2CheckPerformanceLambda_PP_MC"); | |
1660 | TList* v0listMCFD = (TList*)fileMCFD->FindObjectAny("clistV0MC"); | |
1661 | TTree* lTreeMCFD; | |
1662 | TH1F* fHistMultiplicityBeforeTrigSelMCFD; | |
1663 | lTreeMCFD = (TTree*)fileMCFD->FindObjectAny("fTree"); | |
1664 | fHistMultiplicityBeforeTrigSelMCFD = (TH1F*)v0listMCFD->FindObject("fHistMultiplicityBeforeTrigSel"); | |
1665 | ||
1666 | Long_t lNCandidatesMCFD = lTreeMCFD->GetEntries(); | |
1667 | Long_t lNEventsMCFD = fHistMultiplicityBeforeTrigSelMCFD->GetEntries(); | |
1668 | cout<<" MC Events (before trig sel)...: "<<lNEventsMCFD<<endl; | |
1669 | cout<<" Cascade MC Candidates.........: "<<lNCandidatesMCFD<<endl; | |
1670 | cout<<"--------------------------------------------------------"<<endl; | |
1671 | //================================================================ | |
1672 | ||
1673 | //Linking to Tree================================================= | |
1674 | //--- Base Variables ---------------------------------------------- | |
1675 | lTreeMCFD->SetBranchAddress("fTreeVariablePosEta",&lPosEta); | |
1676 | lTreeMCFD->SetBranchAddress("fTreeVariableNegEta",&lNegEta); | |
1677 | lTreeMCFD->SetBranchAddress("fTreeVariablePrimaryStatus",&lPrimaryStatus); | |
1678 | lTreeMCFD->SetBranchAddress("fTreeVariablePrimaryStatusMother",&lPrimaryStatusMother); | |
1679 | lTreeMCFD->SetBranchAddress("fTreeVariablePt",&lPt); | |
1680 | lTreeMCFD->SetBranchAddress("fTreeVariablePtMC",&lPtMC); | |
1681 | lTreeMCFD->SetBranchAddress("fTreeVariablePtXiMother",&lPtMother); | |
1682 | if ( fWhichParticle == "Lambda" ) lTreeMCFD->SetBranchAddress("fTreeVariableInvMassLambda",&lInvariantMass); | |
1683 | if ( fWhichParticle == "AntiLambda" ) lTreeMCFD->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMass); | |
1684 | if ( fWhichParticle == "K0Short" ) lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMass); | |
1685 | if ( fWhichParticle == "Lambda" ){ //For symmetry of computation... | |
1686 | lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne); | |
1687 | lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo); | |
1688 | } | |
1689 | if ( fWhichParticle == "AntiLambda" ){ //For symmetry of computation... | |
1690 | lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne); | |
1691 | lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo); | |
1692 | } | |
1693 | if ( fWhichParticle == "K0Short" ){ | |
1694 | lTreeMCFD->SetBranchAddress("fTreeVariableInvMassLambda" ,&lInvariantMassCompetingOne); | |
1695 | lTreeMCFD->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMassCompetingTwo); | |
1696 | } | |
1697 | //if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ) | |
1698 | // lTreeMCFD->SetBranchAddress("fTreeVariableRapLambda",&lRap); | |
1699 | //if ( fWhichParticle == "K0Short" ) | |
1700 | // lTreeMCFD->SetBranchAddress("fTreeVariableRapK0Short",&lRap); | |
1701 | ||
1702 | lTreeMCFD->SetBranchAddress("fTreeVariableRapMC",&lRap); | |
1703 | lTreeMCFD->SetBranchAddress("fTreeVariableLeastNbrCrossedRows",&lLeastNbrCrossedRows); | |
1704 | lTreeMCFD->SetBranchAddress("fTreeVariableLeastRatioCrossedRowsOverFindable",&lLeastNbrCrossedRowsOverFindable); | |
1705 | //--- Topological selection variables ----------------------------- | |
1706 | lTreeMCFD->SetBranchAddress("fTreeVariableV0Radius",&lV0Radius); | |
1707 | lTreeMCFD->SetBranchAddress("fTreeVariableDcaNegToPrimVertex",&lDcaNegToPrimVertex); | |
1708 | lTreeMCFD->SetBranchAddress("fTreeVariableDcaPosToPrimVertex",&lDcaPosToPrimVertex); | |
1709 | lTreeMCFD->SetBranchAddress("fTreeVariableDcaV0Daughters",&lDcaV0Daughters); | |
1710 | lTreeMCFD->SetBranchAddress("fTreeVariableV0CosineOfPointingAngle",&lV0CosinePointingAngle); | |
1711 | lTreeMCFD->SetBranchAddress("fTreeVariableDistOverTotMom",&lDistOverTotMom); | |
1712 | //---- Extra Information Available in MC only -------------------- | |
1713 | lTreeMCFD->SetBranchAddress("fTreeVariableNegTransvMomentumMC",&lNegTransvMomentumMC); | |
1714 | lTreeMCFD->SetBranchAddress("fTreeVariablePosTransvMomentumMC",&lPosTransvMomentumMC); | |
1715 | lTreeMCFD->SetBranchAddress("fTreeVariablePID",&lPID); | |
1716 | lTreeMCFD->SetBranchAddress("fTreeVariablePIDMother",&lPIDMother); | |
1717 | lTreeMCFD->SetBranchAddress("fTreeVariablePIDPositive",&lPIDPositive); | |
1718 | lTreeMCFD->SetBranchAddress("fTreeVariablePIDNegative",&lPIDNegative); | |
d1d5d7e6 | 1719 | //---- Multiplicity info ------------------------------------------- |
1720 | lTreeMCFD->SetBranchAddress("fTreeVariableMultiplicity",&lMultiplicity); | |
1721 | //--- Armenteros-Podolansky ---------------------------------------- | |
1722 | lTreeMCFD->SetBranchAddress("fTreeVariableAlphaV0",&lArmAlpha); | |
1723 | lTreeMCFD->SetBranchAddress("fTreeVariablePtArmV0",&lArmPt); | |
57a6d361 | 1724 | //================================================================ |
1725 | ||
1726 | //================================================================ | |
1727 | // Define Feeddown matrix | |
1728 | Double_t lFeedDownMatrix [100][100]; | |
1729 | // FeedDownMatrix [Lambda Bin][Xi Bin]; | |
1730 | for(Int_t ilb = 0; ilb<100; ilb++){ | |
1731 | for(Int_t ixb = 0; ixb<100; ixb++){ | |
1732 | lFeedDownMatrix[ilb][ixb]=0; | |
1733 | } | |
1734 | } | |
1735 | //Define Xi Binning | |
1736 | Double_t xibinlimits[] = { | |
1737 | 0.00, 0.20, 0.40, 0.60, 0.80, 0.90, | |
1738 | 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, | |
1739 | 1.70, 1.90, 2.20, 2.60, 3.10, 3.90, | |
1740 | 4.90, 6.00, 7.20, 8.50 ,10.00, 12.00 }; | |
1741 | Long_t xibinnumb = sizeof(xibinlimits)/sizeof(Double_t) - 1; | |
1742 | ||
1743 | TH1F* fHistPtXiReference = | |
1744 | new TH1F("fHistPtXiReference","#Xi candidates uncorrected p_{T};p_{T} (GeV/c);Counts",xibinnumb,xibinlimits); | |
1745 | ||
1746 | //Feeddown: will use a double-coordinate system: | |
1747 | // ( lambda bin, xi bin ) all the time! | |
1748 | lWeAreAtBin = 0; // Lambda bin | |
1749 | Int_t lWeAreAtXiBin = 0; // Xi Bin | |
1750 | ||
1751 | cout<<"--------------- MC Data File Loop, Feeddown -----------"<<endl; | |
1752 | Long_t lOneTenthOfNCandidatesMCFD = ((double)(lNCandidatesMCFD) / 10. ); | |
1753 | for(Long_t icand = 0;icand<lNCandidatesMCFD;icand++){ | |
1754 | lTreeMCFD->GetEntry(icand); | |
1755 | if( icand % lOneTenthOfNCandidatesMCFD == 0 ) | |
1756 | cout<<" Currently at candidate........: "<<icand<<" / "<<lNCandidatesMCFD<<" ( "<<(long)(((double)(icand)/(double)(lNCandidatesMCFD))*(100.+1e-3))<<"% )"<<endl; | |
1757 | if(TMath::Abs(lRap)<fRapidityBoundary && | |
1758 | TMath::Abs(lNegEta) <= fCutDaughterEta && | |
1759 | TMath::Abs(lPosEta) <= fCutDaughterEta && | |
1760 | lV0Radius >= fCutV0Radius && | |
1761 | lDcaNegToPrimVertex >= fCutDCANegToPV && | |
1762 | lDcaPosToPrimVertex >= fCutDCAPosToPV && | |
1763 | lDcaV0Daughters <= fCutDCAV0Daughters && | |
1764 | lV0CosinePointingAngle >= fCutV0CosPA && | |
1765 | lParticleMass*lDistOverTotMom <= fCutProperLifetime && | |
1766 | lLeastNbrCrossedRows >= fCutLeastNumberOfCrossedRows && | |
1767 | lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable && | |
1768 | TMath::Abs(lInvariantMassCompetingOne - lCompetingParticleMass) > fCutCompetingV0Rejection && | |
1769 | TMath::Abs(lInvariantMassCompetingTwo - lCompetingParticleMass) > fCutCompetingV0Rejection && | |
25231937 | 1770 | ( fSpecialArmenterosCutK0s == kFALSE || ( fSpecialArmenterosCutK0s == kTRUE && lArmPt*5>TMath::Abs(lArmAlpha) ) ) && |
57a6d361 | 1771 | ( //perfect PID association, IsPhysicalPrimary association |
1772 | ( fWhichParticle == "Lambda" | |
1773 | && lPID == 3122 //V0 is a Lambda | |
1774 | && lPIDPositive == 2212 //Pos Daughter is p | |
1775 | && lPIDNegative == -211 //Neg Daughter is pi- | |
1776 | && (lPIDMother == 3312 || (lPIDMother == 3322 && fFDSwitch=="UseMCRatio") ) | |
1777 | && lPrimaryStatusMother == 1 //Xi is actually a primary (should matter little) | |
1778 | ) || | |
1779 | ( fWhichParticle == "AntiLambda" | |
1780 | && lPID == -3122 //V0 is an AntiLambda | |
1781 | && lPIDPositive == 211 //Pos Daughter is pi+ | |
1782 | && lPIDNegative == -2212 //Neg Daughter is antiproton | |
1783 | && (lPIDMother ==-3312 || (lPIDMother ==-3322 && fFDSwitch=="UseMCRatio") ) | |
1784 | && lPrimaryStatusMother == 1 //Xi is actually a primary (should matter little) | |
1785 | ) | |
d1d5d7e6 | 1786 | ) && |
1787 | ( //Multiplicity Switch | |
1788 | fPerformMultiplicityStudy == kFALSE || //No mult selection, or... | |
1789 | (fPerformMultiplicityStudy == kTRUE && //inside mult bin | |
1790 | lMultiplicity>=fLoMultBound && | |
1791 | lMultiplicity<=fHiMultBound) | |
57a6d361 | 1792 | ) |
1793 | ){ // Start Entry Loop | |
1794 | lWeAreAtBin = fHistPt->FindBin(lPtMC)-1; | |
1795 | if(lWeAreAtBin == -1) lWeAreAtBin = 99; //UnderFlow, special treatment | |
1796 | lWeAreAtXiBin = fHistPtXiReference->FindBin(lPtMother)-1; | |
1797 | if(lWeAreAtXiBin == -1) lWeAreAtXiBin = 99; //UnderFlow, special treatment | |
1798 | //Populate Feeddown Matrix | |
1799 | //cout<<" Populate at coordinates "<<WeAreAtBin<<" vs "<<WeAreAtXiBin<<" ....."<<endl; | |
1800 | lFeedDownMatrix[lWeAreAtBin][lWeAreAtXiBin]++; | |
1801 | } // End Entry Loop | |
1802 | } | |
1803 | cout<<"--------------- Loop Completed -------------------------"<<endl; | |
1804 | cout<<endl; | |
1805 | ||
1806 | ||
1807 | //Filling histogram with original MC particles==================== | |
1808 | TH3F* f3dHistGenPtVsYVsMultV0FeedDown = 0x0; | |
1809 | ||
1810 | if (fWhichParticle == "Lambda" ) | |
1811 | f3dHistGenPtVsYVsMultV0FeedDown = (TH3F*)v0listMCFD->FindObject("f3dHistGenPtVsYVsMultXiMinus"); | |
1812 | if (fWhichParticle == "AntiLambda" ) | |
1813 | f3dHistGenPtVsYVsMultV0FeedDown = (TH3F*)v0listMCFD->FindObject("f3dHistGenPtVsYVsMultXiPlus"); | |
1814 | ||
1815 | TH1D *fHistDummyV0FeedDown = f3dHistGenPtVsYVsMultV0FeedDown->ProjectionX("fHistDummyV0FeedDown", | |
1816 | f3dHistGenPtVsYVsMultV0FeedDown->GetYaxis()->FindBin(-fRapidityBoundary+1e-2), //avoid taking previous bin | |
1817 | f3dHistGenPtVsYVsMultV0FeedDown->GetYaxis()->FindBin(+fRapidityBoundary-1e-2), //avoid taking next bin | |
1818 | f3dHistGenPtVsYVsMultV0FeedDown->GetZaxis()->FindBin(-1), | |
1819 | f3dHistGenPtVsYVsMultV0FeedDown->GetZaxis()->FindBin(101) | |
1820 | ); | |
1821 | ||
1822 | TH1F *fHistMCCountbyptXiFeedDown = new TH1F("fHistMCCountbyptXiFeedDown","#Xi MC count;p_{T} (GeV/c);Counts",xibinnumb,xibinlimits); | |
1823 | //Double_t temppt, y; --- declared before | |
1824 | for(long i = 1; i<fHistDummyV0FeedDown->GetNbinsX()+1;i++){ | |
1825 | temppt = fHistDummyV0FeedDown->GetXaxis()->GetBinCenter(i); | |
1826 | for(long filling = 0; filling<fHistDummyV0FeedDown->GetBinContent(i);filling++){ | |
1827 | fHistMCCountbyptXiFeedDown->Fill(temppt); | |
1828 | } | |
1829 | } | |
1830 | if(fWhichParticle == "Lambda") | |
1831 | cout<<"--------------- Generated Xi- Dump ---------------------"<<endl; | |
1832 | if(fWhichParticle == "AntiLambda") | |
1833 | cout<<"--------------- Generated Xi+ Dump ---------------------"<<endl; | |
1834 | for(Int_t ipoint = 0; ipoint<xibinnumb; ipoint++){ | |
1835 | cout<<"---> ["<<xibinlimits[ipoint]<<" - "<<xibinlimits[ipoint+1]<<" GeV/c]\tSignal: "<<fHistMCCountbyptXiFeedDown->GetBinContent(ipoint+1)<<" +/- "<<TMath::Sqrt(fHistMCCountbyptXiFeedDown->GetBinContent(ipoint+1))<<endl; | |
1836 | } | |
1837 | cout<<"--------------------------------------------------------"<<endl; | |
1838 | cout<<endl; | |
1839 | cout<<"Computing actual feeddown matrix"<<endl; | |
1840 | Double_t lFeedDownEfficiency[100][100]; | |
1841 | Double_t lFeedDownEfficiencyGFCorrected[100][100]; | |
1842 | Double_t lFeedDownEfficiencyError[100][100]; | |
1843 | Double_t lFeedDownEfficiencyGFCorrectedError[100][100]; | |
1844 | ||
1845 | for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){ | |
1846 | for(Int_t ixb = 0; ixb<xibinnumb; ixb++){ | |
1847 | if( fHistMCCountbyptXiFeedDown->GetBinContent(ixb+1)!= 0 ){ //avoid div by zero | |
1848 | lFeedDownEfficiency[ilb][ixb]=((double)(lFeedDownMatrix[ilb][ixb])) / ((double)( fHistMCCountbyptXiFeedDown->GetBinContent(ixb+1) )) ; | |
1849 | lFeedDownEfficiencyError[ilb][ixb]=ErrorInRatio( | |
1850 | ((double)(lFeedDownMatrix[ilb][ixb])), | |
1851 | TMath::Sqrt((double)(lFeedDownMatrix[ilb][ixb])), | |
1852 | ((double)( fHistMCCountbyptXiFeedDown->GetBinContent(ixb+1) )), | |
1853 | TMath::Sqrt((double)( fHistMCCountbyptXiFeedDown->GetBinContent(ixb+1) )) | |
1854 | ); | |
1855 | }else{ | |
1856 | lFeedDownEfficiency[ilb][ixb] = 0; | |
1857 | lFeedDownEfficiencyError[ilb][ixb] = 0; | |
1858 | } | |
1859 | } | |
1860 | } | |
1861 | ||
1862 | //Beware: the feeddown correction efficiency is computed with MC, | |
1863 | //Which has the geant3/fluka problem. Thus, we actually geant3/fluka | |
1864 | //correct the Feeddown efficiencies, too... | |
1865 | for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){ | |
1866 | for(Int_t ixb = 0;ixb<xibinnumb; ixb++){ | |
1867 | lFeedDownEfficiencyGFCorrected[ipoint][ixb] = lFeedDownEfficiency[ipoint][ixb]/fitGeant3FlukaCorr->Eval( lProtonMomentum[ipoint]->GetMean() ) ; | |
1868 | lFeedDownEfficiencyGFCorrectedError[ipoint][ixb] = lFeedDownEfficiencyError[ipoint][ixb]/fitGeant3FlukaCorr->Eval( lProtonMomentum[ipoint]->GetMean() ) ; | |
1869 | } | |
1870 | } | |
1871 | ||
1872 | //Create FD TH2Fs for storing | |
1873 | f2dFeedDownMatrix = new TH2F("f2dFeedDownMatrix","",fptbinnumb,fptbinlimits,xibinnumb,xibinlimits); | |
1874 | f2dFeedDownEfficiency = new TH2F("f2dFeedDownEfficiency","",fptbinnumb,fptbinlimits,xibinnumb,xibinlimits); | |
1875 | f2dFeedDownEfficiencyGFCorrected = new TH2F("f2dFeedDownEfficiencyGFCorrected","",fptbinnumb,fptbinlimits,xibinnumb,xibinlimits); | |
1876 | f2dFeedDownMatrix->SetDirectory(0); | |
1877 | f2dFeedDownEfficiency->SetDirectory(0); | |
1878 | f2dFeedDownEfficiencyGFCorrected->SetDirectory(0); | |
1879 | for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){ | |
1880 | for(Int_t ixb = 0; ixb<xibinnumb; ixb++){ | |
1881 | f2dFeedDownMatrix->SetBinContent(ilb+1,ixb+1,lFeedDownMatrix[ilb][ixb]); | |
1882 | f2dFeedDownEfficiency->SetBinContent(ilb+1,ixb+1,lFeedDownEfficiency[ilb][ixb]); | |
1883 | f2dFeedDownEfficiency->SetBinError(ilb+1,ixb+1,lFeedDownEfficiencyError[ilb][ixb]); | |
1884 | f2dFeedDownEfficiencyGFCorrected->SetBinContent(ilb+1,ixb+1,lFeedDownEfficiencyGFCorrected[ilb][ixb]); | |
1885 | f2dFeedDownEfficiencyGFCorrected->SetBinError(ilb+1,ixb+1,lFeedDownEfficiencyGFCorrectedError[ilb][ixb]); | |
1886 | } | |
1887 | } | |
1888 | ||
1889 | Double_t lProducedXi[100]; | |
1890 | ||
d1d5d7e6 | 1891 | TF1 *lLevyFitXiFeedDown = new TF1("LevyFitXiFeedDown",this ,&AliV0Module::MyLevyPtXi, 0.0,15,3, "AliV0Module","MyLevyPXi"); |
57a6d361 | 1892 | |
1893 | //Set Parameters Adequately, as in paper | |
1894 | //FIXME: These are the 7TeV Xi- parameters!!! | |
1895 | //FIXME: Use points vs fit function (should matter little if Xi- fit is good) | |
1896 | ||
1897 | if(fWhichParticle == "Lambda"){ | |
1898 | lLevyFitXiFeedDown->SetParameter(0, 7.98e-03); | |
1899 | lLevyFitXiFeedDown->SetParameter(1, 3.4429e-1); | |
1900 | lLevyFitXiFeedDown->SetParameter(2, 1.0787e+1); | |
1901 | } | |
1902 | if(fWhichParticle == "AntiLambda"){ | |
1903 | lLevyFitXiFeedDown->SetParameter(0, 7.79e-03); | |
1904 | lLevyFitXiFeedDown->SetParameter(1, 3.3934e-1); | |
1905 | lLevyFitXiFeedDown->SetParameter(2, 1.0428e+1); | |
1906 | } | |
1907 | //If you want me to double charged Xi feeddown, I'll do it here | |
1908 | if( fFDSwitch == "DoubleChargedXi" ){ | |
1909 | lLevyFitXiFeedDown->SetParameter(0, lLevyFitXiFeedDown->GetParameter(0)*2 ); | |
1910 | } | |
1911 | for(Int_t ixb = 0; ixb<xibinnumb; ixb++){ | |
1912 | lProducedXi[ixb] = lNInelasticEvents * lLevyFitXiFeedDown->Integral(xibinlimits[ixb],xibinlimits[ixb+1]); | |
1913 | } | |
1914 | if(fWhichParticle == "Lambda") | |
1915 | cout<<"--------------- Generated Xi- Dump (real-corrected) ----"<<endl; | |
1916 | if(fWhichParticle == "AntiLambda") | |
1917 | cout<<"--------------- Generated Xi+ Dump (real-corrected) ----"<<endl; | |
1918 | for(Int_t ixb = 0; ixb<xibinnumb; ixb++){ | |
1919 | cout<<"Xi bin "<<ixb<<"\t"<<lProducedXi[ixb]<<endl; | |
1920 | } | |
1921 | cout<<"--------------------------------------------------------"<<endl; | |
1922 | cout<<endl; | |
1923 | ||
1924 | //=== Computing actual feeddown numbers... === | |
1925 | Double_t lFeedDownToSubtract[100]; | |
1926 | Double_t lFeedDownToSubtractError[100]; | |
1927 | for(Int_t ilb = 0; ilb < fptbinnumb; ilb++){ | |
1928 | lFeedDownToSubtract[ilb] = 0; | |
1929 | lFeedDownToSubtractError[ilb] = 0; | |
1930 | for(Int_t ixb = 0; ixb < xibinnumb; ixb++){ | |
1931 | lFeedDownToSubtract[ilb] += lProducedXi[ixb]*lFeedDownEfficiencyGFCorrected[ilb][ixb]; | |
1932 | lFeedDownToSubtractError[ilb] += TMath::Power( (lProducedXi[ixb]*lFeedDownEfficiencyGFCorrectedError[ilb][ixb]), 2 ); | |
1933 | } | |
1934 | lFeedDownToSubtractError[ilb] = TMath::Sqrt(lFeedDownToSubtractError[ilb]); | |
1935 | } | |
1936 | fHistFeeddownSubtraction = new TH1F("fHistFeeddownSubtraction","FD Subtraction;p_{T} (GeV/c);Ratio removed",fptbinnumb,fptbinlimits); | |
1937 | fHistFeeddownSubtraction->SetDirectory(0); | |
1938 | for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){ | |
1939 | fHistFeeddownSubtraction->SetBinContent( ilb+1 , ((double)(lFeedDownToSubtract[ilb])) / ((double)(lSigRealV0[ilb]) ) ); | |
1940 | fHistFeeddownSubtraction->SetBinError ( ilb+1 , | |
1941 | ErrorInRatio( | |
1942 | ((double)(lFeedDownToSubtract[ilb])), | |
1943 | ((double)(lFeedDownToSubtractError[ilb])), | |
1944 | ((double)(lSigRealV0[ilb]) ), | |
1945 | ((double)(lSigErrRealV0[ilb]) ) ) | |
1946 | ); | |
1947 | } | |
1948 | ||
1949 | ||
1950 | ||
1951 | cout<<"--------------- FD Subtraction Fraction ----------------"<<endl; | |
1952 | for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){ | |
1953 | cout<<"---> ["<<fptbinlimits[ilb]<<" - "<<fptbinlimits[ilb+1]<<" GeV/c]\t"<<fHistFeeddownSubtraction->GetBinContent(ilb+1)<<"\t+/-\t"<<fHistFeeddownSubtraction->GetBinError(ilb+1)<<endl; | |
1954 | } | |
1955 | cout<<"--------------------------------------------------------"<<endl; | |
1956 | cout<<endl; | |
1957 | cout<<"--------------------------------------------------------"<<endl; | |
1958 | cout<<" Performing Actual Correction... "<<endl; | |
1959 | for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){ | |
1960 | lSigRealV0 [ilb] = lSigRealV0[ilb] - lFeedDownToSubtract[ilb]; | |
1961 | lSigErrRealV0[ilb] = TMath::Sqrt( TMath::Power(lSigErrRealV0[ilb],2) + TMath::Power(lFeedDownToSubtractError[ilb],2) ); | |
1962 | } | |
1963 | cout<<"--------------------------------------------------------"<<endl; | |
1964 | cout<<endl; | |
1965 | cout<<"--------------- Memory Cleanup -------------------------"<<endl; | |
1966 | lLevyFitXiFeedDown->Delete(); | |
1967 | v0listMCFD->Delete(); | |
1968 | delete v0listMCFD; | |
1969 | fileMCFD->Close("R"); | |
1970 | fileMCFD->Delete(); | |
1971 | delete fileMCFD; | |
1972 | cout<<endl; | |
1973 | }//--- End Feeddown Correction section | |
1974 | //========================================================================= | |
1975 | ||
1976 | //========================================================================= | |
1977 | //---> At this stage, everthing's just ready for the actual spectrum | |
1978 | //---> computation to occur! | |
1979 | Double_t lSpectrum[100]; | |
1980 | Double_t lSpectrumError[100]; | |
1981 | ||
1982 | for(Int_t ibin=0;ibin<fptbinnumb;ibin++){ | |
1983 | lSpectrum[ibin] = lSigRealV0[ibin] / lEfficiency [ibin]; | |
1984 | lSpectrumError[ibin] = ErrorInRatio( | |
1985 | lSigRealV0 [ibin], | |
1986 | lSigErrRealV0 [ibin], | |
1987 | lEfficiency [ibin], | |
1988 | lEfficiencyError [ibin]); | |
1989 | } | |
1990 | ||
1991 | //Divide by: Bin Width, Rapidity Window, N_{Inel} | |
1992 | for(Int_t ibin=0;ibin<fptbinnumb;ibin++){ | |
1993 | lSpectrum[ibin] /= (fptbinlimits[ibin+1]-fptbinlimits[ibin]); | |
1994 | lSpectrum[ibin] /= lNInelasticEvents; | |
1995 | lSpectrum[ibin] /= 2*fRapidityBoundary; | |
1996 | lSpectrumError[ibin] /= (fptbinlimits[ibin+1]-fptbinlimits[ibin]); | |
1997 | lSpectrumError[ibin] /= lNInelasticEvents; | |
1998 | lSpectrumError[ibin] /= 2*fRapidityBoundary; | |
1999 | } | |
2000 | ||
2001 | TH1F* fHistPtLambda = new TH1F("fHistPtLambda","#Lambda Corrected Spectrum;p_{T} (GeV/c);1/N_{INEL} #frac{d^{2}N}{dydp_{t}}",fptbinnumb,fptbinlimits); | |
2002 | TH1F* fHistPtAntiLambda = new TH1F("fHistPtAntiLambda","#bar{#Lambda} Corrected Spectrum;p_{T} (GeV/c);1/N_{INEL} #frac{d^{2}N}{dydp_{t}}",fptbinnumb,fptbinlimits); | |
2003 | TH1F* fHistPtK0Short = new TH1F("fHistPtK0Short","K^{0}_{S} Corrected Spectrum;p_{T} (GeV/c);1/N_{INEL} #frac{d^{2}N}{dydp_{t}}",fptbinnumb,fptbinlimits); | |
2004 | ||
2005 | //Copy to Histogram | |
2006 | for(Int_t ibin=0;ibin<fptbinnumb;ibin++){ | |
2007 | if(fWhichParticle == "Lambda" ){ | |
2008 | fHistPtLambda->SetBinContent( ibin+1, lSpectrum[ibin] ); | |
2009 | fHistPtLambda->SetBinError ( ibin+1, lSpectrumError[ibin] ); | |
2010 | } | |
2011 | if(fWhichParticle == "AntiLambda" ){ | |
2012 | fHistPtAntiLambda->SetBinContent( ibin+1, lSpectrum[ibin] ); | |
2013 | fHistPtAntiLambda->SetBinError ( ibin+1, lSpectrumError[ibin] ); | |
2014 | } | |
2015 | if(fWhichParticle == "K0Short" ){ | |
2016 | fHistPtK0Short->SetBinContent( ibin+1, lSpectrum[ibin] ); | |
2017 | fHistPtK0Short->SetBinError ( ibin+1, lSpectrumError[ibin] ); | |
2018 | } | |
2019 | } | |
2020 | ||
2021 | ||
2022 | //========================================================================= | |
2023 | ||
2024 | cout<<"--------------- Result Output --------------------------"<<endl; | |
2025 | cout<<" ---> Writing information to "<<fOutputDataFile<<endl; | |
2026 | // Open an output file | |
2027 | TFile* lResultsFile = TFile::Open(fOutputDataFile, "RECREATE"); | |
2028 | if (!lResultsFile || !lResultsFile->IsOpen()){ | |
2029 | cout<<"Error! Couldn't open file!"<<endl; | |
2030 | return; | |
2031 | } | |
2032 | ||
2033 | //Preparing Signal Extraction Range Canvas | |
2034 | TCanvas *cSigExtRange = new TCanvas("cSigExtRange","Extraction Range",900,600); | |
2035 | cSigExtRange->SetFillColor(kWhite); | |
2036 | cSigExtRange->SetLeftMargin(0.17); | |
2037 | cSigExtRange->SetRightMargin(0.17); | |
2038 | cSigExtRange->cd(); | |
2039 | fHistSignalExtractionRange->SetFillColor(18); | |
2040 | fHistSignalExtractionRange->SetMarkerStyle(20); | |
2041 | if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
2042 | fHistSignalExtractionRange->GetYaxis()->SetRangeUser(1.115-0.08,1.115+0.08); | |
2043 | } | |
2044 | if( fWhichParticle == "K0Short" ){ | |
2045 | fHistSignalExtractionRange->GetYaxis()->SetRangeUser(0.498-0.15,0.498+0.15); | |
2046 | } | |
2047 | fHistSignalExtractionRange->Draw("E2"); | |
2048 | if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
2049 | fLDataUpper->Draw("same"); | |
2050 | fLDataLower->Draw("same"); | |
2051 | } | |
2052 | if( fWhichParticle == "K0Short" ){ | |
2053 | fKDataUpper->Draw("same"); | |
2054 | fKDataLower->Draw("same"); | |
2055 | } | |
2056 | //Saving Invariant Mass Plots (real data) | |
2057 | lResultsFile->cd(); | |
2058 | TDirectoryFile *lInvMassReal = new TDirectoryFile("lInvMassReal","Invariant Mass Plots (Real Data)"); | |
2059 | lInvMassReal->cd(); | |
2060 | cSigExtRange->Write(); | |
2061 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) { | |
2062 | lCanvasHistoFullV0[ibin] -> Write(); | |
2063 | } | |
2064 | ||
2065 | TDirectoryFile *lInvMassRealRawData = new TDirectoryFile("lInvMassRealRawData","Objects for Inv Mass Plots (Real Data)"); | |
2066 | lInvMassRealRawData->cd(); | |
2067 | fHistSignalExtractionRange->Write(); | |
2068 | fHistPeakPosition->Write(); | |
2069 | fHistSigToNoise->Write(); | |
2070 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) { | |
2071 | lHistoFullV0[ibin] -> Write(); | |
2072 | fgausPt[ibin] -> Write(); | |
2073 | if( fFitBackgroundSwitch ) lfitNoise[ibin] -> Write(); | |
2074 | } | |
2075 | if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
2076 | fLDataUpper->Write(); | |
2077 | fLDataLower->Write(); | |
2078 | } | |
2079 | if( fWhichParticle == "K0Short" ){ | |
2080 | fKDataUpper->Write(); | |
2081 | fKDataLower->Write(); | |
2082 | } | |
2083 | ||
2084 | //Saving Invariant Mass Plots (MC) | |
2085 | lResultsFile->cd(); | |
2086 | TDirectoryFile *lInvMassMC = new TDirectoryFile("lInvMassMC","Invariant Mass Plots (Monte Carlo)"); | |
2087 | lInvMassMC->cd(); | |
2088 | ||
2089 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) { | |
2090 | lCanvasHistoFullV0MC[ibin] -> Write(); | |
2091 | } | |
2092 | ||
2093 | TDirectoryFile *lInvMassMCRawData = new TDirectoryFile("lInvMassMCRawData","Objects for Inv Mass Plots (MC)"); | |
2094 | lInvMassMCRawData->cd(); | |
2095 | fHistPeakPositionMC->Write(); | |
2096 | fHistSigToNoiseMC->Write(); | |
2097 | f2dHistPtResolution->Write(); | |
d1d5d7e6 | 2098 | f2dHistPtBlur->Write(); |
2099 | f2dHistPtSharpen->Write(); | |
57a6d361 | 2100 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) { |
2101 | lHistoFullV0MC[ibin] ->Write(); | |
2102 | fgausPtMC[ibin] ->Write(); | |
2103 | if( fFitBackgroundSwitch ) lfitNoiseMC[ibin] -> Write(); | |
2104 | } | |
2105 | ||
2106 | //Saving Geant3/Fluka Correction Data (MC) | |
2107 | if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
2108 | lResultsFile->cd(); | |
2109 | TDirectoryFile *lGFCorrection = new TDirectoryFile("lGFCorrection","Geant3/Fluka Correction Histograms"); | |
2110 | lGFCorrection->cd(); | |
2111 | h1dPanosCorrections->Write(); | |
2112 | fitGeant3FlukaCorr->Write(); | |
2113 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) { | |
2114 | lProtonMomentum[ibin] ->Write(); | |
2115 | } | |
2116 | } | |
2117 | ||
2118 | //Saving Feeddown Correction information, if needed | |
2119 | if( (fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda")&&fFDSwitch!="NoFD" ){ | |
2120 | lResultsFile->cd(); | |
2121 | TDirectoryFile *lFeeddown = new TDirectoryFile("lFeeddown","Feeddown Subtraction Information"); | |
2122 | lFeeddown->cd(); | |
2123 | f2dFeedDownMatrix->Write(); | |
2124 | f2dFeedDownEfficiency->Write(); | |
2125 | f2dFeedDownEfficiencyGFCorrected->Write(); | |
2126 | fHistFeeddownSubtraction->Write(); | |
2127 | } | |
2128 | ||
d1d5d7e6 | 2129 | //Saving Resolution Information |
2130 | ||
2131 | //preparing... | |
2132 | for(Int_t ibin=0; ibin<fptbinnumb; ibin++){ | |
2133 | fHistResolutionVsPt->SetBinContent(ibin+1,lHistResolution[ibin]->GetMean()); | |
2134 | fHistResolutionVsPt->SetBinError(ibin+1,lHistResolution[ibin]->GetRMS()); | |
2135 | fHistResolutionVsPtDivByBinWidth->SetBinContent(ibin+1,lHistResolution[ibin]->GetMean() / (fptbinlimits[ibin+1]-fptbinlimits[ibin]) ); | |
2136 | fHistResolutionVsPtDivByBinWidth->SetBinError(ibin+1,lHistResolution[ibin]->GetRMS() / (fptbinlimits[ibin+1]-fptbinlimits[ibin]) ); | |
2137 | ||
2138 | fHistResolutionVsPtWithGaussians->SetBinContent(ibin+1,lHistResolutionGaussian[ibin]->GetParameter(1)); | |
2139 | fHistResolutionVsPtWithGaussians->SetBinError(ibin+1,lHistResolutionGaussian[ibin]->GetParameter(2)); | |
2140 | fHistResolutionVsPtDivByBinWidthWithGaussians->SetBinContent(ibin+1,lHistResolutionGaussian[ibin]->GetParameter(1) / (fptbinlimits[ibin+1]-fptbinlimits[ibin]) ); | |
2141 | fHistResolutionVsPtDivByBinWidthWithGaussians->SetBinError(ibin+1,lHistResolutionGaussian[ibin]->GetParameter(2) / (fptbinlimits[ibin+1]-fptbinlimits[ibin]) ); | |
2142 | } | |
2143 | ||
2144 | ||
2145 | lResultsFile->cd(); | |
2146 | TDirectoryFile *lDirResolution = new TDirectoryFile("lInfoResolution","Resolution Information"); | |
2147 | lDirResolution->cd(); | |
2148 | fHistResolutionVsPt->Write(); | |
2149 | fHistResolutionVsPtDivByBinWidth->Write(); | |
2150 | fHistResolutionVsPtWithGaussians->Write(); | |
2151 | fHistResolutionVsPtDivByBinWidthWithGaussians->Write(); | |
2152 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) { | |
2153 | lHistResolution[ibin] ->Write(); | |
2154 | lHistResolutionGaussian[ibin]->Write(); | |
2155 | } | |
2156 | ||
2157 | ||
57a6d361 | 2158 | lResultsFile->cd(); |
2159 | ||
2160 | if( fWhichParticle == "K0Short") fHistPureEfficiency->Write(); | |
2161 | if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
2162 | fHistPureEfficiency->Write(); | |
2163 | fHistEfficiency->Write(); | |
2164 | } | |
2165 | if(fWhichParticle == "Lambda" ){ | |
2166 | fHistPtLambda -> Write(); | |
2167 | } | |
2168 | if(fWhichParticle == "AntiLambda" ){ | |
2169 | fHistPtAntiLambda -> Write(); | |
2170 | } | |
2171 | if(fWhichParticle == "K0Short" ){ | |
2172 | fHistPtK0Short -> Write(); | |
2173 | } | |
2174 | ||
2175 | lResultsFile->Write(); | |
2176 | lResultsFile->Close(); | |
2177 | delete lResultsFile; | |
2178 | ||
2179 | //================================================ | |
2180 | //Manual Cleanup of all created Histograms | |
2181 | //================================================ | |
2182 | // needed if you want to re-run the whole thing without | |
2183 | // memory leaks (systematics, etc) | |
2184 | ||
2185 | //switch on if you want large amounts of printout | |
2186 | Bool_t lDebugCleaningProcess = kFALSE; | |
2187 | ||
2188 | if (lDebugCleaningProcess) cout<<"fHistPt->Delete()"<<endl; | |
2189 | fHistPt->Delete(); | |
2190 | if (lDebugCleaningProcess) cout<<"fHistPeakPosition->Delete()"<<endl; | |
2191 | fHistPeakPosition->Delete(); | |
2192 | if (lDebugCleaningProcess) cout<<"fHistPeakPositionMC->Delete()"<<endl; | |
2193 | fHistPeakPositionMC->Delete(); | |
2194 | if (lDebugCleaningProcess) cout<<"fHistSigToNoise->Delete()"<<endl; | |
2195 | fHistSigToNoise->Delete(); | |
2196 | if (lDebugCleaningProcess) cout<<"fHistSigToNoiseMC->Delete()"<<endl; | |
2197 | fHistSigToNoiseMC->Delete(); | |
2198 | if (lDebugCleaningProcess) cout<<"fHistSignalExtractionRange->Delete()"<<endl; | |
2199 | fHistSignalExtractionRange->Delete(); | |
2200 | if (lDebugCleaningProcess) cout<<"fKData*->Delete()"<<endl; | |
2201 | fKDataUpper->Delete(); | |
2202 | fKDataLower->Delete(); | |
2203 | fLDataUpper->Delete(); | |
2204 | fLDataLower->Delete(); | |
2205 | if(lDebugCleaningProcess) cout<<"f2dHistPtResolution->Delete()"<<endl; | |
2206 | f2dHistPtResolution->Delete(); | |
2207 | ||
2208 | if(lDebugCleaningProcess) cout<<"lfitNoise*[*]->Delete(); lSampleNoise*[*]->Delete()"<<endl; | |
2209 | if( fFitBackgroundSwitch ){ | |
2210 | for(long i=0; i<fptbinnumb; i++){ | |
2211 | lfitNoise[i] -> Delete(); | |
2212 | lSampleNoise[i] -> Delete(); | |
2213 | lfitNoiseMC[i] -> Delete(); | |
2214 | lSampleNoiseMC[i] -> Delete(); | |
2215 | } | |
2216 | } | |
2217 | ||
2218 | ||
2219 | //pt-by-pt histos | |
d1d5d7e6 | 2220 | fHistResolutionVsPt->Delete(); |
2221 | fHistResolutionVsPtDivByBinWidth->Delete(); | |
2222 | fHistResolutionVsPtWithGaussians->Delete(); | |
2223 | fHistResolutionVsPtDivByBinWidthWithGaussians->Delete(); | |
2224 | f2dHistPtBlur->Delete(); | |
2225 | f2dHistPtSharpen->Delete(); | |
57a6d361 | 2226 | if (lDebugCleaningProcess) cout<<"lHistoFullV0*[*]->Delete()"<<endl; |
2227 | for(Int_t ihist=0;ihist<100;ihist++){ | |
2228 | lHistoFullV0[ihist]->Delete(); | |
2229 | lHistoFullV0MC[ihist]->Delete(); | |
d1d5d7e6 | 2230 | lHistResolution[ihist]->Delete(); |
2231 | lHistResolutionGaussian[ihist]->Delete(); | |
57a6d361 | 2232 | } |
2233 | ||
2234 | if (lDebugCleaningProcess) cout<<"lLine*[*]->Delete()"<<endl; | |
2235 | //gaussian fit functions, drawing lines | |
2236 | for(Int_t ibin = 0; ibin<fptbinnumb; ibin++){ | |
2237 | lLineLeftMost[ibin] ->Delete(); | |
2238 | lLineLeft[ibin] ->Delete(); | |
2239 | lLineRight[ibin] ->Delete(); | |
2240 | lLineRightMost[ibin]->Delete(); | |
2241 | lLineLeftMostMC[ibin] ->Delete(); | |
2242 | lLineLeftMC[ibin] ->Delete(); | |
2243 | lLineRightMC[ibin] ->Delete(); | |
2244 | lLineRightMostMC[ibin]->Delete(); | |
2245 | ||
2246 | if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
2247 | fgausPt[ibin]->Delete(); | |
2248 | fgausPtMC[ibin]->Delete(); | |
2249 | } | |
2250 | if ( fWhichParticle == "K0Short"){ | |
2251 | fgausPt[ibin]->Delete(); | |
2252 | fgausPtMC[ibin]->Delete(); | |
2253 | } | |
2254 | } | |
2255 | if (lDebugCleaningProcess) cout<<"lCanvasHistoFullV0*[*]->Delete()"<<endl; | |
2256 | for(Int_t ihist=0;ihist<100;ihist++){ | |
2257 | lCanvasHistoFullV0[ihist]->Close(); | |
2258 | lCanvasHistoFullV0MC[ihist]->Close(); | |
2259 | delete lCanvasHistoFullV0[ihist]; | |
2260 | delete lCanvasHistoFullV0MC[ihist]; | |
2261 | } | |
2262 | ||
2263 | if (lDebugCleaningProcess) cout<<"cSigExtRange->Delete()"<<endl; | |
2264 | cSigExtRange->Close(); | |
2265 | delete cSigExtRange; | |
2266 | ||
2267 | if (lDebugCleaningProcess) cout<<"lProtonMomentum[*]->Delete()"<<endl; | |
2268 | //histograms for G3/F correction | |
2269 | for(Int_t ibin=0;ibin<100;ibin++){ | |
2270 | lProtonMomentum[ibin]->Delete(); | |
2271 | } | |
2272 | ||
2273 | if (lDebugCleaningProcess) cout<<"fHist*Efficiency->Delete()"<<endl; | |
2274 | //MC info, efficiencies | |
2275 | fHistPureEfficiency->Delete(); | |
2276 | fHistEfficiency->Delete(); | |
2277 | //data for G3/F, continued | |
2278 | if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){ | |
2279 | h1dPanosCorrections->Delete(); | |
2280 | fitGeant3FlukaCorr->Delete(); | |
2281 | } | |
2282 | ||
2283 | //histograms, feeddown | |
2284 | if (lDebugCleaningProcess) cout<<"f2dFeedDown*->Delete()"<<endl; | |
2285 | if( fFDSwitch != "NoFD" && fWhichParticle != "K0Short"){ | |
2286 | //fHistPtXiReference->Delete(); | |
2287 | //fHistMCCountbyptXiFeedDown->Delete(); | |
2288 | f2dFeedDownMatrix->Delete(); | |
2289 | f2dFeedDownEfficiency->Delete(); | |
2290 | f2dFeedDownEfficiencyGFCorrected->Delete(); | |
2291 | fHistFeeddownSubtraction->Delete(); | |
2292 | } | |
2293 | if (lDebugCleaningProcess) cout<<"fHistPt*->Delete()"<<endl; | |
2294 | //Corrected Spectra Histograms | |
2295 | fHistPtLambda->Delete(); | |
2296 | fHistPtAntiLambda->Delete(); | |
2297 | fHistPtK0Short->Delete(); | |
2298 | ||
2299 | //Exit Batch Mode | |
2300 | gROOT->SetBatch (kFALSE); | |
2301 | ||
2302 | cout<<"--------------------------------------------------------"<<endl; | |
2303 | cout<<" There, done! "<<endl; | |
2304 | cout<<endl; | |
2305 | ||
2306 | ||
2307 | ||
2308 | ||
2309 | } |