]>
Commit | Line | Data |
---|---|---|
1af4558e | 1 | |
88d02ae4 | 2 | class AliSpectraBothHistoManager; |
3 | class AliSpectraBothEventCuts; | |
4 | class AliSpectraBothTrackCuts; | |
5 | TString Charge[]={"Pos","Neg"}; | |
6 | TString Sign[]={"Plus","Minus"}; | |
7 | TString Particle[]={"Pion","Kaon","Proton"}; | |
8 | AliSpectraBothHistoManager* managerdata=0x0; | |
9 | AliSpectraBothEventCuts* ecutsdata=0x0; | |
10 | AliSpectraBothTrackCuts* tcutsdata=0x0; | |
11 | ||
12 | AliSpectraBothHistoManager* managermc=0x0; | |
13 | AliSpectraBothEventCuts* ecutsmc=0x0; | |
14 | AliSpectraBothTrackCuts* tcutsmc=0x0; | |
15 | ||
16 | Float_t TOFMatchingScalling[2]={-1,-1}; | |
17 | Int_t Color[3]={1,2,4}; | |
18 | Int_t Marker[6]={20,21,22,24,25,26}; | |
19 | Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p | |
20 | ||
a4ca64e3 | 21 | |
22 | Double_t FitRange[2]={-2.0,2.0}; | |
23 | Double_t CutRange[2]={-3.0,3.0}; | |
24 | Double_t minptformaterial[6]={0.0,0.2,0.0,0.0,0.2,0.0}; | |
25 | Double_t maxptformaterial[6]={0.0,0.6,1.3,0.0,0.6,0.0}; | |
26 | Double_t minptforWD[6]={0.2,100.0,0.3,0.2,100.0,0.3}; | |
27 | Double_t maxptforWD[6]={1.5,-100.0,2.0,1.5,-100.0,2.0}; | |
1af4558e | 28 | Double_t minRanges[3]={0.3,0.3,0.45}; |
29 | Double_t maxRanges[3]={1.5,1.2,2.2}; | |
cd28560b | 30 | Double_t TOFPIDsignalmatching[]={-1.0,-1.0,-1.0}; |
1af4558e | 31 | Double_t fMaxContaminationPIDMC=0.2; |
a4ca64e3 | 32 | |
33 | ||
88d02ae4 | 34 | enum ECharge_t { |
35 | kPositive, | |
36 | kNegative, | |
37 | kNCharges | |
38 | }; | |
b3ea73e1 | 39 | enum { |
a4ca64e3 | 40 | kdodca=0x1, //dca fits are made |
41 | kgeantflukaKaon=0x2,// geant fluka correction is used for kaons | |
42 | kgeantflukaProton=0x4, // geant fluka correction is used for protons and antiprotons | |
43 | knormalizationtoeventspassingPhySel=0x8,// spectra are divided by number of events passing physic selection | |
44 | kveretxcorrectionandbadchunkscorr=0x10, // correction for difference in z vertex distribution in data and MC and correction for bad chunks is applied | |
45 | kmcisusedasdata=0x20, // the result of the looping over MC is used as data input | |
46 | kdonotusedcacuts=0x40, // allows to use the constant dca cut for all pt bins not the pt dependet defined in stardrad track cuts 2011 | |
47 | kuseprimaryPIDcont=0x80, //pid contamination is calculated using only primiary particle in this case K should use dca fits | |
48 | knormalizationwithbin0integralsdata=0x100, // the normalization factor is calcualte using integral over z vertex distributions (in this case reconstructed vertex disitrbution uses z vertex for data) | |
49 | knormalizationwithbin0integralsMC=0x200, //in this case reconstructed vertex disitrbution uses z vertex for data, those to options will be use only if knormalizationtoeventspassingPhySel is not set | |
b6a74f8a | 50 | kuserangeonfigfile=0x400, // use of config file for dca fit settings |
1af4558e | 51 | kskipconcutonspectra=0x800, //do not use conPID<02 cut useful for syst. studies |
cd28560b | 52 | kuseTOFmatchingcorrection=0x1000, // if set tof matching correction is applied. |
53 | kuseTOFcorrforPIDsignalmatching=0x2000 // rescale the for spectra by the factor given in config files | |
1af4558e | 54 | |
b3ea73e1 | 55 | }; |
88d02ae4 | 56 | |
0c1898ba | 57 | Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false); |
a4ca64e3 | 58 | void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" ) |
88d02ae4 | 59 | { |
b6a74f8a | 60 | gStyle->SetOptStat(0); |
88d02ae4 | 61 | TH1::AddDirectory(kFALSE); |
62 | gSystem->Load("libCore.so"); | |
63 | gSystem->Load("libPhysics.so"); | |
64 | gSystem->Load("libTree"); | |
65 | gSystem->Load("libMatrix"); | |
66 | gSystem->Load("libSTEERBase"); | |
67 | gSystem->Load("libESD"); | |
68 | gSystem->Load("libAOD"); | |
69 | gSystem->Load("libANALYSIS"); | |
70 | gSystem->Load("libOADB"); | |
71 | gSystem->Load("libANALYSISalice"); | |
72 | gSystem->Load("libTENDER"); | |
73 | gSystem->Load("libCORRFW"); | |
74 | gSystem->Load("libPWGTools"); | |
75 | gSystem->Load("libPWGLFspectra"); | |
76 | ||
77 | gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/QAPlotsBoth.C"); | |
78 | Double_t mass[3]; | |
79 | mass[0] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); | |
80 | mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass(); | |
81 | mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass(); | |
82 | ||
8cec8e63 | 83 | TFormula* dcacutxy=0x0; |
84 | if(!(options&kdonotusedcacuts)) | |
85 | { | |
86 | ||
87 | AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE); | |
a4ca64e3 | 88 | TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep()); |
8cec8e63 | 89 | formulastring.ReplaceAll("pt","x"); |
90 | dcacutxy=new TFormula("dcacutxy",formulastring.Data()); | |
91 | } | |
a4ca64e3 | 92 | if(options&kuserangeonfigfile) |
93 | if(!ReadConfigFile(configfile)) | |
94 | return; | |
88d02ae4 | 95 | TList* lout=new TList(); |
a4ca64e3 | 96 | |
88d02ae4 | 97 | |
98 | TString indirname=Form("/output/train%s",outdate.Data()); | |
99 | //TString indirname("/output/train24072012"); | |
100 | if(outnamemc.Length()==0) | |
101 | outnamemc=outnamedata; | |
102 | cout<<indirname.Data()<<" "<<outnamemc.Data()<<endl; | |
103 | // Do the job | |
104 | ||
105 | ||
106 | OpenFile(indirname,outnamemc,true); | |
0c1898ba | 107 | OpenFile(indirname,outnamedata,false,((Bool_t)(options&kmcisusedasdata))); |
88d02ae4 | 108 | if(!managermc||!managerdata) |
109 | { | |
110 | cout<<managermc<<" "<<managerdata<<endl; | |
111 | return; | |
112 | } | |
113 | TH1F* rawspectradata[6]; | |
114 | TH1F* rawspectramc[6]; | |
115 | TH1F* MCTruth[6]; | |
116 | TH1F* eff[6]; | |
117 | TH1F* contallMC[6]; | |
118 | TH1F* contPID[6]; | |
119 | TH1F* contWD[6]; | |
120 | TH1F* contMat[6]; | |
121 | TH1F* confinal[6]; | |
2245e1bd | 122 | TH1F* contPIDpri[6]; |
b6a74f8a | 123 | TH1F* contSecMC[6]; |
88d02ae4 | 124 | |
125 | TH1F* contfit[12]; | |
126 | TH1F* contWDfit[12]; | |
127 | TH1F* contMatfit[12]; | |
128 | TH1F* primaryfit[12]; | |
129 | ||
130 | ||
131 | ||
132 | TH1F* spectra[6]; | |
133 | TH1F* spectraLeonardo[6]; | |
134 | ||
135 | TH1F* corrLeonardo[6]; | |
136 | //GetSpectra(managerdata,rawspectradata,true); | |
137 | //GetSpectra(managermc,rawspectramc,true,true); | |
138 | ||
139 | GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraMC",managermc,rawspectramc,dcacutxy); | |
140 | GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraDATA",managerdata,rawspectradata,dcacutxy); | |
141 | GetPtHistFromPtDCAhisto("hHistPtRecTruePrimary","eff",managermc,eff,dcacutxy); | |
142 | GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy); | |
143 | GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy); | |
144 | GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy); | |
2245e1bd | 145 | GetPtHistFromPtDCAhisto("hHistPtRecSigmaPrimary","conPIDprimary",managermc,contPIDpri,dcacutxy); |
146 | ||
88d02ae4 | 147 | |
b3ea73e1 | 148 | Double_t neventsmcall = 1 ; //if loop over MC is done after or befor events cuts this will be changed |
149 | Double_t neventsdata = 1; | |
150 | Double_t neventsmc = 1; | |
151 | ||
3a71f081 | 152 | //Normaliztion of MCtruth depends if the loop was done after of before ESD event cuts. |
153 | //In currect code this cannot be check on the level of macro. | |
154 | //If the loop was done before MC should be done to all processed events (NumberOfProcessedEvents()) | |
155 | //If loop was done after MC should be normalized to all accepted events (NumberOfEvents()) | |
156 | // The option one will be alaways use. | |
157 | ||
158 | neventsmcall= ecutsmc->NumberOfProcessedEvents(); | |
b3ea73e1 | 159 | if(options&knormalizationtoeventspassingPhySel) |
160 | { | |
3a71f081 | 161 | //neventsmcall= ecutsmc->NumberOfProcessedEvents(); |
b3ea73e1 | 162 | neventsdata=ecutsdata->NumberOfPhysSelEvents(); |
163 | neventsmc=ecutsmc->NumberOfPhysSelEvents(); | |
164 | } | |
3a71f081 | 165 | else if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC)) |
166 | { | |
167 | neventsdata=Normaliztionwithbin0integrals(options); | |
168 | neventsmc=ecutsmc->NumberOfPhysSelEvents(); | |
169 | } | |
b3ea73e1 | 170 | else |
171 | { | |
50c532c4 | 172 | neventsdata=ecutsdata->NumberOfEvents(); //number of accepted events |
b3ea73e1 | 173 | neventsmc=ecutsmc->NumberOfEvents(); |
174 | neventsmcall= ecutsmc->NumberOfEvents(); | |
175 | ||
b3ea73e1 | 176 | } |
88d02ae4 | 177 | GetMCTruth(MCTruth); |
50c532c4 | 178 | cout<<neventsdata<<" Events"<<endl; |
88d02ae4 | 179 | |
180 | ||
181 | TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone(); | |
182 | allgen->SetName("AllGen"); | |
183 | TH1F* allrecMC=GetOneHistFromPtDCAhisto("hHistPtRec","rawallMC",managermc,dcacutxy); | |
184 | TH1F* alleff=GetOneHistFromPtDCAhisto("hHistPtRecPrimary","effall",managermc,dcacutxy); | |
185 | TH1F* allrecdata=GetOneHistFromPtDCAhisto("hHistPtRec","rawalldata",managerdata,dcacutxy); | |
186 | ||
187 | ||
188 | ||
189 | ||
190 | TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch"); | |
a4ca64e3 | 191 | spectraall->SetTitle("recNch"); |
88d02ae4 | 192 | TH1F* contall=(TH1F*)allrecMC->Clone("contall"); |
a4ca64e3 | 193 | contall->SetTitle("contall"); |
50c532c4 | 194 | //contall->Add(alleff,-1); |
195 | SubHistWithFullCorr(contall,alleff); | |
88d02ae4 | 196 | alleff->Divide(alleff,allgen,1,1,"B"); |
b3ea73e1 | 197 | contall->Divide(contall,allrecMC,1,1,"B"); |
88d02ae4 | 198 | |
199 | GetCorrectedSpectra(spectraall,allrecdata,alleff,contall); | |
200 | Divideby2pipt(spectraall); | |
201 | ||
202 | allrecdata->Scale(1./neventsdata,"width"); | |
203 | allgen->Scale(1./neventsmcall,"width"); | |
204 | allrecMC->Scale(1./neventsmc,"width"); | |
205 | spectraall->Scale(1./neventsdata,"width"); | |
206 | ||
207 | ||
208 | lout->Add(allgen); | |
209 | lout->Add(allrecMC); | |
210 | lout->Add(alleff); | |
211 | lout->Add(allrecdata); | |
212 | lout->Add(spectraall); | |
1e6bff53 | 213 | lout->Add(contall); |
214 | ||
88d02ae4 | 215 | for (int i=0;i<6;i++) |
216 | { | |
217 | ||
218 | ||
219 | TString tmpname(rawspectramc[i]->GetTitle()); | |
220 | tmpname.ReplaceAll("SpectraMC","%s"); | |
a4ca64e3 | 221 | contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC")); |
222 | contallMC[i]->SetTitle(Form(tmpname.Data(),"contallMC")); | |
88d02ae4 | 223 | contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit")); |
a4ca64e3 | 224 | contfit[i]->SetTitle(Form(tmpname.Data(),"contfit")); |
88d02ae4 | 225 | contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit")); |
a4ca64e3 | 226 | contWDfit[i]->SetTitle(Form(tmpname.Data(),"contWDfit")); |
88d02ae4 | 227 | contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit")); |
a4ca64e3 | 228 | contMatfit[i]->SetTitle(Form(tmpname.Data(),"contMatfit")); |
88d02ae4 | 229 | primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit")); |
a4ca64e3 | 230 | primaryfit[i]->SetTitle(Form(tmpname.Data(),"primaryfit")); |
88d02ae4 | 231 | contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC")); |
a4ca64e3 | 232 | contfit[i+6]->SetTitle(Form(tmpname.Data(),"contfitonMC")); |
88d02ae4 | 233 | contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC")); |
a4ca64e3 | 234 | contWDfit[i+6]->SetTitle(Form(tmpname.Data(),"contWDfitonMC")); |
88d02ae4 | 235 | contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC")); |
a4ca64e3 | 236 | contMatfit[i+6]->SetTitle(Form(tmpname.Data(),"contMatfitonMC")); |
88d02ae4 | 237 | primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC")); |
a4ca64e3 | 238 | primaryfit[i+6]->SetTitle(Form(tmpname.Data(),"primaryfitMC")); |
88d02ae4 | 239 | |
b6a74f8a | 240 | |
241 | ||
88d02ae4 | 242 | contfit[i]->Reset(); |
243 | contWDfit[i]->Reset(); | |
244 | contMatfit[i]->Reset(); | |
245 | primaryfit[i]->Reset(); | |
246 | ||
247 | ||
248 | contfit[i+6]->Reset(); | |
249 | contWDfit[i+6]->Reset(); | |
250 | contMatfit[i+6]->Reset(); | |
251 | primaryfit[i+6]->Reset(); | |
252 | ||
253 | SetBintoOne(primaryfit[i]); | |
254 | SetBintoOne(primaryfit[i+6]); | |
255 | spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal")); | |
a4ca64e3 | 256 | spectra[i]->SetTitle(Form(tmpname.Data(),"SpectraFinal")); |
88d02ae4 | 257 | spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo")); |
a4ca64e3 | 258 | spectraLeonardo[i]->SetTitle(Form(tmpname.Data(),"SpectraFinalLeonardo")); |
259 | ||
88d02ae4 | 260 | corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo")); |
a4ca64e3 | 261 | corrLeonardo[i]->SetTitle(Form(tmpname.Data(),"CorrFactLeonardo")); |
88d02ae4 | 262 | corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B"); |
263 | ||
264 | ||
265 | ||
50c532c4 | 266 | //contallMC[i]->Add(eff[i],-1.0); |
267 | SubHistWithFullCorr(contallMC[i],eff[i]); | |
268 | //RecomputeErrors(contallMC[i]); | |
88d02ae4 | 269 | contallMC[i]->Sumw2(); |
270 | contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B"); | |
2245e1bd | 271 | |
272 | // contamintaion from PID but only primaries | |
50c532c4 | 273 | //contPIDpri[i]->Add(eff[i],-1.0); |
274 | SubHistWithFullCorr(contPIDpri[i],eff[i]); | |
275 | ||
276 | //RecomputeErrors(contPIDpri[i]); | |
2245e1bd | 277 | contPIDpri[i]->Divide(contPIDpri[i],rawspectramc[i],1,1,"B"); |
2245e1bd | 278 | |
88d02ae4 | 279 | eff[i]->Divide(eff[i],MCTruth[i],1,1,"B"); |
280 | ||
281 | ||
282 | contPID[i]->Sumw2(); | |
283 | rawspectramc[i]->Sumw2(); | |
50c532c4 | 284 | //contPID[i]->Add(contPID[i],rawspectramc[i],-1,1); |
285 | SubHistWithFullCorr(contPID[i],rawspectramc[i]); | |
286 | contPID[i]->Scale(-1.0); | |
287 | ||
288 | //RecomputeErrors(contPID[i]); | |
88d02ae4 | 289 | contPID[i]->ResetStats(); |
290 | contPID[i]->Sumw2(); | |
291 | contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B"); | |
2245e1bd | 292 | |
293 | if(options&kuseprimaryPIDcont) | |
294 | confinal[i]=(TH1F*)contPIDpri[i]->Clone(Form(tmpname.Data(),"confinal")); | |
295 | else | |
296 | confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal")); | |
a4ca64e3 | 297 | confinal[i]->SetTitle(Form(tmpname.Data(),"confinal")); |
88d02ae4 | 298 | |
b6a74f8a | 299 | contSecMC[i]=(TH1F*)contWD[i]->Clone(Form(tmpname.Data(),"conSecMC")); |
300 | contSecMC[i]->Add(contMat[i]); | |
88d02ae4 | 301 | contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B"); |
302 | contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B"); | |
b6a74f8a | 303 | contSecMC[i]->Divide(contSecMC[i],rawspectramc[i],1,1,"B"); |
88d02ae4 | 304 | |
305 | ||
306 | rawspectradata[i]->Scale(1./neventsdata,"width"); | |
307 | rawspectramc[i]->Scale(1./neventsmc,"width"); | |
308 | MCTruth[i]->Scale(1./neventsmcall,"width"); | |
309 | spectraLeonardo[i]->Scale(1./neventsdata,"width"); | |
310 | ||
311 | ||
312 | ||
313 | lout->Add(rawspectradata[i]); | |
314 | lout->Add(rawspectramc[i]); | |
315 | lout->Add(MCTruth[i]); | |
316 | lout->Add(eff[i]); | |
317 | lout->Add(contallMC[i]); | |
318 | lout->Add(contPID[i]); | |
319 | lout->Add(contWD[i]); | |
320 | lout->Add(contMat[i]); | |
321 | lout->Add(contfit[i]); | |
322 | lout->Add(contWDfit[i]); | |
323 | lout->Add(contMatfit[i]); | |
1e6bff53 | 324 | lout->Add(primaryfit[i]); |
88d02ae4 | 325 | lout->Add(contfit[i+6]); |
326 | lout->Add(contWDfit[i+6]); | |
327 | lout->Add(contMatfit[i+6]); | |
1e6bff53 | 328 | lout->Add(primaryfit[i+6]); |
88d02ae4 | 329 | lout->Add(spectra[i]); |
330 | lout->Add(spectraLeonardo[i]); | |
331 | lout->Add(confinal[i]); | |
2245e1bd | 332 | lout->Add(contPIDpri[i]); |
b6a74f8a | 333 | lout->Add(contSecMC[i]); |
88d02ae4 | 334 | } |
2245e1bd | 335 | outdate.ReplaceAll("/","_"); |
a4ca64e3 | 336 | configfile.ReplaceAll(".","_"); |
337 | TFile* fout=0x0; | |
338 | if(configfile.Length()>0&&(options&kuserangeonfigfile)) | |
339 | fout=new TFile(Form("./results/ResMY_%s_%s_%#X_%s.root",outnamemc.Data(),outdate.Data(),options,configfile.Data()),"RECREATE"); | |
340 | else | |
341 | fout=new TFile(Form("./results/ResMY_%s_%s_%#X.root",outnamemc.Data(),outdate.Data(),options),"RECREATE"); | |
b3ea73e1 | 342 | if (options&kdodca) |
88d02ae4 | 343 | DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit); |
344 | for (int i=0;i<6;i++) | |
345 | { | |
b3ea73e1 | 346 | if(options&kdodca) |
88d02ae4 | 347 | { |
2245e1bd | 348 | if((options&kuseprimaryPIDcont)||(i!=1&&i!=4)) //if we do not use cont PId only for primary and this is a kaon that do not use fit |
349 | confinal[i]->Add(contfit[i]); | |
88d02ae4 | 350 | GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]); |
351 | } | |
352 | else | |
353 | { | |
354 | GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]); | |
355 | } | |
356 | GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]); | |
b6a74f8a | 357 | if(options&kskipconcutonspectra) |
358 | continue; | |
2245e1bd | 359 | if(options&kuseprimaryPIDcont) |
360 | { | |
361 | CleanHisto(spectra[i],-1,100,contPIDpri[i]); | |
362 | CleanHisto(spectraLeonardo[i],-1,100,contPIDpri[i]); | |
363 | } | |
364 | else | |
365 | { | |
366 | CleanHisto(spectra[i],-1,100,contPID[i]); | |
367 | CleanHisto(spectraLeonardo[i],-1,100,contPID[i]); | |
cd28560b | 368 | } |
369 | // Apply correction for wrongly simulated TOF signal in MC | |
370 | if(options&kuseTOFcorrforPIDsignalmatching) | |
371 | TOFPIDsignalmatchingApply(spectra[i],TOFPIDsignalmatching[i%3]); | |
372 | ||
88d02ae4 | 373 | } |
374 | ||
b3ea73e1 | 375 | GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options); |
376 | GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options); | |
ab0298c1 | 377 | |
378 | Double_t ycut=tcutsdata->GetY(); | |
379 | Double_t etacut=tcutsdata->GetEtaMax()-tcutsdata->GetEtaMin(); | |
380 | if(etacut<0.00001) | |
381 | { | |
382 | cout<<"using eta window 1.6"<<endl; | |
383 | etacut=1.6; | |
384 | ||
385 | } | |
cd28560b | 386 | |
ab0298c1 | 387 | TH1F* allch=GetSumAllCh(spectra,mass,etacut); |
88d02ae4 | 388 | lout->Add(allch); |
1af4558e | 389 | if(options&kuseTOFmatchingcorrection) |
390 | { | |
391 | MatchingTOFEff(spectra,lout); | |
392 | MatchingTOFEff(spectraLeonardo); | |
393 | TOFMatchingForNch(spectraall); | |
394 | } | |
88d02ae4 | 395 | // lout->ls(); |
396 | fout->cd(); | |
397 | TList* listqa=new TList(); | |
33d374d5 | 398 | TList* canvaslist=new TList(); |
b3ea73e1 | 399 | Float_t vertexcorrection=1.0; |
3a71f081 | 400 | Float_t corrbadchunksvtx=1.0; |
401 | if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC)) | |
402 | corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,0); | |
403 | else | |
404 | corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,1); | |
b3ea73e1 | 405 | if (options&kveretxcorrectionandbadchunkscorr) |
406 | vertexcorrection=corrbadchunksvtx; | |
88d02ae4 | 407 | cout<<" VTX corr="<<vertexcorrection<<endl; |
b3ea73e1 | 408 | if(TMath::Abs(ycut)>0.0) |
409 | vertexcorrection=vertexcorrection/(2.0*ycut); | |
88d02ae4 | 410 | for (int i=0;i<6;i++) |
411 | { | |
412 | spectra[i]->Scale(vertexcorrection); | |
413 | spectraLeonardo[i]->Scale(vertexcorrection); | |
890fd855 | 414 | if(TMath::Abs(ycut)>0.0) |
415 | { | |
416 | rawspectradata[i]->Scale(1.0/(2.0*ycut)); | |
417 | rawspectramc[i]->Scale(1.0/(2.0*ycut)); | |
418 | MCTruth[i]->Scale(1.0/(2.0*ycut)); | |
419 | } | |
88d02ae4 | 420 | } |
421 | allch->Scale(vertexcorrection); | |
ab0298c1 | 422 | spectraall->Scale(vertexcorrection/etacut); |
88d02ae4 | 423 | |
424 | //spectraall->Scale(1.0/1.6); | |
425 | lout->Write("output",TObject::kSingleKey); | |
426 | listqa->Write("outputQA",TObject::kSingleKey); | |
33d374d5 | 427 | canvaslist->Write("outputcanvas",TObject::kSingleKey); |
428 | ||
88d02ae4 | 429 | fout->Close(); |
3a71f081 | 430 | //Normaliztionwithbin0integrals(); |
88d02ae4 | 431 | |
432 | } | |
433 | ||
0c1898ba | 434 | Bool_t OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcasdata) |
88d02ae4 | 435 | { |
0c1898ba | 436 | |
88d02ae4 | 437 | |
438 | TString nameFile = Form("./%s/AnalysisResults%s.root",dirname.Data(),(mcflag?"MC":"DATA")); | |
439 | TFile *file = TFile::Open(nameFile.Data()); | |
440 | if(!file) | |
441 | { | |
442 | cout<<"no file"<<endl; | |
443 | return false; | |
444 | } | |
445 | TString sname=Form("OutputBothSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data()); | |
0c1898ba | 446 | if(mcasdata) |
447 | { | |
448 | cout<<"using MC as data "<<endl; | |
449 | sname=Form("OutputBothSpectraTask_%s_%s","MC",outputname.Data()); | |
450 | } | |
88d02ae4 | 451 | file->ls(); |
452 | TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data()); | |
453 | if(!dir) | |
454 | { | |
ab13cf18 | 455 | // cout<<"no dir "<<sname.Data()<<endl; |
456 | if(mcasdata) | |
0c1898ba | 457 | { |
458 | cout<<"using MC as data "<<endl; | |
459 | sname=Form("OutputAODSpectraTask_%s_%s","MC",outputname.Data()); | |
460 | } | |
461 | else | |
462 | sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data()); | |
88d02ae4 | 463 | // cout<<"trying "<<sname.Data()<<endl; |
464 | dir=(TDirectoryFile*)file->Get(sname.Data()); | |
465 | if(!dir) | |
466 | { | |
467 | cout<<"no dir "<<sname.Data()<<endl; | |
468 | return false; | |
469 | } | |
470 | } | |
471 | cout << " -- Info about " <<(mcflag?"MC":"DATA") <<" -- "<< endl; | |
472 | if(mcflag) | |
473 | { | |
474 | managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos"); | |
475 | ecutsmc = (AliSpectraBothEventCuts*) dir->Get("Event Cuts"); | |
476 | tcutsmc = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts"); | |
477 | ecutsmc->PrintCuts(); | |
478 | tcutsmc->PrintCuts(); | |
479 | if(!managermc||!ecutsmc||!tcutsmc) | |
480 | return false; | |
b6a74f8a | 481 | if(managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsmc->GetHistoCuts()->GetBinContent(3)) |
482 | cout<<"Please check MC file possible problem with merging"<<" "<<managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsmc->GetHistoCuts()->GetBinContent(3)<<endl; | |
88d02ae4 | 483 | } |
484 | else | |
485 | { | |
486 | managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos"); | |
487 | ecutsdata = (AliSpectraBothEventCuts*) dir->Get("Event Cuts"); | |
488 | tcutsdata = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts"); | |
489 | ecutsdata->PrintCuts(); | |
490 | tcutsdata->PrintCuts(); | |
491 | if(!managerdata||!ecutsdata||!tcutsdata) | |
492 | return false; | |
b6a74f8a | 493 | if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3)) |
494 | cout<<"Please check DATA file possible problem with merging"<<" "<<anagerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl; | |
a4ca64e3 | 495 | |
88d02ae4 | 496 | } |
497 | return true; | |
498 | } | |
499 | ||
500 | void GetMCTruth(TH1F** MCTruth) | |
501 | { | |
502 | for(Int_t icharge=0;icharge<2;icharge++) | |
503 | { | |
504 | for(Int_t ipart=0;ipart<3;ipart++) | |
505 | { | |
506 | Int_t index=ipart+3*icharge; | |
507 | TString hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); | |
508 | MCTruth[index]=(TH1F*)((TH1F*)managermc->GetPtHistogram1D(hname.Data(),1,1))->Clone(); | |
509 | MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); | |
510 | MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); | |
511 | MCTruth[index]->Sumw2(); | |
512 | } | |
513 | } | |
514 | } | |
515 | ||
516 | TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHistoManager* hman,TFormula* dcacutxy) | |
517 | { | |
518 | histo =(TH1F*)((TH1F*) hman->GetPtHistogram1D(name.Data(),-1,-1))->Clone(); | |
519 | histo->SetName(hnameout.Data()); | |
520 | histo->SetTitle(hnameout.Data()); | |
521 | ||
522 | if(dcacutxy) | |
523 | { | |
524 | for(int ibin=1;ibin<histo->GetNbinsX();ibin++) | |
525 | { | |
526 | Double_t lowedge=histo->GetBinLowEdge(ibin); | |
527 | Float_t cut=dcacutxy->Eval(lowedge); | |
50c532c4 | 528 | TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge); |
1af4558e | 529 | // Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut)); |
50c532c4 | 530 | Float_t testyield=0.0; |
531 | Float_t testerror=0.0; | |
532 | for (int itest=dcahist->GetXaxis()->FindBin(-1.0*cut);itest<=dcahist->GetXaxis()->FindBin(cut);itest++) | |
533 | { | |
534 | testyield+=dcahist->GetBinContent(itest); | |
535 | testerror+=dcahist->GetBinError(itest)*dcahist->GetBinError(itest); | |
536 | } | |
1af4558e | 537 | // cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl; |
538 | // cout<<"test dca "<<lowedge<<" "<<dcacutxy->Eval(lowedge)<<" "<<dcacutxy->Eval(histo->GetXaxis()->GetBinUpEdge(ibin))<<" "<<dcahist->GetBinLowEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<" "<<dcahist->GetXaxis()->GetBinUpEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<endl; | |
50c532c4 | 539 | |
1af4558e | 540 | // cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl; |
50c532c4 | 541 | |
542 | //histo->SetBinContent(ibin,inyield); | |
543 | //histo->SetBinError(ibin,TMath::Sqrt(inyield)); | |
544 | histo->SetBinContent(ibin,testyield); | |
545 | histo->SetBinError(ibin,TMath::Sqrt(testerror)); | |
546 | ||
88d02ae4 | 547 | } |
548 | } | |
549 | histo->Sumw2(); | |
550 | return histo; | |
551 | } | |
552 | ||
553 | ||
554 | ||
555 | ||
556 | void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy) | |
557 | { | |
1af4558e | 558 | //Float_t min[3]={0.3,0.3,0.45}; |
559 | //Float_t max[3]={1.5,1.2,2.2}; | |
88d02ae4 | 560 | for(Int_t icharge=0;icharge<2;icharge++) |
561 | { | |
562 | for(Int_t ipart=0;ipart<3;ipart++) | |
563 | { | |
564 | Int_t index=ipart+3*icharge; | |
a4ca64e3 | 565 | Printf("Getting %s",hnamein.Data()); |
88d02ae4 | 566 | TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()); |
567 | TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data()); | |
568 | ||
569 | ||
570 | histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy); | |
1af4558e | 571 | CleanHisto(histo[index],minRanges[ipart],maxRanges[ipart]); |
88d02ae4 | 572 | } |
573 | } | |
574 | } | |
575 | void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0) | |
576 | { | |
577 | for (int i=0;i<=h->GetNbinsX();i++) | |
578 | { | |
579 | if(h->GetXaxis()->GetBinCenter(i)<minV||h->GetXaxis()->GetBinCenter(i)>maxV) | |
580 | { | |
581 | h->SetBinContent(i,0); | |
582 | h->SetBinError(i,0); | |
583 | } | |
584 | if(contpid) | |
585 | { | |
1af4558e | 586 | if(contpid->GetBinContent(i)>fMaxContaminationPIDMC) |
88d02ae4 | 587 | { |
588 | h->SetBinContent(i,0); | |
589 | h->SetBinError(i,0); | |
590 | } | |
591 | } | |
592 | } | |
593 | } | |
594 | ||
595 | ||
596 | void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary) | |
597 | { | |
598 | printf("\n\n-> DCA Correction"); | |
a4ca64e3 | 599 | |
2245e1bd | 600 | |
88d02ae4 | 601 | Printf("\DCACorr"); |
88d02ae4 | 602 | TString sample[2]={"data","mc"}; |
603 | ofstream debug("debugDCA.txt"); | |
604 | TList* listofdcafits=new TList(); | |
605 | for(Int_t icharge=0;icharge<2;icharge++) | |
606 | { | |
607 | for(Int_t ipart=0;ipart<3;ipart++) | |
608 | { | |
609 | Int_t index=ipart+3*icharge; | |
610 | for(Int_t isample=0;isample<2;isample++) | |
611 | { | |
612 | ||
88d02ae4 | 613 | for(Int_t ibin_data=7;ibin_data<40;ibin_data++) |
614 | { | |
615 | ||
616 | Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data); | |
617 | Double_t binwidth=hcon[index]->GetBinWidth(ibin_data); | |
8cec8e63 | 618 | debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl; |
619 | if(fun) | |
620 | { | |
621 | CutRange[1]=fun->Eval(lowedge); | |
622 | CutRange[0]=-1.0*CutRange[1]; | |
623 | } | |
88d02ae4 | 624 | debug<<"cut "<<CutRange[1]<<" "<<CutRange[0]<<endl; |
2245e1bd | 625 | Short_t fitsettings=DCAfitsettings(lowedge+0.5*binwidth,index); |
626 | debug<<"settings "<< fitsettings<<endl; | |
627 | if(fitsettings==0) | |
628 | continue; | |
629 | ||
88d02ae4 | 630 | TCanvas *cDCA=new TCanvas(Form("cDCA%d%s%s%sbin%d",index,sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),ibin_data),Form("cDCA%d%s%s%sbin%d",index,sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),ibin_data),1700,1500); |
b6a74f8a | 631 | TLegend* Leg1 = new TLegend(0.6,0.6,0.85,0.85,"","NDC"); |
632 | Leg1->SetFillStyle(kFALSE); | |
633 | Leg1->SetLineColor(kWhite); | |
634 | Leg1->SetBorderSize(0); | |
635 | ||
636 | ||
88d02ae4 | 637 | if(isample==0) |
638 | TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
639 | if(isample==1) | |
640 | TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
641 | debug<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl; | |
642 | TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
643 | TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
644 | TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone(); | |
645 | Double_t minentries=1; | |
2245e1bd | 646 | debug<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl; |
647 | debug<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl; | |
648 | if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || ((fitsettings&0x1)&&hmc2->GetEntries()<=minentries) || ((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)) | |
88d02ae4 | 649 | continue; |
88d02ae4 | 650 | hToFit->Sumw2(); |
651 | hmc1->Sumw2(); | |
652 | hmc2->Sumw2(); | |
653 | hmc3->Sumw2(); | |
654 | ||
655 | Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0}; | |
2245e1bd | 656 | Int_t binCutRange[]={hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])}; |
657 | ||
88d02ae4 | 658 | |
88d02ae4 | 659 | if(hmc3->GetNbinsX()>300) |
660 | { | |
661 | ||
2245e1bd | 662 | corrforrebinning[0]=hToFit->Integral(binCutRange[0],binCutRange[1]); |
663 | corrforrebinning[1]=hmc1->Integral(binCutRange[0],binCutRange[1]); | |
664 | corrforrebinning[2]=hmc2->Integral(binCutRange[0],binCutRange[1]); | |
665 | corrforrebinning[3]=hmc3->Integral(binCutRange[0],binCutRange[1]); | |
88d02ae4 | 666 | |
667 | hToFit->Rebin(30); | |
668 | hmc1->Rebin(30); | |
669 | hmc2->Rebin(30); | |
670 | hmc3->Rebin(30); | |
2245e1bd | 671 | |
672 | binCutRange[0]=hmc1->GetXaxis()->FindBin(CutRange[0]); | |
673 | binCutRange[1]=hmc1->GetXaxis()->FindBin(CutRange[1]); | |
674 | ||
675 | ||
676 | //after rebbing we lose resolution of the dca this correction also us to do obtain inside used dca | |
677 | ||
678 | if(hToFit->Integral(binCutRange[0],binCutRange[1])>0.0) | |
679 | corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(binCutRange[0],binCutRange[1]); | |
88d02ae4 | 680 | else |
681 | corrforrebinning[0]=1.0; | |
2245e1bd | 682 | if(hmc1->Integral(binCutRange[0],binCutRange[1])>0.0) |
683 | corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(binCutRange[0],binCutRange[1]); | |
88d02ae4 | 684 | else |
685 | corrforrebinning[1]=1.0; | |
2245e1bd | 686 | if(hmc2->Integral(binCutRange[0],binCutRange[1])>0.0) |
687 | corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(binCutRange[0],binCutRange[1]); | |
88d02ae4 | 688 | else |
689 | corrforrebinning[2]=1.0; | |
2245e1bd | 690 | if(hmc3->Integral(binCutRange[0],binCutRange[1])>0.0) |
691 | corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(binCutRange[0],binCutRange[1]); | |
88d02ae4 | 692 | else |
693 | corrforrebinning[3]=1.0; | |
694 | ||
88d02ae4 | 695 | |
696 | ||
697 | } | |
698 | ||
699 | cDCA->cd(); | |
700 | gPad->SetGridy(); | |
701 | gPad->SetGridx(); | |
702 | gPad->SetLogy(); | |
703 | ||
704 | TObjArray *mc=0x0; | |
1af4558e | 705 | Int_t Npar=3; |
2245e1bd | 706 | if(fitsettings==3) |
88d02ae4 | 707 | mc = new TObjArray(3); // MC histograms are put in this array |
708 | else | |
1af4558e | 709 | { |
88d02ae4 | 710 | mc = new TObjArray(2); |
1af4558e | 711 | Npar=2; |
712 | } | |
2245e1bd | 713 | |
88d02ae4 | 714 | mc->Add(hmc1); |
2245e1bd | 715 | if(fitsettings&0x1) |
716 | mc->Add(hmc2); | |
717 | if(fitsettings&0x2) | |
88d02ae4 | 718 | mc->Add(hmc3); |
719 | TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise | |
720 | fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1 | |
1af4558e | 721 | Double_t defaultStep = 0.01; |
722 | for (int iparm = 0; iparm < Npar; ++iparm) | |
723 | { | |
724 | TString name("frac"); name += iparm; | |
725 | if(iparm==0) | |
726 | fit->GetFitter()->SetParameter(iparm,name.Data(),0.85,0.01,0.0,1.0); | |
727 | else if (Npar==2) | |
728 | fit->GetFitter()->SetParameter(iparm,name.Data(),0.15,0.01,0.0,1.0); | |
729 | else | |
730 | fit->GetFitter()->SetParameter(iparm,name.Data(),0.075,0.01,0.0,1.0); | |
731 | ||
732 | } | |
2245e1bd | 733 | if(fitsettings&0x1) |
734 | fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1 | |
735 | if(fitsettings&0x2) | |
736 | fit->Constrain(1+(fitsettings&0x1),0.0,1.0); // constrain fraction 1 to be between 0 and 1 | |
737 | ||
738 | Int_t binFitRange[]={hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])}; | |
739 | fit->SetRangeX(binFitRange[0],binFitRange[1]); | |
740 | hToFit->GetXaxis()->SetRange(binFitRange[0],binFitRange[1]); | |
88d02ae4 | 741 | hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge)); |
742 | Int_t status = fit->Fit(); // perform the fit | |
743 | cout << "fit status: " << status << endl; | |
744 | debug<<"fit status: " << status << endl; | |
745 | ||
2245e1bd | 746 | if (status == 0) |
747 | { | |
748 | Double_t v1=0.0,v2=0.0,v3=0.0; | |
749 | Double_t ev1=0.0,ev2=0.0,ev3=0.0; | |
750 | Double_t cov=0.0; | |
751 | ||
752 | // check on fit status | |
88d02ae4 | 753 | TH1F* result = (TH1F*) fit->GetPlot(); |
754 | TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0); | |
2245e1bd | 755 | |
756 | TH1F* secStMCPred=0X0; | |
88d02ae4 | 757 | |
758 | TH1F* secMCPred=0x0; | |
88d02ae4 | 759 | |
88d02ae4 | 760 | //first method, use directly the fit result |
761 | fit->GetResult(0,v1,ev1); | |
2245e1bd | 762 | |
763 | if(fitsettings&0x1) | |
764 | { | |
765 | fit->GetResult(1,v2,ev2); | |
766 | secStMCPred=(TH1F*)fit->GetMCPrediction(1); | |
767 | ||
768 | } | |
769 | if(fitsettings&0x2) | |
88d02ae4 | 770 | { |
2245e1bd | 771 | fit->GetResult(1+(fitsettings&0x1),v3,ev3); |
772 | secMCPred=(TH1F*)fit->GetMCPrediction(1+(fitsettings&0x1)); | |
773 | if(fitsettings&0x1) | |
774 | cov=fit->GetFitter()->GetCovarianceMatrixElement(1,2); | |
775 | ||
776 | ||
88d02ae4 | 777 | } |
2245e1bd | 778 | |
88d02ae4 | 779 | |
2245e1bd | 780 | debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<" "<<cov<<endl; |
781 | ||
782 | // becuase dca cut range is not a fit range the results from TFractionFitter should be rescale | |
783 | // This can be done in two ways or use input histograms or output histograms | |
784 | // The difference between those two methods should be on the level of statistical error | |
785 | // I use output histograms | |
88d02ae4 | 786 | |
2245e1bd | 787 | // Method 1 input histo |
788 | ||
789 | Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(binFitRange[0],binFitRange[1]); | |
790 | normalizationdata*=corrforrebinning[0]; | |
791 | ||
792 | Float_t normalizationmc1=(hmc1->Integral(binCutRange[0],binCutRange[1])/hmc1->Integral(binFitRange[0],binFitRange[1]))/normalizationdata; | |
793 | Float_t normalizationmc2=0.0; | |
794 | if(fitsettings&0x1) | |
795 | normalizationmc2=(hmc2->Integral(binCutRange[0],binCutRange[1])/hmc2->Integral(binFitRange[0],binFitRange[1]))/normalizationdata; | |
796 | Float_t normalizationmc3=0.0; | |
797 | if(fitsettings&0x2) | |
798 | normalizationmc3=(hmc3->Integral(binCutRange[0],binCutRange[1])/hmc3->Integral(binFitRange[0],binFitRange[1]))/normalizationdata; | |
799 | ||
800 | normalizationmc1*=corrforrebinning[1]; | |
801 | normalizationmc2*=corrforrebinning[2]; | |
802 | normalizationmc3*=corrforrebinning[3]; | |
803 | ||
88d02ae4 | 804 | debug<<"After Nor"<<endl; |
805 | debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl; | |
806 | debug<<1.0-v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2+v3*normalizationmc3<<" "<<TMath::Sqrt(ev2*ev2*normalizationmc2*normalizationmc2+ev3*ev3*normalizationmc3*normalizationmc3+cov*normalizationmc3*normalizationmc2)<<endl; | |
8cec8e63 | 807 | debug<<"addtional info"<<endl; |
8cec8e63 | 808 | |
809 | ||
2245e1bd | 810 | //Method 2 output histo |
8cec8e63 | 811 | |
2245e1bd | 812 | Float_t normalizationdata1=result->Integral(binCutRange[0],binCutRange[1])/result->Integral(binFitRange[0],binFitRange[1]); |
88d02ae4 | 813 | |
814 | ||
a4ca64e3 | 815 | // if the cut range is bigger the fit range we should calculate the normalization factor for data using the data histogram |
816 | // because result histogram has entries only in fits range | |
817 | if(FitRange[0]>CutRange[0]||FitRange[1]<CutRange[1]) | |
818 | normalizationdata1=normalizationdata; | |
819 | ||
88d02ae4 | 820 | normalizationdata1*=corrforrebinning[0]; |
821 | ||
822 | ||
2245e1bd | 823 | Float_t normalizationmc11=(PrimMCPred->Integral(binCutRange[0],binCutRange[1])/PrimMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1; |
824 | Float_t normalizationmc21=0.0; | |
825 | if(fitsettings&0x1) | |
826 | normalizationmc21=(secStMCPred->Integral(binCutRange[0],binCutRange[1])/secStMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1; | |
827 | Float_t normalizationmc31=0.0; | |
828 | if(fitsettings&0x2) | |
829 | normalizationmc31=(secMCPred->Integral(binCutRange[0],binCutRange[1])/secMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1; | |
88d02ae4 | 830 | |
831 | normalizationmc11*=corrforrebinning[1]; | |
832 | normalizationmc21*=corrforrebinning[2]; | |
833 | normalizationmc31*=corrforrebinning[3]; | |
834 | ||
2245e1bd | 835 | debug<<"After Nor 2"<<endl; |
88d02ae4 | 836 | debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl; |
837 | ||
838 | debug<<1.0-v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21+v3*normalizationmc31<<" "<<TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21)<<endl; | |
839 | ||
88d02ae4 | 840 | |
841 | hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21); | |
842 | hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21); | |
843 | hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31); | |
844 | hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31); | |
845 | hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11); | |
1e6bff53 | 846 | hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11); |
2245e1bd | 847 | hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31); |
848 | hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21)); | |
88d02ae4 | 849 | |
850 | ||
851 | ||
852 | //Drawing section | |
853 | result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1]))); | |
854 | hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]))); | |
855 | PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1]))); | |
2245e1bd | 856 | |
88d02ae4 | 857 | hToFit->SetMinimum(0.0001); |
858 | hToFit->DrawClone("E1x0"); | |
859 | result->SetTitle("Fit result"); | |
2245e1bd | 860 | result->SetLineColor(kBlack); |
b6a74f8a | 861 | Leg1->AddEntry(result,"result","lp"); |
88d02ae4 | 862 | result->DrawClone("lhistsame"); |
2245e1bd | 863 | |
864 | PrimMCPred->SetLineColor(kGreen+2); | |
b6a74f8a | 865 | PrimMCPred->SetLineStyle(2); |
866 | PrimMCPred->SetLineWidth(3.0); | |
867 | Leg1->AddEntry(PrimMCPred,"Prmi.","l"); | |
88d02ae4 | 868 | PrimMCPred->DrawClone("lhistsame"); |
2245e1bd | 869 | if(fitsettings&0x1) |
88d02ae4 | 870 | { |
2245e1bd | 871 | |
872 | secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1]))); | |
873 | secStMCPred->SetLineColor(kRed); | |
b6a74f8a | 874 | secStMCPred->SetLineWidth(3.0); |
875 | ||
876 | secStMCPred->SetLineStyle(3); | |
877 | Leg1->AddEntry(secStMCPred,"Sec.WD","l"); | |
2245e1bd | 878 | secStMCPred->DrawClone("lhistsame"); |
879 | ||
880 | } | |
881 | if(fitsettings&0x2) | |
882 | { | |
883 | ||
884 | secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1]))); | |
b6a74f8a | 885 | secMCPred->SetLineColor(kBlue); |
886 | secMCPred->SetLineWidth(3.0); | |
887 | ||
888 | secMCPred->SetLineStyle(4); | |
889 | Leg1->AddEntry(secMCPred,"Sec.Mat","l"); | |
88d02ae4 | 890 | secMCPred->DrawClone("lhistsame"); |
2245e1bd | 891 | |
892 | } | |
893 | } | |
88d02ae4 | 894 | else |
895 | { | |
896 | hconWD[index+6*isample]->SetBinContent(ibin_data,0.0); | |
897 | hconWD[index+6*isample]->SetBinError(ibin_data,0.0); | |
898 | hconMat[index+6*isample]->SetBinContent(ibin_data,0.0); | |
899 | hconMat[index+6*isample]->SetBinError(ibin_data,0.0); | |
900 | hcon[index+6*isample]->SetBinContent(ibin_data,0.0); | |
901 | hcon[index+6*isample]->SetBinError(ibin_data,0.0); | |
902 | hprimary[index+6*isample]->SetBinContent(ibin_data,1.0); | |
903 | hprimary[index+6*isample]->SetBinError(ibin_data,0.0); | |
904 | } | |
b6a74f8a | 905 | Leg1->Draw(); |
88d02ae4 | 906 | listofdcafits->Add(cDCA); |
907 | ||
908 | //cDCA->Write(); | |
909 | delete hToFit; | |
910 | } | |
911 | ||
912 | ||
913 | } | |
914 | ||
915 | } | |
916 | } | |
917 | fout->cd(); | |
918 | listofdcafits->Write("DCAfits",TObject::kSingleKey); | |
919 | } | |
920 | ||
921 | void RecomputeErrors(TH1* h) | |
922 | { | |
923 | for (int i=0; i<=h->GetXaxis()->GetNbins(); i++) | |
50c532c4 | 924 | { |
925 | cout<<h->GetBinContent(i)<<" "<<h->GetBinError(i)<<" error "<<TMath::Sqrt(h->GetBinContent(i))<<endl; | |
88d02ae4 | 926 | h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i))); |
50c532c4 | 927 | } |
88d02ae4 | 928 | h->Sumw2(); |
929 | } | |
930 | void SetBintoOne(TH1* h) | |
931 | { | |
932 | for (int i=0;i<=h->GetXaxis()->GetNbins(); i++) | |
933 | { | |
934 | h->SetBinContent(i,1); | |
935 | h->SetBinError(i,0); | |
936 | } | |
937 | } | |
938 | ||
939 | ||
940 | void GetCorrectedSpectra(TH1F* corr,TH1F* raw,TH1F* eff, TH1F* con) | |
941 | { | |
942 | for (int i=0;i<=corr->GetXaxis()->GetNbins(); i++) | |
943 | { | |
944 | corr->SetBinContent(i,1); | |
945 | corr->SetBinError(i,0); | |
946 | } | |
947 | corr->Sumw2(); | |
948 | corr->Add(con,-1); | |
949 | corr->Sumw2(); | |
950 | corr->Divide(eff); | |
951 | corr->Sumw2(); | |
952 | corr->Multiply(raw); | |
953 | corr->Sumw2(); | |
954 | } | |
955 | void GetCorrectedSpectraLeonardo(TH1F* spectra,TH1F* correction, TH1F* hprimaryData,TH1F* hprimaryMC) | |
956 | { | |
957 | spectra->Sumw2(); | |
958 | spectra->Multiply(correction); | |
959 | spectra->Sumw2(); | |
960 | hprimaryData->Sumw2(); | |
961 | spectra->Multiply(hprimaryData); | |
962 | hprimaryMC->Sumw2(); | |
963 | spectra->Divide(hprimaryMC); | |
964 | } | |
965 | ||
b3ea73e1 | 966 | void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options) |
88d02ae4 | 967 | { |
b3ea73e1 | 968 | if (options&kgeantflukaKaon) |
969 | { | |
970 | //Geant/Fluka Correction | |
971 | Printf("\nGF correction for Kaons"); | |
972 | //Getting GF For Kaons in TPC | |
973 | TGraph *gGFCorrectionKaonPlus=new TGraph(); | |
974 | gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus"); | |
975 | gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus"); | |
976 | TGraph *gGFCorrectionKaonMinus=new TGraph(); | |
977 | gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus"); | |
978 | gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus"); | |
979 | TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root"; | |
980 | TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data()); | |
2245e1bd | 981 | if (!fGeanFlukaK) |
982 | { | |
ab0298c1 | 983 | fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.321.root"; |
2245e1bd | 984 | fGeanFlukaK= new TFile(fnameGeanFlukaK.Data()); |
985 | if (!fGeanFlukaK) | |
986 | return; | |
987 | } | |
b3ea73e1 | 988 | TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles"); |
989 | TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles"); | |
990 | //getting GF func for Kaons with TOF | |
991 | TF1 *fGFKPosTracking; | |
992 | fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive); | |
993 | TF1 *fGFKNegTracking; | |
994 | fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative); | |
995 | TF1 *fGFKPosMatching; | |
996 | fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive); | |
997 | TF1 *fGFKNegMatching; | |
998 | fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative); | |
999 | for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++) | |
1000 | { | |
88d02ae4 | 1001 | if(Spectra[1]->GetBinCenter(binK)<tofpt) |
1002 | {//use TPC GeantFlukaCorrection | |
1003 | Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK))); | |
1004 | Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK))); | |
1005 | // Printf("TPC Geant/Fluka: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg); | |
1006 | Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos); | |
1007 | Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg); | |
1008 | Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos); | |
1009 | Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg); | |
1010 | gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos); | |
1011 | gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg); | |
1012 | } | |
1013 | else | |
1014 | { | |
1015 | gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0); | |
1016 | gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0); | |
1017 | Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK)); | |
1018 | Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK)); | |
1019 | // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking); | |
1020 | Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking); | |
1021 | Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking); | |
1022 | Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking); | |
1023 | Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking); | |
1024 | Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK)); | |
1025 | Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK)); | |
1026 | // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching); | |
1027 | Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching); | |
1028 | Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching); | |
1029 | Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching); | |
1030 | Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching); | |
1031 | } | |
b3ea73e1 | 1032 | } |
88d02ae4 | 1033 | } |
b3ea73e1 | 1034 | if(!(options&kgeantflukaProton)) |
1035 | return; | |
88d02ae4 | 1036 | //Geant Fluka for P in TPC |
1037 | Printf("\nGF correction for Protons"); | |
1038 | const Int_t kNCharge=2; | |
1039 | Int_t kPos=0; | |
1040 | Int_t kNeg=1; | |
50c532c4 | 1041 | TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root"; |
2245e1bd | 1042 | TFile* fGFProtons = new TFile (fnameGFProtons.Data()); |
1043 | if (!fGFProtons) | |
1044 | { | |
ab0298c1 | 1045 | fnameGFProtons=="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.root"; |
2245e1bd | 1046 | TFile* fGFProtons = new TFile (fnameGFProtons.Data()); |
1047 | if (!fGFProtons) | |
1048 | return; | |
1049 | } | |
1050 | ||
1051 | ||
1052 | ||
88d02ae4 | 1053 | TH2D * hCorrFluka[kNCharge]; |
1054 | TH2D * hCorrFluka[2]; | |
1055 | hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons"); | |
1056 | hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons"); | |
1057 | //getting GF func for Kaons with TPCTOF | |
1058 | TF1 *fGFpPosTracking; | |
1059 | fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive); | |
1060 | TF1 *fGFpNegTracking; | |
1061 | fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative); | |
1062 | TF1 *fGFpPosMatching; | |
1063 | fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive); | |
1064 | TF1 *fGFpNegMatching; | |
1065 | fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative); | |
1066 | ||
1067 | ||
1068 | Int_t nbins = Spectra[2]->GetNbinsX(); | |
1069 | ||
1070 | for(Int_t ibin = 0; ibin < nbins; ibin++) | |
1071 | { | |
1072 | if(Spectra[2]->GetBinCenter(ibin)<tofpt) | |
1073 | {//use TPC GeantFlukaCorrection | |
1074 | for(Int_t icharge = 0; icharge < kNCharge; icharge++) | |
1075 | { | |
1076 | Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY(); | |
1077 | Float_t pt = Spectra[2]->GetBinCenter(ibin); | |
1078 | Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1); | |
1079 | Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1); | |
1080 | if (pt < minPtCorrection) | |
1081 | pt = minPtCorrection+0.0001; | |
1082 | if (pt > maxPtCorrection) | |
1083 | pt = maxPtCorrection; | |
1084 | Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt)); | |
1085 | ||
1086 | if (correction > 0.0) | |
1087 | {// If the bin is empty this is a 0 | |
1088 | // cout<<icharge<<" "<<ibin<<" "<<correction<<endl; | |
1089 | Spectra[icharge*3+2]->SetBinContent(ibin,Spectra[icharge*3+2]->GetBinContent(ibin)*correction); | |
1090 | Spectra[icharge*3+2]->SetBinError(ibin,Spectra[icharge*3+2]->GetBinError(ibin)*correction); | |
1091 | } | |
1092 | else if (Spectra[icharge*3+2]->GetBinContent(ibin) > 0.0) | |
1093 | { | |
1094 | // If we are skipping a non-empty bin, we notify the user | |
1095 | cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for " <<"protons Pt:"<< pt<< endl; | |
1096 | cout << " Bin content: " << Spectra[icharge*3+2]->GetBinContent(ibin) << endl; | |
1097 | } | |
1098 | } | |
1099 | } | |
1100 | else | |
1101 | { | |
1102 | Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin)); | |
1103 | Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin)); | |
1104 | // Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking); | |
1105 | Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking); | |
1106 | Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking); | |
1107 | Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking); | |
1108 | Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking); | |
1109 | Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin)); | |
1110 | Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin)); | |
1111 | // Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching); | |
1112 | Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching); | |
1113 | Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching); | |
1114 | Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching); | |
1115 | Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching); | |
1116 | } | |
1117 | } | |
1118 | fGeanFlukaK->Close(); | |
1119 | delete fGeanFlukaK; | |
1120 | } | |
1121 | ||
1122 | ||
1123 | /////////// | |
1124 | TF1 * | |
1125 | TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge) | |
1126 | { | |
1127 | ||
1128 | if (ipart == 3 && icharge == kNegative) { | |
1129 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.); | |
1130 | return f; | |
1131 | } | |
1132 | else if (ipart == 4 && icharge == kNegative) { | |
1133 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.); | |
1134 | } | |
1135 | else | |
1136 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.); | |
1137 | ||
1138 | return f; | |
1139 | } | |
1140 | ||
1141 | Double_t | |
1142 | TrackingPtGeantFlukaCorrectionNull(Double_t pTmc) | |
1143 | { | |
1144 | return 1.; | |
1145 | } | |
1146 | ||
1147 | Double_t | |
1148 | TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc) | |
1149 | { | |
1150 | return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612)); | |
1151 | } | |
1152 | ||
1153 | Double_t | |
1154 | TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc) | |
1155 | { | |
1156 | return TMath::Min((0.972865 + 0.0117093*pTmc), 1.); | |
1157 | } | |
1158 | /////////////////////////////////////////// | |
1159 | TF1 * | |
1160 | TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge) | |
1161 | { | |
1162 | ||
1163 | if (ipart == 3 && icharge == kNegative) { | |
1164 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.); | |
1165 | return f; | |
1166 | } | |
1167 | else if (ipart == 4 && icharge == kNegative) { | |
1168 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.); | |
1169 | } | |
1170 | else | |
1171 | TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.); | |
1172 | ||
1173 | return f; | |
1174 | } | |
1175 | ||
1176 | ||
1177 | Double_t | |
1178 | MatchingPtGeantFlukaCorrectionNull(Double_t pTmc) | |
1179 | { | |
1180 | return 1.; | |
1181 | } | |
1182 | ||
1183 | Double_t | |
1184 | MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc) | |
1185 | { | |
1186 | Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094)); | |
1187 | return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471)); | |
1188 | } | |
1189 | ||
1190 | Double_t | |
1191 | MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc) | |
1192 | { | |
1193 | Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc); | |
1194 | return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.); | |
1195 | } | |
1196 | void MatchingTOFEff(TH1F** Spectra, TList* list=0x0) | |
1197 | { | |
1198 | if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0) | |
1199 | { | |
50c532c4 | 1200 | TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos(); |
1201 | hMatcEffPos_data->Sumw2(); | |
1202 | //hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos()); | |
1203 | hMatcEffPos_data->Divide(hMatcEffPos_data,(TH1F*)tcutsdata->GetHistoNSelectedPos(),1,1,"B"); | |
88d02ae4 | 1204 | hMatcEffPos_data->SetTitle("Matching Eff Pos - data"); |
1205 | TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg(); | |
50c532c4 | 1206 | hMatcEffNeg_data->Sumw2(); |
1207 | //hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg()); | |
1208 | hMatcEffNeg_data->Divide(hMatcEffNeg_data,(TH1F*)tcutsdata->GetHistoNSelectedNeg(),1,1,"B"); | |
88d02ae4 | 1209 | hMatcEffNeg_data->SetTitle("Matching Eff Neg - data"); |
1210 | TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos(); | |
50c532c4 | 1211 | hMatcEffPos_mc->Sumw2(); |
1212 | //hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos()); | |
1213 | hMatcEffPos_mc->Divide(hMatcEffPos_mc,(TH1F*)tcutsmc->GetHistoNSelectedPos(),1,1,"B"); | |
88d02ae4 | 1214 | hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc"); |
1215 | TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg(); | |
50c532c4 | 1216 | hMatcEffNeg_mc->Sumw2(); |
1217 | //hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg()); | |
1218 | hMatcEffNeg_mc->Divide(hMatcEffNeg_mc,(TH1F*)tcutsmc->GetHistoNSelectedNeg(),1,1,"B"); | |
88d02ae4 | 1219 | hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc"); |
1220 | ||
1221 | ||
1222 | hMatcEffPos_data->Divide(hMatcEffPos_mc); | |
1223 | hMatcEffNeg_data->Divide(hMatcEffNeg_mc); | |
1224 | hMatcEffPos_data->SetName("MatchingTOFPos"); | |
1225 | hMatcEffNeg_data->SetName("MatchingTOFNeg"); | |
1226 | ||
1227 | ||
1228 | TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",.6,5); | |
1229 | hMatcEffPos_data->Fit("pol0MatchPos_data","MNR"); | |
1230 | TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",.6,5); | |
1231 | hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR"); | |
1232 | ||
1233 | list->Add(hMatcEffPos_data); | |
1234 | list->Add(hMatcEffNeg_data); | |
1235 | ||
1236 | ||
1237 | TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0); | |
1238 | TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0); | |
50c532c4 | 1239 | |
88d02ae4 | 1240 | } |
1241 | //Correction spectra for matching efficiency | |
1242 | //For the moment I'm using the inclusive correction | |
1243 | for(Int_t ipart=0;ipart<3;ipart++) | |
1244 | { | |
1245 | ||
1246 | for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++) | |
1247 | { | |
1248 | Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin); | |
1249 | if(ptspectra<tcutsdata->GetPtTOFMatching()) | |
1250 | continue; | |
1251 | //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra)))); | |
1252 | //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra)))); | |
1253 | Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/TOFMatchingScalling[0])); | |
1254 | Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/TOFMatchingScalling[1])); | |
1255 | } | |
1256 | } | |
1257 | } | |
1258 | Double_t eta2y(Double_t pt, Double_t mass, Double_t eta) | |
1259 | { | |
1260 | Double_t mt = TMath::Sqrt(mass * mass + pt * pt); | |
1261 | return TMath::ASinH(pt / mt * TMath::SinH(eta)); | |
1262 | } | |
1263 | ||
ab0298c1 | 1264 | TH1* GetSumAllCh(TH1F** spectra, Double_t* mass,Double_t etacut) |
88d02ae4 | 1265 | { |
1266 | TH1F* allch=(((TH1F*))spectra[0]->Clone("allCh")); | |
1267 | allch->Reset(); | |
1268 | for (int i=0;i<6;i++) | |
1269 | { | |
1270 | Double_t masstmp=mass[i%3]; | |
1271 | for (int j=1;j<spectra[i]->GetXaxis()->GetNbins();j++) | |
1272 | { | |
1273 | Float_t value=spectra[i]->GetBinContent(j); | |
1274 | Float_t error=spectra[i]->GetBinError(j); | |
1275 | if(value>0.0) | |
1276 | { | |
1277 | Float_t pt=spectra[i]->GetXaxis()->GetBinCenter(j); | |
1278 | Float_t mt2=masstmp*masstmp+pt*pt; | |
1279 | Float_t mt=TMath::Sqrt(mt2); | |
ab0298c1 | 1280 | Float_t maxy=eta2y(pt,masstmp,etacut/2.0); |
88d02ae4 | 1281 | Float_t conver=maxy*(TMath::Sqrt(1-masstmp*masstmp/(mt2*TMath::CosH(maxy)*TMath::CosH(maxy)))+TMath::Sqrt(1-masstmp*masstmp/(mt2*TMath::CosH(0.0)*TMath::CosH(0.0)))); |
ab0298c1 | 1282 | conver=conver/etacut; |
a4ca64e3 | 1283 | //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl; |
88d02ae4 | 1284 | Float_t bincontent=allch->GetBinContent(j); |
1285 | Float_t binerror=allch->GetBinError(j); | |
1286 | bincontent+=conver*value; | |
1287 | binerror=TMath::Sqrt(binerror*binerror+conver*conver*error*error); | |
1288 | allch->SetBinContent(j,bincontent); | |
1289 | allch->SetBinError(j,binerror); | |
1290 | } | |
1291 | ||
1292 | } | |
1293 | } | |
1294 | Divideby2pipt(allch); | |
1295 | allch->SetTitle("N_{ch};p_{T} (GeV/c);1/(2 #pi p_{T})dN/p_{T} |#eta|<0.8"); | |
1296 | return allch; | |
1297 | } | |
1298 | ||
1299 | void Divideby2pipt(TH1* hist) | |
1300 | { | |
1301 | ||
1302 | for (int i=1;i<hist->GetXaxis()->GetNbins();i++) | |
1303 | { | |
1304 | Float_t value=hist->GetBinContent(i); | |
1305 | Float_t error=hist->GetBinError(i); | |
1306 | Float_t pt=hist->GetXaxis()->GetBinCenter(i); | |
1307 | hist->SetBinContent(i,value/(2.0*TMath::Pi()*pt)); | |
1308 | hist->SetBinError(i,error/(2.0*TMath::Pi()*pt)); | |
1309 | ||
1310 | } | |
1311 | } | |
2245e1bd | 1312 | |
1313 | Short_t DCAfitsettings (Float_t pt, Int_t type) | |
1314 | { | |
2245e1bd | 1315 | Short_t value=0x0; |
1316 | if (pt<maxptformaterial[type]&&pt>minptformaterial[type]) | |
1317 | value=value+2; | |
1318 | if (pt<maxptforWD[type]&&pt>minptforWD[type]) | |
1319 | value=value+1; | |
1320 | return value; | |
1321 | ||
3a71f081 | 1322 | } |
1323 | ||
1324 | Float_t Normaliztionwithbin0integrals(UInt_t options) | |
1325 | { | |
1326 | ||
1327 | TH1F* bin0mcRec=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_rec"); | |
1328 | TH1F* bin0mcMC=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_MC"); | |
1329 | ||
1330 | TH1F* vertexmc=ecutsmc->GetHistoVtxAftSelwithoutZvertexCut(); | |
1331 | TH1F* vertexmcMCz=ecutsmc->GetHistoVtxAftSelwithoutZvertexCutusingMCz(); | |
1332 | TH1F* vertexdata=ecutsdata->GetHistoVtxAftSelwithoutZvertexCut(); | |
1333 | ||
1334 | TH1I* histodata=ecutsdata->GetHistoCuts(); | |
1335 | TH1I* histomc=ecutsmc->GetHistoCuts(); | |
1336 | ||
1337 | Float_t dataevents=(Float_t)histodata->GetBinContent(3); | |
a4ca64e3 | 1338 | //cout<<histodata->GetBinContent(2)<<endl; |
3a71f081 | 1339 | Float_t databin0events=((Float_t)histodata->GetBinContent(2))-((Float_t)histodata->GetBinContent(4)); |
1340 | ||
1341 | bin0mcRec->Sumw2(); | |
1342 | bin0mcMC->Sumw2(); | |
1343 | ||
1344 | bin0mcRec->Add(vertexmc,-1); | |
1345 | bin0mcMC->Add(vertexmcMCz,-1); | |
1346 | ||
1347 | bin0mcRec->Divide(vertexmc); | |
1348 | bin0mcMC->Divide(vertexmcMCz); | |
1349 | ||
1350 | bin0mcRec->Multiply(vertexdata); | |
1351 | bin0mcMC->Multiply(vertexdata); | |
1352 | ||
1353 | Float_t bin0mcRecN=0.0; | |
1354 | Float_t bin0mcMCN=0.0; | |
1355 | ||
1356 | for (int i=0;i<=bin0mcRec->GetXaxis()->GetNbins();i++) | |
1357 | { | |
1358 | bin0mcRecN+=bin0mcRec->GetBinContent(i); | |
1359 | bin0mcMCN+=bin0mcMC->GetBinContent(i); | |
1360 | ||
1361 | } | |
1362 | bin0mcRec->Scale(databin0events/bin0mcRecN); | |
1363 | bin0mcMC->Scale(databin0events/bin0mcMCN); | |
1364 | ||
1365 | Int_t binmin=bin0mcRec->GetXaxis()->FindBin(-10); | |
1366 | Int_t binmax=bin0mcRec->GetXaxis()->FindBin(10)-1; | |
1367 | cout<< bin0mcRec->GetXaxis()->GetBinLowEdge(binmin)<<" "<<bin0mcRec->GetXaxis()->GetBinUpEdge(binmax)<<endl; | |
1368 | cout<<bin0mcRecN<<" "<<bin0mcMCN<<" "<<databin0events<<endl; | |
1369 | cout<<dataevents<<" normalization "<<dataevents+bin0mcRec->Integral(binmin,binmax)<<" "<<dataevents+bin0mcMC->Integral(binmin,binmax)<<endl; | |
1370 | cout<<histodata->GetBinContent(2)<<" "<<histodata->GetBinContent(4)<<endl; | |
1371 | if ((options&knormalizationwithbin0integralsdata)==knormalizationwithbin0integralsdata) | |
1372 | return dataevents+bin0mcRec->Integral(binmin,binmax); | |
1af4558e | 1373 | else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsMC) |
3a71f081 | 1374 | return dataevents+bin0mcMC->Integral(binmin,binmax) ; |
1375 | else | |
1376 | return 1; | |
1377 | } | |
1378 | ||
a4ca64e3 | 1379 | |
1380 | Bool_t ReadConfigFile(TString configfile) | |
1381 | { | |
1382 | ifstream infile(configfile.Data()); | |
1383 | if(infile.is_open()==false) | |
1384 | return false; | |
cd28560b | 1385 | TString namesofSetting[38]={"CutRangeMin","CutRangeMax","FitRangeMin","FitRangeMax","MinMatPionPlus","MaxMatPionPlus","MinMatKaonPlus","MaxMatKaonPlus","MinMatProtonPlus","MaxMatProtonPlus","MinMatPionMinus","MaxMatPionMinus","MinMatKaonMinus","MaxMatKaonMinus","MinMatProtonMinus","MaxMatProtonMinus","MinWDPionPlus","MaxWDPionPlus","MinWDKaonPlus","MaxWDKaonPlus","MinWDProtonPlus","MaxWDProtonPlus","MinWDPionMinus","MaxWDPionMinus","MinWDKaonMinus","MaxWDKaonMinus","MinWDProtonMinus","MaxWDProtonMinus","MaxContaminationPIDMC","MinPions","MaxPions","MinKaons","MaxKaons","MinProtons","MaxProtons","TOFPIDsignalmatchPion","TOFPIDsignalmatchKaon","TOFPIDsignalmatchProton"}; |
a4ca64e3 | 1386 | |
1387 | char buffer[256]; | |
1388 | while (infile.eof()==false) | |
1389 | { | |
1390 | buffer[0]='#'; | |
1391 | while (buffer[0]=='#'&&infile.eof()==false) | |
1392 | infile.getline(buffer,256); | |
1393 | TString tmpstring(buffer); | |
1394 | cout<<buffer<<endl; | |
1395 | if(tmpstring.Contains(namesofSetting[0])) | |
1396 | CutRange[0]=(tmpstring.Remove(0,namesofSetting[0].Length()+1)).Atof(); | |
1397 | else if (tmpstring.Contains(namesofSetting[1])) | |
1398 | CutRange[1]=(tmpstring.Remove(0,namesofSetting[1].Length()+1)).Atof(); | |
1399 | else if (tmpstring.Contains(namesofSetting[2])) | |
1400 | FitRange[0]=(tmpstring.Remove(0,namesofSetting[2].Length()+1)).Atof(); | |
1401 | else if (tmpstring.Contains(namesofSetting[3])) | |
1402 | FitRange[1]=(tmpstring.Remove(0,namesofSetting[3].Length()+1)).Atof(); | |
1403 | else if (tmpstring.Contains(namesofSetting[4])) | |
1404 | minptformaterial[0]=(tmpstring.Remove(0,namesofSetting[4].Length()+1)).Atof(); | |
1405 | else if (tmpstring.Contains(namesofSetting[5])) | |
1406 | maxptformaterial[0]=(tmpstring.Remove(0,namesofSetting[5].Length()+1)).Atof(); | |
1407 | else if (tmpstring.Contains(namesofSetting[6])) | |
1408 | minptformaterial[1]=(tmpstring.Remove(0,namesofSetting[6].Length()+1)).Atof(); | |
1409 | else if (tmpstring.Contains(namesofSetting[7])) | |
1410 | maxptformaterial[1]=(tmpstring.Remove(0,namesofSetting[7].Length()+1)).Atof(); | |
1411 | else if (tmpstring.Contains(namesofSetting[8])) | |
1412 | minptformaterial[2]=(tmpstring.Remove(0,namesofSetting[8].Length()+1)).Atof(); | |
1413 | else if (tmpstring.Contains(namesofSetting[9])) | |
1414 | maxptformaterial[2]=(tmpstring.Remove(0,namesofSetting[9].Length()+1)).Atof(); | |
1415 | else if (tmpstring.Contains(namesofSetting[10])) | |
1416 | minptformaterial[3]=(tmpstring.Remove(0,namesofSetting[10].Length()+1)).Atof(); | |
1417 | else if (tmpstring.Contains(namesofSetting[11])) | |
1418 | maxptformaterial[3]=(tmpstring.Remove(0,namesofSetting[11].Length()+1)).Atof(); | |
1419 | else if (tmpstring.Contains(namesofSetting[12])) | |
1420 | minptformaterial[4]=(tmpstring.Remove(0,namesofSetting[12].Length()+1)).Atof(); | |
1421 | else if (tmpstring.Contains(namesofSetting[13])) | |
1422 | maxptformaterial[4]=(tmpstring.Remove(0,namesofSetting[13].Length()+1)).Atof(); | |
1423 | else if (tmpstring.Contains(namesofSetting[14])) | |
1424 | minptformaterial[5]=(tmpstring.Remove(0,namesofSetting[14].Length()+1)).Atof(); | |
1425 | else if (tmpstring.Contains(namesofSetting[15])) | |
1426 | maxptformaterial[5]=(tmpstring.Remove(0,namesofSetting[15].Length()+1)).Atof(); | |
1427 | else if (tmpstring.Contains(namesofSetting[16])) | |
1428 | minptforWD[0]=(tmpstring.Remove(0,namesofSetting[16].Length()+1)).Atof(); | |
1429 | else if (tmpstring.Contains(namesofSetting[17])) | |
1430 | maxptforWD[0]=(tmpstring.Remove(0,namesofSetting[17].Length()+1)).Atof(); | |
1431 | else if (tmpstring.Contains(namesofSetting[18])) | |
1432 | minptforWD[1]=(tmpstring.Remove(0,namesofSetting[18].Length()+1)).Atof(); | |
1433 | else if (tmpstring.Contains(namesofSetting[19])) | |
1434 | maxptforWD[1]=(tmpstring.Remove(0,namesofSetting[19].Length()+1)).Atof(); | |
1435 | else if (tmpstring.Contains(namesofSetting[20])) | |
1436 | minptforWD[2]=(tmpstring.Remove(0,namesofSetting[20].Length()+1)).Atof(); | |
1437 | else if (tmpstring.Contains(namesofSetting[21])) | |
1438 | maxptforWD[2]=(tmpstring.Remove(0,namesofSetting[21].Length()+1)).Atof(); | |
1439 | else if (tmpstring.Contains(namesofSetting[22])) | |
1440 | minptforWD[3]=(tmpstring.Remove(0,namesofSetting[22].Length()+1)).Atof(); | |
1441 | else if (tmpstring.Contains(namesofSetting[23])) | |
1442 | maxptforWD[3]=(tmpstring.Remove(0,namesofSetting[23].Length()+1)).Atof(); | |
1443 | else if (tmpstring.Contains(namesofSetting[24])) | |
1444 | minptforWD[4]=(tmpstring.Remove(0,namesofSetting[24].Length()+1)).Atof(); | |
1445 | else if (tmpstring.Contains(namesofSetting[25])) | |
1446 | maxptforWD[4]=(tmpstring.Remove(0,namesofSetting[25].Length()+1)).Atof(); | |
1447 | else if (tmpstring.Contains(namesofSetting[26])) | |
1448 | minptforWD[5]=(tmpstring.Remove(0,namesofSetting[26].Length()+1)).Atof(); | |
1449 | else if (tmpstring.Contains(namesofSetting[27])) | |
1450 | maxptforWD[5]=(tmpstring.Remove(0,namesofSetting[27].Length()+1)).Atof(); | |
1af4558e | 1451 | else if (tmpstring.Contains(namesofSetting[28])) |
1452 | fMaxContaminationPIDMC=(tmpstring.Remove(0,namesofSetting[28].Length()+1)).Atof(); | |
1453 | else if (tmpstring.Contains(namesofSetting[29])) | |
1454 | minRanges[0]=(tmpstring.Remove(0,namesofSetting[29].Length()+1)).Atof(); | |
1455 | else if (tmpstring.Contains(namesofSetting[30])) | |
1456 | maxRanges[0]=(tmpstring.Remove(0,namesofSetting[30].Length()+1)).Atof(); | |
1457 | else if (tmpstring.Contains(namesofSetting[31])) | |
1458 | minRanges[1]=(tmpstring.Remove(0,namesofSetting[31].Length()+1)).Atof(); | |
1459 | else if (tmpstring.Contains(namesofSetting[32])) | |
1460 | maxRanges[1]=(tmpstring.Remove(0,namesofSetting[32].Length()+1)).Atof(); | |
1461 | else if (tmpstring.Contains(namesofSetting[33])) | |
1462 | minRanges[2]=(tmpstring.Remove(0,namesofSetting[33].Length()+1)).Atof(); | |
1463 | else if (tmpstring.Contains(namesofSetting[34])) | |
1464 | maxRanges[2]=(tmpstring.Remove(0,namesofSetting[34].Length()+1)).Atof(); | |
cd28560b | 1465 | else if (tmpstring.Contains(namesofSetting[35])) |
1466 | TOFPIDsignalmatching[0]=(tmpstring.Remove(0,namesofSetting[35].Length()+1)).Atof(); | |
1467 | else if (tmpstring.Contains(namesofSetting[36])) | |
1468 | TOFPIDsignalmatching[1]=(tmpstring.Remove(0,namesofSetting[36].Length()+1)).Atof(); | |
1469 | else if (tmpstring.Contains(namesofSetting[37])) | |
1470 | TOFPIDsignalmatching[2]=(tmpstring.Remove(0,namesofSetting[37].Length()+1)).Atof(); | |
1471 | else | |
a4ca64e3 | 1472 | continue; |
1473 | ||
1af4558e | 1474 | |
1475 | // Double_t minRanges[3]={0.3,0.3,0.45}; | |
1476 | //Double_t maxRanges[3]={1.5,1.2,2.2}; | |
1477 | //Double_t fMaxContaminationPIDMC=0.2; | |
1478 | ||
1479 | ||
1480 | ||
a4ca64e3 | 1481 | } |
1482 | for(int i=0;i<6;i++) | |
1483 | cout<<minptformaterial[i]<<" "<<maxptformaterial[i]<<" "<<minptforWD[i]<<" "<<maxptforWD[i]<<endl; | |
1484 | cout<<FitRange[0]<<" "<<FitRange[1]<<" "<<CutRange[0]<<CutRange[1]<<endl; | |
1485 | if(FitRange[0]>=FitRange[1]) | |
1486 | { | |
1487 | cout<<"A"<<endl; | |
1488 | return false; | |
1489 | } | |
1490 | if(CutRange[0]>=CutRange[1]) | |
1491 | { | |
1492 | cout<<"B"<<endl; | |
1493 | return false; | |
1494 | } | |
1495 | for(int i=0;i<6;i++) | |
1496 | { | |
1497 | if((minptformaterial[i]>maxptformaterial[i]&&minptformaterial[i]>0.0)||minptformaterial[i]<0.0||maxptformaterial[i]<0.0) | |
1498 | { | |
1499 | cout<<"C"<<endl; | |
1500 | return false; | |
1501 | } | |
1502 | if((minptforWD[i]>maxptforWD[i]&&minptforWD[i]>0.0)||minptforWD[i]<0.0||maxptforWD[i]<0.0) | |
1503 | { | |
1504 | cout<<"D"<<endl; | |
1505 | return false; | |
1506 | } | |
1507 | } | |
1af4558e | 1508 | for(int i=0;i<3;i++) |
1509 | if(minRanges[i]>maxRanges[i]) | |
1510 | return false; | |
a4ca64e3 | 1511 | |
1512 | ||
1513 | return true; | |
1514 | } | |
50c532c4 | 1515 | |
1516 | void SubHistWithFullCorr(TH1F* h1, TH1F* h2, Float_t factor1=1.0, Float_t factor2=1.0) | |
1517 | { | |
1518 | if(h1->GetNbinsX()!=h2->GetNbinsX()) | |
1519 | return; | |
1520 | for (int i=0;i<=h1->GetNbinsX();i++) | |
1521 | { | |
1522 | Float_t tmpvalue=factor1*h1->GetBinContent(i)-factor2*h2->GetBinContent(i); | |
1523 | Float_t tmperror=TMath::Abs(factor1*factor1*h1->GetBinError(i)*h1->GetBinError(i)-factor2*factor2*h2->GetBinError(i)*h2->GetBinError(i)); | |
1524 | h1->SetBinContent(i,tmpvalue); | |
1525 | h1->SetBinError(i,TMath::Sqrt(tmperror)); | |
1526 | } | |
1527 | ||
1af4558e | 1528 | } |
1529 | ||
1530 | void TOFMatchingForNch(TH1* h) | |
1531 | { | |
1532 | if(TOFMatchingScalling[0]>0.0&&TOFMatchingScalling[1]>0.0) | |
1533 | { | |
1534 | Float_t factor=0.5*TOFMatchingScalling[0]+0.5*TOFMatchingScalling[1]; | |
1535 | for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++) | |
1536 | { | |
1537 | Float_t ptspectra=h->GetBinCenter(ibin); | |
1538 | if(ptspectra<tcutsdata->GetPtTOFMatching()) | |
1539 | continue; | |
1540 | h->SetBinContent(ibin,(h->GetBinContent(ibin)/factor)); | |
1541 | } | |
1542 | ||
1543 | } | |
1544 | else | |
1545 | return; | |
1546 | ||
1547 | ||
cd28560b | 1548 | } |
1549 | void TOFPIDsignalmatchingApply(TH1* h, Float_t factor) | |
1550 | { | |
1551 | if(factor<0.0) | |
1552 | return; | |
1553 | for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++) | |
1554 | { | |
1555 | Float_t ptspectra=h->GetBinCenter(ibin); | |
1556 | if(ptspectra<tcutsdata->GetPtTOFMatching()) | |
1557 | continue; | |
1558 | h->SetBinContent(ibin,(h->GetBinContent(ibin)*factor)); | |
1559 | ||
1560 | } | |
1561 | ||
1562 | } |