]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/macros/DrawQAoutput.C
Transition PWG3 --> PWGHF
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / DrawQAoutput.C
CommitLineData
9af24f46 1#include <Riostream.h>
2#include <TFile.h>
3#include <TString.h>
4#include <TH2F.h>
5#include <TH1F.h>
6#include <TF1.h>
1fc8d3c0 7#include <TGraph.h>
9af24f46 8#include <TDirectoryFile.h>
9#include <TList.h>
10#include <TCanvas.h>
11#include <TLegend.h>
12#include <TPaveText.h>
13#include <TStyle.h>
82f89d0d 14#include <TClass.h>
1fc8d3c0 15#include <TDatabasePDG.h>
16#include <TParameter.h>
82f89d0d 17#include <AliCounterCollection.h>
1fc8d3c0 18#include <AliRDHFCuts.h>
9af24f46 19
04fd2703 20TString *periodsname;
21
9af24f46 22//read the file and take list and stat
23
04fd2703 24Bool_t ReadFile(TList* &list,TH1F* &hstat, TString listname,TString partname,TString path="./",TString filename=/*"AnalysisResults.root"*/"PWG3histograms.root");
25Bool_t ReadFileMore(TList* &list,TH1F* &hstat, AliRDHFCuts* &cutobj, TString listname,TString partname,TString path="./",TString filename=/*"AnalysisResults.root"*/"PWG3histograms.root");
26void SuperimposeBBToTPCSignal(Int_t period /*0=LHC10bc, 1=LHC10d, 2=LHC10h*/,TCanvas* cpid, Int_t set);
1fc8d3c0 27void TPCBetheBloch(Int_t set);
9af24f46 28
1fc8d3c0 29Bool_t ReadFile(TList* &list,TH1F* &hstat, TString listname,TString partname,TString path,TString filename){
9af24f46 30
1fc8d3c0 31 TString hstatname="nEntriesQA",dirname="PWG3_D2H_QA", cutobjname="";
9af24f46 32 filename.Prepend(path);
33 listname+=partname;
34 hstatname+=partname;
35
36 TFile* f=new TFile(filename.Data());
37 if(!f){
38 cout<<filename.Data()<<" not found"<<endl;
39 return kFALSE;
40 }
41 TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname);
42 if(!f){
43 cout<<dirname.Data()<<" not found in "<<filename.Data()<<endl;
44 return kFALSE;
45 }
46
47 list=(TList*)dir->Get(listname);
48 if(!list){
49 cout<<"List "<<listname.Data()<<" not found"<<endl;
50 dir->ls();
51 return kFALSE;
52 }
53
54 hstat=(TH1F*)dir->Get(hstatname);
55 if(!hstat){
56 cout<<hstatname.Data()<<" not found"<<endl;
57 return kFALSE;
58 }
1fc8d3c0 59
60 return kTRUE;
61}
62
63Bool_t ReadFileMore(TList* &list,TH1F* &hstat, AliRDHFCuts* &cutobj, TString listname,TString partname,TString path,TString filename){
64
65 TString hstatname="nEntriesQA",dirname="PWG3_D2H_QA", cutobjname="";
66 filename.Prepend(path);
67 listname+=partname;
68 hstatname+=partname;
69
f3e1ad2f 70 if(partname.Contains("Dplus")) cutobjname="AnalysisCuts";//"DplustoKpipiCutsStandard";
1fc8d3c0 71 else{
04fd2703 72 if(partname.Contains("D0")) cutobjname="D0toKpiCuts";//"D0toKpiCutsStandard"
1fc8d3c0 73 else{
74 if(partname.Contains("Dstar")) cutobjname="DStartoKpipiCuts";
75 else{
76 if(partname.Contains("Ds")) cutobjname="DstoKKpiCuts";
77 else{
78 if(partname.Contains("D04")) cutobjname="D0toKpipipiCuts";
79 else{
80 if(partname.Contains("Lc")) cutobjname="LctopKpiAnalysisCuts";
81 }
82 }
83 }
84 }
85 }
86
87 TFile* f=new TFile(filename.Data());
88 if(!f){
89 cout<<filename.Data()<<" not found"<<endl;
90 return kFALSE;
91 }
92 TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname);
93 if(!f){
94 cout<<dirname.Data()<<" not found in "<<filename.Data()<<endl;
95 return kFALSE;
96 }
97
98 list=(TList*)dir->Get(listname);
99 if(!list){
100 cout<<"List "<<listname.Data()<<" not found"<<endl;
101 dir->ls();
102 return kFALSE;
103 }
104
105 hstat=(TH1F*)dir->Get(hstatname);
106 if(!hstat){
107 cout<<hstatname.Data()<<" not found"<<endl;
108 return kFALSE;
109 }
110
111 cutobj=(AliRDHFCuts*)dir->Get(cutobjname);
112 if(!cutobj){
113 cout<<cutobjname.Data()<<" not found"<<endl;
114 return kFALSE;
115 }
116
9af24f46 117 return kTRUE;
118}
119
120//draw "track related" histograms (list "outputTrack")
1fc8d3c0 121void DrawOutputTrack(TString partname="D0",TString textleg="",TString path="./", Bool_t superimpose=kFALSE){
9af24f46 122 gStyle->SetCanvasColor(0);
123 gStyle->SetTitleFillColor(0);
124 gStyle->SetStatColor(0);
125 gStyle->SetPalette(1);
126
04fd2703 127 TString listname="outputTrack",name1="",name2="",path2="",filename=/*"AnalysisResults.root"*/"PWG3histograms.root",filename2="PWG3histograms.root";
1fc8d3c0 128 TString tmp="y";
129
130 if(superimpose){
131 cout<<"Enter the names:\n>";
132 cin>>name1;
133 cout<<">";
134 cin>>name2;
135 cout<<"Are they in the same output file? (y/n)"<<endl;
136 cin>>tmp;
137 if(tmp=="n"){
04fd2703 138 cout<<"Path: \n";
139 cout<<">";
1fc8d3c0 140 cin>>path2;
04fd2703 141 cout<<"Filename: "<<endl;
142 cout<<">";
1fc8d3c0 143 cin>>filename2;
144 }
145
146 }
9af24f46 147
148 TList* list;
149 TH1F * hstat;
150
1fc8d3c0 151 Bool_t isRead=ReadFile(list,hstat,listname,Form("%s%s",partname.Data(),name1.Data()),path,filename);
9af24f46 152 if(!isRead) return;
153 if(!list || !hstat){
154 cout<<":-( null pointers..."<<endl;
155 return;
156 }
1fc8d3c0 157 TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
158 pvtxt->SetBorderSize(0);
159 pvtxt->SetFillStyle(0);
160 pvtxt->AddText(name1);
161
162 TList* llist;
163 TH1F* hhstat;
164 if(superimpose){
165 isRead=ReadFile(llist,hhstat,listname,Form("%s%s",partname.Data(),name2.Data()),path2,filename2);
166 if(!isRead) return;
167 if(!llist || !hhstat){
168 cout<<":-( null pointers..."<<endl;
169 return;
170 }
171 TText *redtext=pvtxt->AddText(name2);
172 redtext->SetTextColor(kRed);
04fd2703 173 hhstat->Scale(hstat->Integral()/hhstat->Integral());
1fc8d3c0 174
175 }
9af24f46 176
177 for(Int_t i=0;i<list->GetEntries();i++){
178 TH1F* h=(TH1F*)list->At(i);
1fc8d3c0 179 TH1F* hh=0x0;
52bcee7a 180 TH1F* hr=0x0;
1fc8d3c0 181 if(superimpose){
182 hh=(TH1F*)llist->At(i);
52bcee7a 183 hr=(TH1F*)hh->Clone(Form("$s_ratio",hh->GetName()));
04fd2703 184 hh->Scale(h->Integral()/hh->Integral());
1fc8d3c0 185 }
186 if(!h || (superimpose && !hh)){
9af24f46 187 cout<<"Histogram "<<i<<" not found"<<endl;
188 continue;
189 }
1fc8d3c0 190 if(superimpose){
191 hhstat->SetLineColor(kRed);
192 hh->SetLineColor(kRed);
52bcee7a 193 hr->Divide(h);
1fc8d3c0 194 }
195
9af24f46 196 TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
197 c->cd();
198 c->SetGrid();
199 TString hname=h->GetName();
04fd2703 200 if(!hname.Contains("nCls")){
9af24f46 201 c->SetLogy();
1fc8d3c0 202 //h->SetMinimum(1);
9af24f46 203 h->Draw();
52bcee7a 204 if(superimpose) {
205 hh->Draw("sames");
206 TCanvas* c2=new TCanvas(Form("c2%s",h->GetName()),h->GetName());
207 c2->cd();
208 c2->SetGrid();
209 hr->Draw();
210 c2->SaveAs(Form("%s%s%s%sRatio.png",c->GetName(),name1.Data(),name2.Data(),textleg.Data()));
211 }
1fc8d3c0 212 } else {
213 h->Draw("htext0");
214 if(superimpose)hh->Draw("htext0sames");
215 }
52bcee7a 216 c->cd();
1fc8d3c0 217 pvtxt->Draw();
218 c->SaveAs(Form("%s%s%s%s.png",c->GetName(),name1.Data(),name2.Data(),textleg.Data()));
9af24f46 219 }
220
221 TCanvas* cst=new TCanvas("cst","Stat");
222 cst->SetGridy();
223 cst->cd();
224 hstat->Draw("htext0");
1fc8d3c0 225 if(superimpose) {
226 hhstat->Draw("htext0sames");
227 pvtxt->Draw();
228 }
9af24f46 229 cst->SaveAs(Form("%s%s.png",hstat->GetName(),textleg.Data()));
230
9af24f46 231
82f89d0d 232}
9af24f46 233
82f89d0d 234//draw "pid related" histograms (list "outputPID")
1fc8d3c0 235//period=-999 to draw the pull instead of the cut
236void DrawOutputPID(TString partname="D0", Int_t mode=0/*0=with pull, 1=with nsigma*/,TString textleg="",TString path="./"){
82f89d0d 237 gStyle->SetCanvasColor(0);
238 gStyle->SetTitleFillColor(0);
239 gStyle->SetStatColor(0);
240 gStyle->SetPalette(1);
241
1fc8d3c0 242 Int_t period=2 ,set=0;
243 if(mode==1){
04fd2703 244 cout<<"Choose period: \n-LHC10h -> 2;\n-LHC10de -> 1;\n-LHC10bc -> 0"<<endl;
1fc8d3c0 245 cin>>period;
246 if(period>0){
247 cout<<"Choose set: "<<endl;
248 if(period==2) cout<<"-pass1 -> 0;\n-pass2 -> 1"<<endl;
249 cin>>set;
250 }
251 }
252
82f89d0d 253 TString listname="outputPid";
254
255 TList* list;
256 TH1F * hstat;
1fc8d3c0 257 //needed only for mode 1
258 AliRDHFCuts* cutobj;
259 AliAODPidHF* aodpid;
260 Double_t nsigmaTOF=0;
261 Double_t nsigmaTPC[3]={},plimTPC[2]={};
262
263 if(mode==1){
264 Bool_t isRead=ReadFileMore(list,hstat,cutobj,listname,partname,path);
265 if(!isRead) return;
266 if(!list || !hstat){
267 cout<<":-( null pointers..."<<endl;
268 return;
269 }
270 aodpid=(AliAODPidHF*)cutobj->GetPidHF();
04fd2703 271 if(!aodpid){
272 cout<<"PidHF object not found! cannot get the nsigma values"<<endl;
273 return;
274 }
1fc8d3c0 275 nsigmaTOF=aodpid->GetSigma(3);
276
277 nsigmaTPC[0]=aodpid->GetSigma(0);
278 nsigmaTPC[1]=aodpid->GetSigma(1);
279 nsigmaTPC[2]=aodpid->GetSigma(2);
280 aodpid->GetPLimit(plimTPC);
281
282 }else{
283
284 Bool_t isRead=ReadFile(list,hstat,listname,partname,path);
285 if(!isRead) return;
286 if(!list || !hstat){
287 cout<<":-( null pointers..."<<endl;
288 return;
289 }
9af24f46 290 }
291
1fc8d3c0 292
293 TPaveText *txtsigmaTOF=new TPaveText(0.1,0.65,0.5,0.9,"NDC");
294 txtsigmaTOF->SetBorderSize(0);
295 txtsigmaTOF->SetFillStyle(0);
296 txtsigmaTOF->AddText(Form("nsigmacut from cutobj = %.1f",nsigmaTOF));
297 TLine lTOF;
298 lTOF.SetLineColor(kMagenta+1);
299 lTOF.SetLineStyle(2);
300 lTOF.SetLineWidth(3);
301
302 TPaveText *txtsigmaTPC=new TPaveText(0.3,0.6,0.6,0.9,"NDC");
303 txtsigmaTPC->SetBorderSize(0);
304 txtsigmaTPC->SetFillStyle(0);
305 txtsigmaTPC->AddText("nsigmacut from cutobj \n");
306 txtsigmaTPC->AddText(Form("p < %.1f : %.1f \n",plimTPC[0],nsigmaTPC[0]));
307 txtsigmaTPC->AddText(Form("%.1f < p < %.1f : %.1f \n",plimTPC[0],plimTPC[1],nsigmaTPC[1]));
308 txtsigmaTPC->AddText(Form("p > %.1f : %.1f \n",plimTPC[1],nsigmaTPC[2]));
309 TLine lTPC;
310 lTPC.SetLineColor(kMagenta+1);
311 lTPC.SetLineStyle(2);
312 lTPC.SetLineWidth(3);
313
314 // TCanvas *ctest=new TCanvas("text","Test text");
315 // ctest->cd();
316 // txtsigmaTPC->Draw();
317 // txtsigmaTOF->Draw();
318
319
82f89d0d 320 for(Int_t i=0;i<list->GetEntries();i++){
321 TClass* objtype=list->At(i)->IsA();
322 TString tpname=objtype->GetName();
323
324 if(tpname=="TH1F"){
325 TH1F* h=(TH1F*)list->At(i);
326
327 if(!h){
328 cout<<"Histogram "<<i<<" not found"<<endl;
329 continue;
330 }
331 //h->Scale(1./h->Integral("width"));
332 TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
333 c->SetLogz();
334 c->cd();
335 h->Draw();
336
337 //write
1ff03d38 338 c->SaveAs(Form("%s%s.png",h->GetName(),textleg.Data()));
82f89d0d 339 TFile* fout=new TFile(Form("%s.root",h->GetName()),"recreate");
340 fout->cd();
341 c->Write();
342 }
343
344 if(tpname=="TH2F"){
345 TH2F* h=(TH2F*)list->At(i);
346
347 if(!h){
348 cout<<"Histogram "<<i<<" not found"<<endl;
349 continue;
350 }
ac0c2841 351 h->Sumw2();
82f89d0d 352 h->Scale(1./h->Integral("width"));
1fc8d3c0 353 TString hname=h->GetName();
354
355 if(hname.Contains("hTOFtimeKaonHyptime")){
356 TCanvas* cz=new TCanvas(Form("c%szoom",hname.Data()),Form("%szoom",hname.Data()));
357 cz->SetLogz();
358 TH2F* hz=(TH2F*)h->Clone(Form("%sz",hname.Data()));
359 hz->Draw("colz");
360 hz->SetAxisRange(-1500,1500,"Y");
361 //write
362 cz->SaveAs(Form("%szoom.png",h->GetName()));
363 }
364
365 TCanvas* c=new TCanvas(Form("c%s",hname.Data()),hname.Data());
82f89d0d 366 c->SetLogz();
ac0c2841 367 //c->SetLogx();
82f89d0d 368 c->cd();
369
370 h->Draw("colz");
1fc8d3c0 371
372 //TCanvas *test=new TCanvas("test","test");
373 if(mode==0){
374 //mean and pull, code from Jens Wiechula
375 TF1 fg("fg","gaus",-2.,2.); // fit range +- 2 sigma
376 TLine l;
377 TObjArray arr;
378
379 //h->Draw("colz");
380 fg.SetParameters(1,0,1);
381 h->FitSlicesY(&fg,0,-1,0,"NQR",&arr);
382
383 TH1 *hM=(TH1*)arr.At(1);
384 hM->SetMarkerStyle(20);
385 hM->SetMarkerSize(.5);
386 hM->DrawClone("sames");
387
388 TH1 *hS=(TH1*)arr.At(2);
389 hS->SetMarkerStyle(20);
390 hS->SetMarkerSize(.5);
391 hS->SetMarkerColor(kRed);
392 hS->SetLineColor(kRed);
393 hS->DrawClone("same");
394
395 l.SetLineColor(kBlack);
396 l.DrawLine(.2,0,20,0);
397 l.SetLineColor(kRed);
398 l.DrawLine(.2,1,20,1);
399
400 }else{ //mode 1
401
402 if(hname.Contains("TOFsigma")) {
403
404 c->cd();
405 txtsigmaTOF->Draw();
406 lTOF.DrawLine(.2,nsigmaTOF,20,nsigmaTOF);
407 lTOF.DrawLine(.2,-1*nsigmaTOF,4.,-1*nsigmaTOF);
ac0c2841 408
1fc8d3c0 409 }
410
ac0c2841 411
1fc8d3c0 412 if(hname.Contains("TPCsigma")){
413
414 c->cd();
415 txtsigmaTPC->Draw();
416 lTPC.DrawLine(0.,nsigmaTPC[0],plimTPC[0],nsigmaTPC[0]);
417 lTPC.DrawLine(plimTPC[0],nsigmaTPC[1],plimTPC[1],nsigmaTPC[1]);
418 lTPC.DrawLine(plimTPC[1],nsigmaTPC[2],4,nsigmaTPC[2]);
419 lTPC.DrawLine(0.,-1*nsigmaTPC[0],plimTPC[0],-1*nsigmaTPC[0]);
420 lTPC.DrawLine(plimTPC[0],-1*nsigmaTPC[1],plimTPC[1],-1*nsigmaTPC[1]);
421 lTPC.DrawLine(plimTPC[1],-1*nsigmaTPC[2],4,-1*nsigmaTPC[2]);
422 }
ac0c2841 423
1fc8d3c0 424 if(hname.Contains("TPCsigvsp")){
425 SuperimposeBBToTPCSignal(period,c,set);
426 }
427 }
428
82f89d0d 429 //write
f3e1ad2f 430 c->SaveAs(Form("%s%d.png",h->GetName(),mode));
431 TFile* fout=new TFile(Form("%s%d.root",h->GetName(),mode),"recreate");
82f89d0d 432 fout->cd();
433 c->Write();
1fc8d3c0 434 }
435 }
436}
437
438void SuperimposeBBToTPCSignal(Int_t period /*0=LHC10bc, 1=LHC10d, 2=LHC10h*/,TCanvas* cpid,Int_t set /*see below*/){
439
440 TFile* fBethe=new TFile("BetheBlochTPC.root");
441 if(!fBethe->IsOpen()){
442 TPCBetheBloch(set);
443 fBethe=new TFile("BetheBlochTPC.root");
444 }
445 const Int_t npart=4;
446 TString partnames[npart]={"Kaon","Pion","Electron","Proton"};
447 for(Int_t ipart=0;ipart<npart;ipart++){
448 TString grname=Form("%sP%d",partnames[ipart].Data(),period);
449 TGraph* gr=(TGraph*)fBethe->Get(grname);
450 cpid->cd();
451 gr->SetLineColor(1);
452 gr->SetLineWidth(2);
453 gr->Draw("L");
454 }
455
456 //cpid->SaveAs(Form("%sBB.png",hname.Data()));
457}
458
459//draw and save Bethe Bloch from TPC in different periods
460void TPCBetheBloch(Int_t set){
461 gStyle->SetOptTitle(0);
462 gStyle->SetCanvasColor(0);
463
464 AliTPCPIDResponse *tpcResp=new AliTPCPIDResponse();
465
466 const Int_t npart=4;
467 Double_t masses[npart]={TDatabasePDG::Instance()->GetParticle(321)->Mass()/*Kaon*/,TDatabasePDG::Instance()->GetParticle(211)->Mass()/*Pion*/,TDatabasePDG::Instance()->GetParticle(11)->Mass()/*Electron*/,TDatabasePDG::Instance()->GetParticle(2212)->Mass()/*Proton*/};
468 TString partnames[npart]={"Kaon","Pion","Electron","Proton"};
469 //printf("%s = %.4f,%s = %.4f,%s = %.4f\n",partnames[0].Data(),masses[0],partnames[1].Data(),masses[1],partnames[2].Data(),masses[2]);
470 TCanvas *cBethe=new TCanvas("cBethe","Bethe Bloch K pi e p");
04fd2703 471 Int_t nperiods=4; //LHC10b+c, LHC10d, LHC10h, MC
1fc8d3c0 472 Double_t alephParameters[5]={};
04fd2703 473 Int_t nsets=1/*LHC10bc*/+2/*LHC10de*/+2/*LHC10h*/+3/*MC*/;
1fc8d3c0 474
04fd2703 475 periodsname=new TString[nsets];
476 cout<<"Creating the file of the Bethe Bloch"<<endl;
1fc8d3c0 477 TFile* fout=new TFile("BetheBlochTPC.root","recreate");
478
479 for(Int_t iperiod=0;iperiod<nperiods;iperiod++){
04fd2703 480 cout<<"Period "<<iperiod<<" : ";
481 if(iperiod==0){ //LHC10bc
482
483 alephParameters[0] = 0.0283086/0.97;
484 alephParameters[1] = 2.63394e+01;
485 alephParameters[2] = 5.04114e-11;
486 alephParameters[3] = 2.12543e+00;
487 alephParameters[4] = 4.88663e+00;
488 periodsname[0]="dataLHC10bc";
1fc8d3c0 489 }
04fd2703 490 if(iperiod==1){ //LHC10de,low energy
1fc8d3c0 491 if(set==0){
492 alephParameters[0] = 1.63246/50.;
493 alephParameters[1] = 2.20028e+01;
494 alephParameters[2] = TMath::Exp(-2.48879e+01);
495 alephParameters[3] = 2.39804e+00;
496 alephParameters[4] = 5.12090e+00;
04fd2703 497 periodsname[1]="dataLHC10deold";
1fc8d3c0 498 }
499 if(set==1){
500 alephParameters[0] = 1.34490e+00/50.;
501 alephParameters[1] = 2.69455e+01;
502 alephParameters[2] = TMath::Exp(-2.97552e+01);
503 alephParameters[3] = 2.35339e+00;
504 alephParameters[4] = 5.98079e+00;
04fd2703 505 periodsname[2]="dataLHC10denew";
1fc8d3c0 506 }
507 }
04fd2703 508
509 if(iperiod==2){//LHC10h
510 if(set==0){//pass1
511 alephParameters[0]=1.25202/50.;
512 alephParameters[1]=2.74992e+01;
513 alephParameters[2]=TMath::Exp(-3.31517e+01);
514 alephParameters[3]=2.46246;
515 alephParameters[4]=6.78938;
516 periodsname[3]="dataLHC10hpass1";
517 }
518 if (set==1){//pass2 (AOD044)
519 alephParameters[0]=1.25202/50.;
520 alephParameters[1]=2.74992e+01;
521 alephParameters[2]=TMath::Exp(-3.31517e+01);
522 alephParameters[3]=2.46246;
523 alephParameters[4]=6.78938;
524 periodsname[4]="dataLHC10hpass2";
525 }
1fc8d3c0 526 }
04fd2703 527 if(iperiod==3){ //MC
f3e1ad2f 528 if(set==0){
529 alephParameters[0] = 2.15898e+00/50.;
530 alephParameters[1] = 1.75295e+01;
531 alephParameters[2] = 3.40030e-09;
532 alephParameters[3] = 1.96178e+00;
533 alephParameters[4] = 3.91720e+00;
04fd2703 534 periodsname[5]="MCold";
f3e1ad2f 535 }
536 if(set==1){ //new
537 alephParameters[0] = 1.44405/50;
538 alephParameters[1] = 2.35409e+01;
539 alephParameters[2] = TMath::Exp(-2.90330e+01);
540 alephParameters[3] = 2.10681;
541 alephParameters[4] = 4.62254;
04fd2703 542 periodsname[6]="MCnew";
543 }
544
545 if(set==2){ //new BB from Francesco
546 alephParameters[0] = 0.029021;
547 alephParameters[1] = 25.4181;
548 alephParameters[2] = 4.66596e-08;
549 alephParameters[3] = 1.90008;
550 alephParameters[4] = 4.63783;
551 periodsname[7]="MCBBFrancesco";
f3e1ad2f 552 }
04fd2703 553
554 if(set==3){ //low energy 2011
555 alephParameters[0] = 0.0207667;
556 alephParameters[1] = 29.9936;
557 alephParameters[2] = 3.87866e-11;
558 alephParameters[3] = 2.17291;
559 alephParameters[4] = 7.1623;
560 //periodsname[8]="MClowen2011";
561 }
562
563
f3e1ad2f 564 }
04fd2703 565 //cout<<periodsname[iperiod]<<endl;
1fc8d3c0 566 tpcResp->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
04fd2703 567 cout<<"here"<<endl;
1fc8d3c0 568 for(Int_t ipart=0;ipart<npart;ipart++){
ac0c2841 569
1fc8d3c0 570 const Int_t n=1000;
571 Double_t p[n],bethe[n];
572
573 for(Int_t k=0;k<n;k++){ //loop on the momentum steps
574 p[k]=0.0001+k*4./n; //limits 0.-4. GeV/c
575 //cout<<p[k]<<"\t";
576 //bethe[k]=-tpcResp->Bethe(p[k]/masses[ipart]);
577 AliPID::EParticleType ptype=AliPID::kKaon;
578 if(ipart==1) ptype=AliPID::kPion;
579 if(ipart==2) ptype=AliPID::kElectron;
580 if(ipart==3) ptype=AliPID::kProton;
581 bethe[k]=tpcResp->GetExpectedSignal(p[k],ptype);
582 }
583 //cout<<endl;
584 TGraph *gr=new TGraph(n,p,bethe);
585 gr->SetName(Form("%sP%d",partnames[ipart].Data(),iperiod));
586 gr->SetTitle(Form("%sP%d;p (GeV/c);",partnames[ipart].Data(),iperiod));
587 gr->SetLineColor(ipart+1);
588 gr->SetMarkerColor(ipart+1);
589 gr->GetYaxis()->SetRangeUser(35,100);
590 cBethe->cd();
591 if(iperiod==0 && ipart==0)gr->DrawClone("AL");
592 else gr->DrawClone("L");
593
594 fout->cd();
595 gr->Write();
82f89d0d 596 }
1fc8d3c0 597
9af24f46 598 }
1fc8d3c0 599 TParameter<int> sett;
600 sett.SetVal(set);
601 fout->cd();
602 sett.Write();
603
604 fout->Close();
9af24f46 605}
606
1fc8d3c0 607void DrawOutputCentrality(TString partname="D0",TString textleg="",TString path="./", Bool_t superimpose=kFALSE){
9af24f46 608 gStyle->SetCanvasColor(0);
609 gStyle->SetTitleFillColor(0);
610 gStyle->SetStatColor(0);
611 gStyle->SetPalette(1);
612
04fd2703 613 TString listname="outputCentrCheck",name1="",name2="",path2="",filename=/*"AnalysisResults.root"*/"PWG3histograms.root",filename2="PWG3histograms.root";
614
615 // if(superimpose){
616 // cout<<"Enter the names:\n>";
617 // cin>>name1;
618 // cout<<">";
619 // cin>>name2;
620 // }
621 // TString listname="outputTrack",name1="",name2="";
622 TString tmp="y";
1fc8d3c0 623
624 if(superimpose){
625 cout<<"Enter the names:\n>";
626 cin>>name1;
627 cout<<">";
628 cin>>name2;
04fd2703 629 cout<<"Are they in the same output file? (y/n)"<<endl;
630 cin>>tmp;
631 if(tmp=="n"){
632 cout<<"Path: \n";
633 cout<<">";
634 cin>>path2;
635 cout<<"Filename: "<<endl;
636 cout<<">";
637 cin>>filename2;
638 }
639
1fc8d3c0 640 }
641 // Int_t nhist=1;
642 // TString *name=0x0;
643 // if(superimpose){
644 // cout<<"Number of histogram to superimpose: ";
645 // cin>>nhist;
646 // name=new TString[nhist];
647 // for (Int_t j=0;j<nhist;j++){
648 // cout<<">";
649 // cin>>name[j];
650 // }
651 // }
9af24f46 652
653 TList* list;
654 TH1F * hstat;
655
1fc8d3c0 656 Bool_t isRead=ReadFile(list,hstat,listname,Form("%s%s",partname.Data(),name1.Data()),path);
9af24f46 657 if(!isRead) return;
658 if(!list || !hstat){
659 cout<<":-( null pointers..."<<endl;
660 return;
661 }
1ff03d38 662
1fc8d3c0 663 TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
664 pvtxt->SetBorderSize(0);
665 pvtxt->SetFillStyle(0);
666 pvtxt->AddText(name1);
667
668 TList* llist;
669 TH1F* hhstat;
670 if(superimpose){
04fd2703 671 isRead=ReadFile(llist,hhstat,listname,Form("%s%s",partname.Data(),name2.Data()),path2);
1fc8d3c0 672 if(!isRead) return;
673 if(!llist || !hhstat){
674 cout<<":-( null pointers..."<<endl;
675 return;
676 }
677 TText *redtext=pvtxt->AddText(name2);
678 redtext->SetTextColor(kRed);
679
680 }
681
682
82f89d0d 683 TCanvas* cst=new TCanvas("cst","Stat");
684 cst->SetGridy();
685 cst->cd();
1ff03d38 686 Int_t nevents=hstat->Integral(1,2);
82f89d0d 687 hstat->Draw("htext0");
688 cst->SaveAs(Form("%s%s.png",hstat->GetName(),textleg.Data()));
04fd2703 689 Int_t nevents080=1,nnevents080=1;
1ff03d38 690
ac0c2841 691 //TCanvas *spare=new TCanvas("sparecv","Spare");
692
1ff03d38 693 for(Int_t i=0;i<list->GetEntries();i++){
694
695 TClass* objtype=list->At(i)->IsA();
696 TString tpname=objtype->GetName();
697
698 if(tpname=="TH1F"){
699
700 TH1F* h=(TH1F*)list->At(i);
1fc8d3c0 701 TH1F* hh=0x0;
702 if(superimpose){
703 hh=(TH1F*)llist->At(i);
704 }
705 if(!h || (superimpose && !hh)){
1ff03d38 706 cout<<"Histogram "<<i<<" not found"<<endl;
707 continue;
708 }
1fc8d3c0 709 if(superimpose){
710 hhstat->SetLineColor(kRed);
711 hh->SetLineColor(kRed);
712 }
713
1ff03d38 714 TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
1fc8d3c0 715 TPaveText *pvtxt2=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
716 pvtxt2->SetBorderSize(0);
717 pvtxt2->SetFillStyle(0);
1ff03d38 718
719 c->cd();
720 c->SetGrid();
1fc8d3c0 721 c->SetLogy();
1ff03d38 722 Int_t entries=h->Integral();
04fd2703 723 pvtxt2->AddText(Form("%.1f %s of the events",(Double_t)entries/(Double_t)nevents*100,"%"));
1ff03d38 724 h->Draw();
1fc8d3c0 725 if(superimpose) {
726 hh->Draw("sames");
727 pvtxt->Draw();
728 }
729 pvtxt2->Draw();
1ff03d38 730 c->SaveAs(Form("%s%s.png",c->GetName(),textleg.Data()));
731 }
732 if(tpname=="TH2F"){
733 TH2F* h=(TH2F*)list->At(i);
734 if(!h){
735 cout<<"Histogram "<<i<<" not found"<<endl;
736 continue;
737 }
738 TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
739 TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
740 pvtxt->SetBorderSize(0);
741 pvtxt->SetFillStyle(0);
742
743 c->cd();
744 c->SetGrid();
745 Int_t entries=h->Integral();
04fd2703 746 pvtxt->AddText(Form("%.1f %s of the events",(Double_t)entries/(Double_t)nevents*100,"%"));
1ff03d38 747 h->Draw("colz");
748 c->SetLogz();
749 pvtxt->Draw();
750 c->SaveAs(Form("%s%s.png",c->GetName(),textleg.Data()));
751 }
752 }
753
82f89d0d 754
755 listname="countersCentrality";
756
04fd2703 757 isRead=ReadFile(list,hstat,listname,Form("%s%s",partname.Data(),name1.Data()),path);
82f89d0d 758 if(!isRead) return;
759 if(!list || !hstat){
760 cout<<":-( null pointers..."<<endl;
761 return;
762 }
1ff03d38 763
04fd2703 764
765 if(superimpose){
766 isRead=ReadFile(llist,hhstat,listname,Form("%s%s",partname.Data(),name2.Data()),path2);
767 if(!isRead) return;
768 if(!llist || !hhstat){
769 cout<<":-( null pointers..."<<endl;
770 return;
771 }
772 TText *redtext=pvtxt->AddText(name2);
773 redtext->SetTextColor(kRed);
774
775 }
776
1ff03d38 777 TH1F* hallcntr=0x0;
04fd2703 778 TH1F* hhallcntr=0x0;
1ff03d38 779 cout<<"normalizing to 0-80% as a check"<<endl;
780
781 TCanvas *cvnocnt=new TCanvas("cvnocnt","No Centrality estimation",800,400);
782 cvnocnt->Divide(2,1);
783
82f89d0d 784 for(Int_t i=0;i<list->GetEntries();i++){
785 AliCounterCollection* coll=(AliCounterCollection*)list->At(i);
04fd2703 786 AliCounterCollection* colle=0x0;
787 if(superimpose) colle=(AliCounterCollection*)llist->At(i);
82f89d0d 788 coll->SortRubric("run");//sort by run number
789 //coll->PrintKeyWords();
790 Int_t ncentr=10;//check this
791 TH1F* h020=0x0;
792 TH1F* h2080=0x0;
04fd2703 793 TH1F* hh020=0x0;
794 TH1F* hh2080=0x0;
1ff03d38 795 hallcntr=0x0;
04fd2703 796 hhallcntr=0x0;
1ff03d38 797
798 TH1F* hbad=(TH1F*)coll->Get("run",Form("centralityclass:-990_-980"));
799 cvnocnt->cd(i+1);
800 if(hbad) hbad->Draw();
801
82f89d0d 802 TCanvas *ccent=new TCanvas(Form("ccent%s",coll->GetName()),Form("Centrality vs Run (%s)",coll->GetName()),1400,800);
803 ccent->Divide(5,3);
04fd2703 804
805 TH1F* hh=0x0;
82f89d0d 806
1ff03d38 807 for(Int_t ic=0;ic<8/*ncentr*/;ic++){ //normalizing to 0-80% as a check
808
809 TH1F* h=(TH1F*)coll->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
810 h->SetName(Form("h%d%d",i,ic));
811 if(!hallcntr) {
812 hallcntr=(TH1F*)h->Clone("hallcntr");
813 hallcntr->Sumw2();
04fd2703 814 } else {
815 hallcntr->Add(h);
816 }
1fc8d3c0 817
1ff03d38 818 nevents080+=h->Integral();
04fd2703 819
820 if(superimpose){
821 hh=(TH1F*)colle->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
822 hh->SetName(Form("hh%d%d",i,ic));
823 if(!hhallcntr) {
824 hhallcntr=(TH1F*)hh->Clone("hhallcntr");
825 hhallcntr->Sumw2();
826 }else hhallcntr->Add(hh);
827 nnevents080+=hh->Integral();
828
829 }
1ff03d38 830 }
831
82f89d0d 832 for(Int_t ic=0;ic<ncentr;ic++){
833
834 TH1F* h=(TH1F*)coll->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
835 h->SetName(Form("h%d%d",i,ic));
1ff03d38 836 h->Sumw2();
837
82f89d0d 838 if(ic>=0 && ic<=1){ //0-20
839 if(!h020) {
840 h020=(TH1F*)h->Clone(Form("h020%s",coll->GetName()));
841 h020->SetTitle(Form("Centrality 0-20 %s",coll->GetName()));
04fd2703 842 if(superimpose){
843 hh020=(TH1F*)hh->Clone(Form("hh020%s",coll->GetName()));
844 hh020->SetTitle(Form("Centrality 0-20 %s",coll->GetName()));
845 }
846 }
847 else {
848 h020->Add(h);
849 if(superimpose)hh020->Add(hh);
82f89d0d 850 }
82f89d0d 851 }
852 if(ic>=2 && ic<=7){ //20-80
853 if(!h2080) {
854 h2080=(TH1F*)h->Clone(Form("h2080%s",coll->GetName()));
855 h2080->SetTitle(Form("Centrality 20-80 %s",coll->GetName()));
04fd2703 856 if(superimpose){
857 hh2080=(TH1F*)hh->Clone(Form("hh2080%s",coll->GetName()));
858 hh2080->SetTitle(Form("Centrality 20-80 %s",coll->GetName()));
859 }
82f89d0d 860 }
04fd2703 861 else {
862 h2080->Add(h);
863 if(superimpose)hh2080->Add(hh);
864 }
865
82f89d0d 866 }
867
1ff03d38 868 h->Divide(hallcntr);
869
82f89d0d 870 ccent->cd(ic+1);
1ff03d38 871 h->GetYaxis()->SetRangeUser(0.,0.15);
82f89d0d 872 h->DrawClone();
ac0c2841 873 /*
874 if(ic==0&&i==0){
875 spare->cd();
876 h->Draw();
877 }
878 */
82f89d0d 879 // ccent->cd(1);
880 // h->SetLineColor(ic+1);
881 // if(ic==0)h->DrawClone();
882 // else h->DrawClone("sames");
883 }
884
885 ccent->cd(ncentr+1);
1ff03d38 886 h020->Divide(hallcntr);
887 h020->DrawClone();
04fd2703 888 if(superimpose){
889 hh020->Divide(hhallcntr);
890 hh020->SetLineColor(2);
891 hh020->SetMarkerColor(2);
892 hh020->DrawClone("sames");
893 }
1fc8d3c0 894 TCanvas* cv020=new TCanvas(Form("cv020-%d",i),"0-20% vs run number",1400,600);
1ff03d38 895 cv020->cd();
896 h020->GetYaxis()->SetRangeUser(0.,1.);
82f89d0d 897 h020->DrawClone();
04fd2703 898 if(superimpose)hh020->DrawClone("sames");
1ff03d38 899 cv020->SaveAs(Form("cv020-%d.png",i));
82f89d0d 900
901 ccent->cd(ncentr+2);
1ff03d38 902 h2080->Divide(hallcntr);
82f89d0d 903 h2080->DrawClone();
04fd2703 904 if(superimpose){
905 hh2080->Divide(hhallcntr);
906 hh2080->SetLineColor(2);
907 hh2080->SetMarkerColor(2);
908 hh2080->DrawClone("sames");
909 }
82f89d0d 910
1fc8d3c0 911 TCanvas* cv2080=new TCanvas(Form("cv2080-%d",i),"20-80% vs run number",1400,600);
1ff03d38 912 cv2080->cd();
913 h2080->GetYaxis()->SetRangeUser(0.,1.);
914 h2080->DrawClone();
04fd2703 915 if(superimpose)hh2080->DrawClone("sames");
1ff03d38 916 cv2080->SaveAs(Form("cv2080-%d.png",i));
917
82f89d0d 918 ccent->SaveAs(Form("%s%s.png",ccent->GetName(),textleg.Data()));
9af24f46 919 }
920
921}
1ff03d38 922
923void DrawProjections(TString partname="D0",TString h2dname="hMultvsPercentile",Int_t nsteps=0,TString direction="X",TString path="./"){
924 gStyle->SetCanvasColor(0);
925 gStyle->SetTitleFillColor(0);
926 gStyle->SetStatColor(0);
29df3603 927 gStyle->SetPalette(1);
1ff03d38 928
929 TString listname="outputCentrCheck";
930
931 TList* list;
932 TH1F * hstat;
933
934 Bool_t isRead=ReadFile(list,hstat,listname,partname,path);
935 if(!isRead) return;
936 if(!list || !hstat){
937 cout<<":-( null pointers..."<<endl;
938 return;
939 }
1fc8d3c0 940 Double_t nevents=hstat->Integral(5,6); //ev good vertex
1ff03d38 941
942 TH2F* h2=(TH2F*)list->FindObject(h2dname);
943 if(!h2){
944 cout<<h2dname.Data()<<" not found"<<endl;
945 return;
946 }
29df3603 947 TCanvas* cv2d=new TCanvas("cv2d",h2->GetName());
948 cv2d->cd();
949 cv2d->SetLogz();
1fc8d3c0 950 cv2d->SetGrid();
29df3603 951 h2->Draw("colz");
1fc8d3c0 952 TPaveText *pvst=new TPaveText(0.6,0.2,0.9,0.7,"NDC");
953 pvst->SetBorderSize(0);
954 pvst->SetFillStyle(0);
955 pvst->AddText("Bin -> Cont/nEvVtx");
956
29df3603 957
1ff03d38 958 Int_t kbins=1;
959 if(nsteps==0){
960 if(direction=="X") nsteps=h2->GetNbinsY();
961 if(direction=="Y") nsteps=h2->GetNbinsX();
962 }
963 else{
964 if(direction=="X") kbins=h2->GetNbinsY()/nsteps;
965 if(direction=="Y") kbins=h2->GetNbinsX()/nsteps;
966 }
967
968 TCanvas *cvpj=new TCanvas(Form("cvpj%s%s",direction.Data(),h2dname.Data()),Form("cvpj%s",direction.Data()),1000,800);
969 cvpj->Divide((Int_t)(nsteps/3)+1,3);
970 TFile* fout=new TFile(Form("proj%s%s.root",direction.Data(),h2dname.Data()), "recreate");
1fc8d3c0 971 //Float_t maxx[nsteps];
972 Float_t maxx[12]={9000,9000,6000,4000,2000,1400,800,500,200,100,40,25};
973 Double_t integralpernev[nsteps];
974
1ff03d38 975 for(Int_t i=0;i<nsteps;i++){
976 TH1F* h=0x0;
977 if(direction=="X")h=(TH1F*)h2->ProjectionX(Form("px%d",i),i+kbins,i+2*kbins);
978 if(direction=="Y")h=(TH1F*)h2->ProjectionY(Form("py%d",i),i+kbins,i+2*kbins);
1fc8d3c0 979 integralpernev[i]=h->Integral()/nevents;
ac0c2841 980
981 TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
982 pvtxt->SetBorderSize(0);
983 pvtxt->SetFillStyle(0);
1fc8d3c0 984 pvtxt->AddText(Form("%.0f - %.0f",h2->GetYaxis()->GetBinLowEdge((i+kbins)),h2->GetYaxis()->GetBinLowEdge((i+2*kbins))));
985 pvst->AddText(Form("%.0f - %.0f -> %.2f",h2->GetYaxis()->GetBinLowEdge((i+kbins)),h2->GetYaxis()->GetBinLowEdge((i+2*kbins)),integralpernev[i]));
1ff03d38 986 cvpj->cd(i+1);
1fc8d3c0 987 h->GetXaxis()->SetRangeUser(0,maxx[i]);
1ff03d38 988 h->Draw();
ac0c2841 989 pvtxt->Draw();
1ff03d38 990 fout->cd();
991 h->Write();
992 }
993 cvpj->SaveAs(Form("cvpj%s%s.png",direction.Data(),h2dname.Data()));
994
1fc8d3c0 995 cv2d->cd();
996 pvst->Draw();
997 cv2d->SaveAs(Form("%s.png",h2->GetName()));
998
1ff03d38 999}