--- /dev/null
+void QAphi(char* system="pp276",char* name_fin="AnalysisResults.root",char* name_fout="QAphi",char* name_list="RsnHistMini_Phi_PhiNsigma_KTPCnsig30",char* name_hist_base="RsnMini_phi.PhiNsigma_KTPCnsig30"){
+ //This macro was written by Anders Knospe (anders.knospe@cern.ch).
+ //5 November 2013
+
+ //FOR DETAILED INSTRUCTIONS, SEE THE END OF THIS DOCUMENT.
+
+ //Arguments:
+ //system: string giving the collision system ("pp276","pp7", or "PbPb276"), used to figure out what the expected yield should be. If the system you want is unavailable, you will need to calculate your own expected yield.
+ //name_fin: name of input file
+ //name_fout: base name of output files (without suffix)
+ //name_list: name of the list in fin that contains histograms, may be different depending on the cuts applied to the kaon daughters
+ //name_hist_base: base name of the THnSparse histograms, may be different depending on the cuts applied to the kaon daughters
+
+ //This Macro does the following:
+ //1.) plots the standard histograms showing the number of events that pass selection of criteria and the number of events in each multiplicity or centrality bin
+ //2.) plots invariant-mass histograms (unlike-charge, like-charge, background-subtracted)
+ //3.) fits the background-subtracted invariant-mass distribution. In case there is a problem with the "default" fit, it does multiple fits using different residual backgrounds and fit regions. All fits are plotted and all fit results are saved.
+ //4.) prints out relevant information: pT range, phi yield per event, phi mass, and phi width (taken from the default fit)
+
+ gStyle->SetOptStat(0);
+
+ //----- get input histograms -----
+
+ TFile* fin=TFile::Open(name_fin);//open input file
+ if(!fin) return;
+
+ TList* l=(TList*) fin->Get(name_list);
+ if(!l){cerr<<"Error in QAphi(): missing input list "<<name_list<<" in file "<<name_fin<<". Stopping."<<endl; return;}
+
+ TH1D *hevent,*hmult,*hu,*hm,*hp,*hl,*hs,*hs_plot;
+ THnSparse *su,*sm,*sp;
+
+ hevent=(TH1D*) l->FindObject("hEventStat");//standard histogram in resonance package, shows number of events after various selections
+ if(!hevent){cerr<<"Error in QAphi(): missing input histogram hEventStat in file "<<name_fin<<". Stopping."<<endl; return;}
+
+ hmult=(TH1D*) l->FindObject("hAEventsVsMulti");//standard histogram in resonance package, shows number of events in each multiplicity (for pp) or centrality (for Pb-Pb) bin
+ if(!hmult){cerr<<"Error in QAphi(): missing input histogram hAEventsVsMulti in file "<<name_fin<<". Stopping."<<endl; return;}
+ double nevt=hmult->Integral();
+
+ su=(THnSparse*) l->FindObject(Form("%s_Unlike",name_hist_base));//invariant-mass histogram, unlike-charge
+ if(!su){cerr<<"Error in QAphi(): missing input histogram "<<name_hist_base<<"_Unlike in file "<<name_fin<<". Stopping."<<endl; return;}
+
+ sm=(THnSparse*) l->FindObject(Form("%s_LikeMM",name_hist_base));//invariant-mass histogram, like-charge (K-K-)
+ if(!sm){cerr<<"Error in QAphi(): missing input histogram "<<name_hist_base<<"_LikeMM in file "<<name_fin<<". Stopping."<<endl; return;}
+
+ sp=(THnSparse*) l->FindObject(Form("%s_LikePP",name_hist_base));//invariant-mass histogram, like-charge (K+K+)
+ if(!sp){cerr<<"Error in QAphi(): missing input histogram "<<name_hist_base<<"_LikePP in file "<<name_fin<<". Stopping."<<endl; return;}
+
+ TFile* fout=new TFile(Form("%s.root",name_fout),"RECREATE","HistoFile");//open output file
+
+ bool ptOK=SetPtRange(su);//Was the requested pT range set correctly?
+ SetPtRange(sm);
+ SetPtRange(sp);
+
+ double dy;
+ bool yOK=SetRapidityRange(su,dy);//Is the expected rapidity range (|y|<0.5) being used? Fill value of dy.
+ SetRapidityRange(sm,dy);
+ SetRapidityRange(sp,dy);
+
+ //----- plot the event and multiplicity/centrality histograms -----
+
+ TCanvas* c=new TCanvas("c","",10,10,1500,500);
+ c->SetFillColor(0);
+
+ TPad* p1=new TPad("p1","",0.,0.,1./3,1.);
+ p1->SetFillColor(0);
+
+ TPad* p2=new TPad("p2","",1./3,0.,2./3,1.);
+ p2->SetFillColor(0);
+
+ TPad* p3=new TPad("p3","",2./3,0.,1.,1.);
+ p3->SetFillColor(0);
+ p3->SetLogy();
+
+ p1->cd();
+ hevent->SetLineColor(1);
+ hevent->Draw();
+
+ p2->cd();
+ hmult->SetLineColor(1);
+ int j,k,xmax=0;
+ for(j=1;j<=hmult->GetNbinsX();j++) if(hmult->GetBinContent(j)>0.) xmax=j;
+ xmax=(int) (1.1*xmax);
+ if(xmax>hmult->GetNbinsX()) xmax=hmult->GetNbinsX();
+ hmult->GetXaxis()->SetRange(1,xmax);
+ hmult->Draw();
+
+ p3->cd();
+ hmult->Draw();//same as p2, but log scale on y axis
+
+ c->cd();
+ p1->Draw();
+ p2->Draw();
+ p3->Draw();
+
+ c->SaveAs(Form("%s.pdf(",name_fout));
+
+ delete p1;
+ delete p2;
+ delete p3;
+ delete c;
+
+ fout->cd();
+ hevent->Write();
+ hmult->Write();
+
+ //----- get invariant-mass histograms -----
+
+ //Project the THnSparse histograms onto their x axes.
+ hu=(TH1D*) su->Projection(0,"e");
+ hu->SetName("mass_unlike");
+ hm=(TH1D*) sm->Projection(0,"e");
+ hm->SetName("mass_likeMM");
+ hp=(TH1D*) sp->Projection(0,"e");
+ hp->SetName("mass_likePP");
+
+ double A,UA,B,UB;
+
+ hl=(TH1D*) hm->Clone("mass_like");//total like-charge histogram
+ for(j=1;j<=hl->GetNbinsX();j++){
+ A=hm->GetBinContent(j);
+ B=hp->GetBinContent(j);
+ hl->SetBinContent(j,2*sqrt(A*B));//note: 2sqrt(n{K-K-}*n{K+K+})
+ hl->SetBinError(j,sqrt(A+B));
+ }
+ hl->SetTitle("Like-Charge Combinatorial Background");
+
+ hs=(TH1D*) hu->Clone("mass_signal");
+ for(j=1;j<=hs->GetNbinsX();j++){
+ A=hu->GetBinContent(j); UA=hu->GetBinError(j);
+ B=hl->GetBinContent(j); UB=hl->GetBinError(j);
+ hs->SetBinContent(j,A-B);//subtract the combinatorial background
+ hs->SetBinError(j,sqrt(UA*UA+UB*UB));
+ }
+
+ double Imin=1.01,Imax=1.03,dx=hs->GetXaxis()->GetBinWidth(1);
+
+ hb=(TH1D*) hs->Clone("mass_background");//temporary histogram with the peak removed; used to fit the residual background
+ for(j=1;j<=hb->GetNbinsX();j++){
+ A=hb->GetXaxis()->GetBinCenter(j);
+ if(A>Imin && A<Imax){hb->SetBinContent(j,0); hb->SetBinError(j,0);}
+ }
+
+ hs_plot=(TH1D*) hs->Clone("mass_signal_plot");//temporary histogram for plotting
+
+ //----- plot the invariant-mass histograms -----
+
+ c=new TCanvas("c","",10,10,1000,500);
+ c->SetFillColor(0);
+
+ p1=new TPad("p11","",0.,0.,0.5,1.);
+ p1->SetFillColor(0);
+ p1->SetRightMargin(0.05);
+
+ p2=new TPad("p21","",0.5,0.,1.,1.);
+ p2->SetFillColor(0);
+ p2->SetRightMargin(0.05);
+
+ p1->cd();
+ hu->SetLineColor(1); hu->SetMarkerColor(1);
+ hl->SetLineColor(2); hl->SetMarkerColor(2);
+ hu->SetTitle("Invariant Mass (KK) ");
+ hu->SetXTitle("Invariant Mass (KK) [GeV/#it{c}^{2}]");
+ hu->Draw();
+ hl->Draw("same");
+
+ TLatex* t1=new TLatex(0.7,0.96,"Unlike-Charge");
+ t1->SetTextSize(0.04);
+ t1->SetNDC();
+ t1->Draw();
+
+ TLatex* t2=new TLatex(0.7,0.91,"Like-Charge");
+ t2->SetTextSize(0.04);
+ t2->SetNDC();
+ t2->SetTextColor(2);
+ t2->Draw();
+
+ p2->cd();
+ hs->SetLineColor(1); hs->SetMarkerColor(1);
+ hs->SetTitle("Combinatorial Background Subtracted");
+ hs->SetXTitle("Invariant Mass (KK) [GeV/#it{c}^{2}]");
+ hs->Draw();
+
+ c->cd();
+ p1->Draw();
+ p2->Draw();
+
+ c->SaveAs(Form("%s.pdf",name_fout));
+
+ delete p1;
+ delete p2;
+ delete c;
+
+ fout->cd();
+ hu->Write();
+ hl->Write();
+ hs->Write();
+
+ //----- fit the signal histogram -----
+
+ c=new TCanvas("c","",10,10,500,500);
+ c->SetFillColor(0);
+
+ //fit peak and extract yield, mass, and width
+ //The peak is fit using a polynomial for the residual background, plus a Voigtian peak. A Voigtian peak is the convolution of a (non-relativistic) Breit-Wigner peak and a Gaussian (to describe detector effects). The resolution (Gaussian sigma) of the Voigtian peak will be fixed.
+ //Ideally, we would only need to fit once. However, the fits are not always stable so we try a variety of fits and record all results. A similar procedure was followed in the analysis of phi mesons in Pb+Pb collisions (2010 data).
+
+ int jBF;//index controlling the order of the residual background polynomial
+ int nBF=3;//number of different orders used for the residual background polynomial; currently 3 (linear, quadratic, cubic)
+ int jFR;//index controlling the fit region
+ int nFR=4;//number of different fit regions
+ double urb;
+
+ //This TH2D is just an array of numbers to hold the different fit results and associated quantities. The x-axis gives the name of the fit (order of residual background and number of the fit region), the y-axis gives the name of the quantity stored.
+ TH2D* r=new TH2D("fit_results","",nBF*nFR,0,nBF*nFR, 16,0,16);
+ r->GetYaxis()->SetBinLabel(1,"peak integral");
+ r->GetYaxis()->SetBinLabel(2,"mass");
+ r->GetYaxis()->SetBinLabel(3,"width");
+ r->GetYaxis()->SetBinLabel(4,"resolution (fixed)");
+ r->GetYaxis()->SetBinLabel(5,"p0");
+ r->GetYaxis()->SetBinLabel(6,"p1");
+ r->GetYaxis()->SetBinLabel(7,"p2");
+ r->GetYaxis()->SetBinLabel(8,"p3");
+ r->GetYaxis()->SetBinLabel(9,"chi2");
+ r->GetYaxis()->SetBinLabel(10,"NDF");
+ r->GetYaxis()->SetBinLabel(11,"yield (fit)");
+ r->GetYaxis()->SetBinLabel(12,"yield/event (fit)");
+ r->GetYaxis()->SetBinLabel(13,"peak correction factor");
+ r->GetYaxis()->SetBinLabel(14,"yield (bin counting)");
+ r->GetYaxis()->SetBinLabel(15,"yield/event (bin conunting)");
+ r->GetYaxis()->SetBinLabel(16,"yield/event (bin conunting + tails)");
+
+ TF1 *g[3][4],*gp[3][4],*gb[3][4];
+ TFitResultPtr fr;
+ int status;
+ double pc;
+
+ cerr<<"Please wait a moment while the peak fits are performed."<<endl;
+
+ for(jBF=0;jBF<nBF;jBF++) for(jFR=0;jFR<nFR;jFR++){//loops: change the order of the residual background polynomial and the fit region
+ if(!jFR){A=0.995; B=1.06;}//fr0
+ else if(jFR==1){A=1.; B=1.06;}//fr1
+ else if(jFR==2){A=0.995; B=1.07;}//fr2
+ else if(jFR==3){A=1.; B=1.07;}//fr3
+
+ if(!jBF){//pol1
+ g[jBF][jFR]=new TF1(Form("fit_pol1_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])+[4]+[5]*x",A,B);
+ gb[jBF][jFR]=new TF1(Form("fit_back_pol1_fr%i",jFR),"pol1",A,B);
+ gp[jBF][jFR]=new TF1(Form("fit_peak_pol1_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])",A,B);
+ }else if(jBF==1){//pol2
+ g[jBF][jFR]=new TF1(Form("fit_pol2_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])+[4]+[5]*x+[6]*x*x",A,B);
+ gb[jBF][jFR]=new TF1(Form("fit_back_pol2_fr%i",jFR),"pol2",A,B);
+ gp[jBF][jFR]=new TF1(Form("fit_peak_pol2_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])",A,B);
+ }else if(jBF==2){//pol3
+ g[jBF][jFR]=new TF1(Form("fit_pol3_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])+[4]+[5]*x+[6]*x*x+[7]*x*x*x",A,B);
+ gb[jBF][jFR]=new TF1(Form("fit_back_pol3_fr%i",jFR),"pol3",A,B);
+ gp[jBF][jFR]=new TF1(Form("fit_peak_pol3_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])",A,B);
+ }
+
+ for(j=0;j<100;j++){//initial fit of residual background
+ fr=hb->Fit(gb[jBF][jFR],"RSQ");
+ status=fr->Status();
+ if(!status) break;
+ }
+ urb=gb[jBF][jFR]->IntegralError(Imin,Imax)/dx;//estimated uncertainty of residual background integral
+
+ j=hs->GetXaxis()->FindBin(Imin+1e-5); k=hs->GetXaxis()->FindBin(Imax-1e-5);
+ g[jBF][jFR]->SetParameter(0,hs->Integral(j,k)*dx);//initial value of peak integral
+ g[jBF][jFR]->SetParameter(1,1.019455); g[jBF][jFR]->FixParameter(1,1.019455);//fix mass to vacuum value
+ g[jBF][jFR]->SetParameter(2,0.00426); g[jBF][jFR]->FixParameter(2,0.00426);//fix width to vacuum value
+ g[jBF][jFR]->SetParameter(3,0.0011); g[jBF][jFR]->FixParameter(3,0.0011);//fix resolution to 1.1 MeV/c^2 (will not be changed)
+ for(j=0;j<=jBF+1;j++){//fix residual background to the initial fit
+ A=gb[jBF][jFR]->GetParameter(j);
+ g[jBF][jFR]->SetParameter(j+4,A);
+ g[jBF][jFR]->SetParError(j+4,gb[jBF][jFR]->GetParError(j));
+ g[jBF][jFR]->FixParameter(j+4,A);
+ }
+
+ cerr<<"Fitting pol"<<jBF+1<<"_fr"<<jFR<<endl;
+
+ for(j=0;j<100;j++){//fit with only the peak integral free
+ fr=hs->Fit(g[jBF][jFR],"RSQ");
+ status=fr->Status();
+ if(!status) break;
+ }
+
+ g[jBF][jFR]->ReleaseParameter(2);//release width parameter
+ for(j=0;j<100;j++){
+ fr=hs->Fit(g[jBF][jFR],"RSQ");
+ status=fr->Status();
+ if(!status) break;
+ }
+
+ g[jBF][jFR]->ReleaseParameter(1);//release mass parameter
+ for(j=0;j<100;j++){
+ fr=hs->Fit(g[jBF][jFR],"RSQ");
+ status=fr->Status();
+ if(!status) break;
+ }
+
+ g[jBF][jFR]->ReleaseParameter(4);//release residual background constant parameter
+ for(j=0;j<100;j++){
+ fr=hs->Fit(g[jBF][jFR],"RSQ");
+ status=fr->Status();
+ if(!status) break;
+ }
+
+ for(j=1;j<=gb[jBF][jFR]->GetNpar();j++) g[jBF][jFR]->ReleaseParameter(4+j);//release other residual background parameters
+ for(j=0;j<100;j++){
+ fr=hs->Fit(g[jBF][jFR],"RSQ");
+ status=fr->Status();
+ if(!status) break;
+ }
+
+ for(j=0;j<100;j++){//redo fit using "I" option, typically does not affect the result very much, but best to do it anyway. Comment this out if you are debugging and want to save time.
+ fr=hs->Fit(g[jBF][jFR],"RSQI");
+ status=fr->Status();
+ if(!status) break;
+ }
+
+ fout->cd();
+ g[jBF][jFR]->Write();//save the combined fit
+
+ for(j=0;j<gb[jBF][jFR]->GetNpar();j++) gb[jBF][jFR]->SetParameter(j,g[jBF][jFR]->GetParameter(j+4));//put the background parameters into the function gb
+ for(j=0;j<4;j++) gp[jBF][jFR]->SetParameter(j,g[jBF][jFR]->GetParameter(j));//put the peak parameters into the funciton gp
+
+ j=nFR*jBF+jFR+1;//bin for the x-axis of results histogram
+ r->GetXaxis()->SetBinLabel(j,Form("pol%i_fr%i",jBF+1,jFR));
+ for(k=0;k<g[jBF][jFR]->GetNpar();k++){//store the values of the peak parameters
+ r->SetBinContent(j,k+1,g[jBF][jFR]->GetParameter(k));
+ r->SetBinError(j,k+1,g[jBF][jFR]->GetParError(k));
+ }
+ r->SetBinContent(j,9,g[jBF][jFR]->GetChisquare());//store chi^2
+ r->SetBinContent(j,10,g[jBF][jFR]->GetNDF());//store number of degrees of freedom
+
+ A=r->GetBinContent(j,1)-gp[jBF][jFR]->Integral(0.,2*0.493667);//subtract integral of peak below kinematic cutoff
+ UA=A*r->GetBinError(j,1)/r->GetBinContent(j,1);
+ r->SetBinContent(j,11,A/dx); r->SetBinError(j,11,UA/dx);//peak yield extracted from the fit
+ r->SetBinContent(j,12,A/dx/nevt/dy); r->SetBinError(j,12,UA/dx/nevt/dy);//peak yield per event extracted from the fit, corrected for dy
+
+ B=gp[jBF][jFR]->Integral(Imin,Imax);
+ pc=B/A;//peak correction factor: the yield in the interval (Imin,Imax) divided by the total integral above the kinematic cutoff, used to correct the yield from bin counting to account for the yield in the tails
+ r->SetBinContent(j,13,pc);//store peak correction factor
+
+ A=UA=0;
+ for(k=hs->GetXaxis()->FindBin(Imin+1e-5);k<=hs->GetXaxis()->FindBin(Imax-1e-5);k++){//yield from bin counting
+ A+=hs->GetBinContent(k);
+ UA+=pow(hs->GetBinError(k),2);
+ }
+ A-=gb[jBF][jFR]->Integral(Imin,Imax)/dx;//subtract residual background integral
+ UA+=urb*urb;
+
+ r->SetBinContent(j,14,A); r->SetBinError(j,14,sqrt(UA));//yield from bin counting
+ r->SetBinContent(j,15,A/nevt/dy); r->SetBinError(j,15,sqrt(UA)/nevt/dy);//yield per event from bin counting, corrected for dy
+ r->SetBinContent(j,16,A/nevt/pc/dy); r->SetBinError(j,16,sqrt(UA)/nevt/pc/dy);//yield per event from bin counting, corrected for dy, corrected to account for yield in tails
+ }
+
+ fout->cd();
+ r->Write();//save the results histogram in case it is needed later
+
+ delete c;
+
+ //----- plot the fits -----
+
+ c=new TCanvas("c","",10,10,1000,1000);
+ c->SetFillColor(0);
+
+ p1=new TPad("p1_2","",0,0.5,0.5,1);
+ p1->SetFillColor(0);
+ p1->SetRightMargin(0.05);
+
+ p2=new TPad("p2_2","",0.5,0.5,1,1);
+ p2->SetFillColor(0);
+ p2->SetRightMargin(0.05);
+
+ p3=new TPad("p3_2","",0,0,0.5,0.5);
+ p3->SetFillColor(0);
+ p3->SetRightMargin(0.05);
+
+ TPad* p4=new TPad("p4_2","",0.5,0,1,0.5);
+ p4->SetFillColor(0);
+
+ hs_plot->GetXaxis()->SetRangeUser(0.995,1.07-1e-5);
+ hs_plot->SetLineColor(1); hs_plot->SetMarkerColor(1);
+ hs_plot->SetTitle(hs->GetTitle());
+ hs_plot->SetXTitle(hs->GetXaxis()->GetTitle());
+
+ for(jBF=0;jBF<nBF;jBF++) for(jFR=0;jFR<nFR;jFR++){
+ //assign styles and colors to the functions
+ g[jBF][jFR]->SetNpx(300);
+ if(!jFR) g[jBF][jFR]->SetLineColor(2);
+ else if(jFR==1) g[jBF][jFR]->SetLineColor(TColor::GetColor("#ffaa00"));
+ else if(jFR==2) g[jBF][jFR]->SetLineColor(TColor::GetColor("#238e23"));
+ else if(jFR==3) g[jBF][jFR]->SetLineColor(4);
+ gb[jBF][jFR]->SetLineColor(g[jBF][jFR]->GetLineColor());
+ gb[jBF][jFR]->SetLineStyle(2);
+ }
+
+ p1->cd();
+ hs_plot->Draw();
+ for(jFR=0;jFR<nFR;jFR++) gb[0][jFR]->Draw("same");
+ for(jFR=0;jFR<nFR;jFR++) g[0][jFR]->Draw("same");
+ t1=new TLatex(0.94,0.87,"linear RB (pol1)"); t1->SetTextAlign(32);
+ t1->SetNDC();
+ t1->Draw();
+
+ p2->cd();
+ hs_plot->Draw();
+ for(jFR=0;jFR<nFR;jFR++) gb[1][jFR]->Draw("same");
+ for(jFR=0;jFR<nFR;jFR++) g[1][jFR]->Draw("same");
+ t2=new TLatex(0.94,0.87,"quadratic RB (pol2)"); t2->SetTextAlign(32);
+ t2->SetNDC();
+ t2->Draw();
+
+ p3->cd();
+ hs_plot->Draw();
+ for(jFR=0;jFR<nFR;jFR++) gb[2][jFR]->Draw("same");
+ for(jFR=0;jFR<nFR;jFR++) g[2][jFR]->Draw("same");
+ TLatex* t3=new TLatex(0.94,0.87,"cubic RB (pol3)"); t3->SetTextAlign(32);
+ t3->SetNDC();
+ t3->Draw();
+
+ //additional information in pad p4
+
+ p4->cd();
+ TLine* dummy1=new TLine(0,0,1,0);
+ dummy1->SetLineColor(1); dummy1->SetLineWidth(2);
+ TLine* dummy2=new TLine(0,0,1,0);
+ dummy2->SetLineColor(1); dummy2->SetLineWidth(2); dummy2->SetLineStyle(gb[0][0]->GetLineStyle());
+ TLegend* L=new TLegend(0.3,0.69,0.7,0.99);
+ L->SetFillColor(0);
+ L->AddEntry(g[0][0],"Fit Region 0 (0.995,1.06)","l");
+ L->AddEntry(g[0][1],"Fit Region 1 (1,1.06)","l");
+ L->AddEntry(g[0][2],"Fit Region 2 (0.995,1.07)","l");
+ L->AddEntry(g[0][3],"Fit Region 3 (1,1.07)","l");
+ L->AddEntry(dummy1,"Combined Fit","l");
+ L->AddEntry(dummy2,"Residual Background Fit","l");
+ L->Draw();
+
+ //----- check results -----
+
+ double ex_yield=0;
+ if(!strcmp(system,"pp276")) ex_yield=1.416e-3;
+ else if(!strcmp(system,"pp7")) ex_yield=1.706e-3;
+ else if(!strcmp(system,"PbPb276")) ex_yield=0.245;
+ else{
+ cerr<<"Warning in QAphi(): Unknown collision system "<<system<<". Expected yield not known."<<endl;
+ ex_yield=1e10;
+ }
+ double yield=r->GetBinContent(8,16);
+ int status_yield;
+ if(yield/ex_yield<1.5 && yield/ex_yield>0.5) status_yield=0;//OK
+ else if(yield/ex_yield<5 && yield/ex_yield>0.2) status_yield=1;//problem
+ else status_yield=2;//big problem
+
+ double ex_mass=1.019455;
+ double mass=r->GetBinContent(8,2);
+ int status_mass;
+ if(fabs(mass-ex_mass)<1) status_mass=0;//OK
+ else if(fabs(mass-ex_mass)<3) status_mass=1;//problem
+ else status_mass=2;//big problem
+
+ double ex_width=0.00426;
+ double width=r->GetBinContent(8,3);
+ int status_width;
+ if(fabs(width-ex_width)<1) status_width=0;//OK
+ else if(fabs(width-ex_width)<2) status_width=1;//problem
+ else status_width=2;//big problem;
+
+ double tx=0.005,ty=0.64,dty=0.06;
+
+ TString s;
+
+ j=su->GetAxis(1)->GetFirst(); k=su->GetAxis(1)->GetLast();
+ s.Form("%1.2f < #it{p}_{T} < %1.2f GeV/#it{c}, ",su->GetAxis(1)->GetBinLowEdge(j),su->GetAxis(1)->GetBinLowEdge(k+1));
+ j=su->GetAxis(2)->GetFirst(); k=su->GetAxis(2)->GetLast();
+ s.Append(Form("%1.2f < #it{y} < %1.2f",su->GetAxis(2)->GetBinLowEdge(j),su->GetAxis(2)->GetBinLowEdge(k+1)));
+ if(!ptOK || !yOK) s.Append(" [PROBLEM]");
+ TLatex* a1=new TLatex(tx,ty,s.Data());
+ SetText(a1,ptOK);
+ a1->Draw();
+
+ TLatex* a2=new TLatex(tx,ty-1*dty,Form("#phi Yield/event = %1.5e #pm %1.5e (%s)",r->GetBinContent(8,16),r->GetBinError(8,16),r->GetXaxis()->GetBinLabel(8)));
+ SetText(a2,status_yield);
+ a2->Draw();
+
+ s.Form("Expected Yield = %1.5e",ex_yield);
+ if(status_yield) s.Append(" [PROBLEM]");
+ TLatex* a3=new TLatex(tx+0.1,ty-2*dty,s.Data());
+ SetText(a3,status_yield);
+ a3->Draw();
+
+ TLatex* a4=new TLatex(tx,ty-3*dty,Form("#phi Mass = %1.5f #pm %1.5f GeV/#it{c}^{2}",r->GetBinContent(8,2),r->GetBinError(8,2)));
+ SetText(a4,status_mass);
+ a4->Draw();
+
+ s.Form("PDG Mass = 1.019455");
+ if(status_mass) s.Append(" [PROBLEM]");
+ TLatex* a5=new TLatex(tx+0.1,ty-4*dty,s.Data());
+ SetText(a5,status_mass);
+ a5->Draw();
+
+ TLatex* a6=new TLatex(tx,ty-5*dty,Form("#phi Width = %1.5f #pm %1.5f MeV/#it{c}^{2}",r->GetBinContent(8,3)*1000.,r->GetBinError(8,3)*1000.));
+ SetText(a6,status_width);
+ a6->Draw();
+
+ s.Form("PDG Width = 4.26");
+ if(status_width) s.Append(" [PROBLEM]");
+ TLatex* a7=new TLatex(tx+0.1,ty-6*dty,s.Data());
+ SetText(a7,status_width);
+ a7->Draw();
+
+ if(!ptOK || !yOK || status_yield || status_mass || status_width){
+ status=2;
+ s.Form("THERE ARE PROBLEMS!!!");
+ }else{
+ status=0;
+ s.Form("Things look OK.");
+ }
+
+ TLatex* a8=new TLatex(tx,ty-7*dty,s.Data());
+ SetText(a8,status); a8->SetTextSize(0.05);
+ a8->Draw();
+
+ if(status) s.Form("PLEASE CHECK THE ERROR MESSAGES.");
+ else s.Form("");
+
+ TLatex* a9=new TLatex(tx,ty-8*dty,s.Data());
+ SetText(a9,status); a9->SetTextSize(0.05);
+ a9->Draw();
+
+ c->cd();
+ p1->Draw();
+ p2->Draw();
+ p3->Draw();
+ p4->Draw();
+
+ c->SaveAs(Form("%s.pdf)",name_fout));
+
+ //----- print out results -----
+
+ j=su->GetAxis(1)->GetFirst(); k=su->GetAxis(1)->GetLast();
+ printf("%1.2f < pT < %1.2f GeV/c, ",su->GetAxis(1)->GetBinLowEdge(j),su->GetAxis(1)->GetBinLowEdge(k+1));
+ j=su->GetAxis(2)->GetFirst(); k=su->GetAxis(2)->GetLast();
+ printf("%1.2f < y < %1.2f\n",su->GetAxis(2)->GetBinLowEdge(j),su->GetAxis(2)->GetBinLowEdge(k+1));
+ printf("phi Yield/event = %1.5e +/- %1.5e (%s)\n",r->GetBinContent(8,16),r->GetBinError(8,16),r->GetXaxis()->GetBinLabel(8));
+ if(status_yield) printf("*** PROBLEM: phi yield too far from expected value (%1.5e).\n",ex_yield);
+ printf("phi Mass = %1.5f +/- %1.5f GeV/c^2 (PDG Value = 1.019455)\n",r->GetBinContent(8,2),r->GetBinError(8,2));
+ if(status_mass) printf("*** PROBLEM: phi mass too far from expected value.\n");
+ printf("phi Width = %1.5f +/- %1.5f MeV/c^2 (PDG Value = 4.26)\n",r->GetBinContent(8,3)*1000.,r->GetBinError(8,3)*1000.);
+ if(status_width) printf("*** PROBLEM: phi width too far from expected value.\n");
+ if(status) printf("\n*** One or more parameters does not have the expected value. Please check the error messages and read the instructions at the bottom of the macro file.\n");
+
+ //----- finish -----
+
+ cerr<<"\nMacro has finished successfully."<<endl;
+
+ fin->Close();
+ fout->Close();
+
+ return;
+}
+
+
+bool SetPtRange(THnSparse* h){
+ TAxis* a=h->GetAxis(1);
+ double min=0.5,max=1.5;
+
+ int b1=a->FindBin(1.00001*min);
+ int b2=a->FindBin(0.99999*max);
+ a->SetRange(b1,b2);
+ b1=a->GetFirst();
+ b2=a->GetLast();
+
+ //min and/or max are not bin boundaries
+ if(fabs(a->GetBinLowEdge(b1)-min)>1e-5*min || fabs(a->GetBinLowEdge(b2+1)-max)>1e-5*max){
+ cerr<<"Error in QAphi(): pT range cannot be set to ("<<min<<","<<max<<"). The yield may therefore be different from the expected value."<<endl;
+ return false;
+ }
+
+ return true;
+}
+
+bool SetRapidityRange(THnSparse* h,double& dy){
+ TAxis* a=h->GetAxis(2);
+ double min=-0.5,max=0.5;
+
+ int b1=a->FindBin(min+1e-5);
+ int b2=a->FindBin(max-1e-5);
+ a->SetRange(b1,b2);
+ b1=a->GetFirst();
+ b2=a->GetLast();
+ dy=a->GetBinLowEdge(b2+1)-a->GetBinLowEdge(b1);
+
+ //min and/or max are not bin boundaries
+ if(fabs(a->GetBinLowEdge(b1)-min)>1e-5 || fabs(a->GetBinLowEdge(b2+1)-max)>1e-5){
+ cerr<<"Error in QAphi(): rapidity range cannot be set to |y|<0.5. Although this macro does correct for the rapidity range dy, using a different rapidity range may cause the yield to be different from the expected value."<<endl;
+ return false;
+ }
+
+ return true;
+}
+
+
+void SetText(TLatex* t,int flag){
+ t->SetTextSize(0.04);
+ if(!flag) t->SetTextColor(TColor::GetColor("#238e23"));
+ else if(flag==1) t->SetTextColor(TColor::GetColor("#cc5500"));
+ else t->SetTextColor(2);
+ return;
+}
+
+
+void SetText(TLatex* t,bool flag){
+ if(flag) SetText(t,0);
+ else SetText(t,2);
+ return;
+}
+
+
+/*-----------------------------------------------
+
+INSTRUCTIONS:
+
+*Introduction:
+
+This macro is designed to read the output of the analysis task to check that the yield per event, mass, and width of the phi meson in the production are within acceptable ranges. The analysis task constructs invariant-mass distributions of charged-kaon pairs in the vicinity of the mass of the phi meson (1.019455 GeV/c^2). The branching ratio of the phi->K-K+ decay is 0.489. In this macro, a combinatorial background is constructed using the like-charge distributions (K-K- and K+K+) and subtracted from the unlike-charge (K-K+) distribution. This leaves a phi peak sitting on top of a residual background. The background-subtracted distribution is fit with a Voigtian peak (a convolution of a non-relativistic Breit-Wigner peak with a Gaussian to account for detector effects) plus a polynomial to describe the residual background. Because these fits sometimes fail, a total of 12 fits are performed. There are 3 choices of residual background polynomial (linear, quadratic, or cubic) and 4 different fit regions (defined in the macro). For more information on the fitting procedure, see refs. [2] and [3]. The yield, mass, and width of the phi meson extracted from the default fit (quadratic residual background, fit region 3) are compared to the expected values.
+
+*Output:
+
+This macro produces an output PDF file and an output ROOT file.
+
+**Output PDF:
+
+The output PDF consists of three canvases, organized into panels as follows:
+
+Canvas 1: [[1a][1b][1c]]
+Canvas 2: [[2a][2b]]
+Canvas 3: [[3a][3b]]
+ [[3c][3d]]
+
+The contents of the panels is:
+1a.) histogram hEventStat (see description below)
+1b.) histogram hAEventsVsMulti (see description below)
+1c.) histogram hAEventsVsMulti, logarithmic y-axis
+2a.) unlike-charge invariant-mass distribution and like-charge combinatorial background
+2b.) background-subtracted invariant-mass distribution
+3a.) background-subtracted invariant-mass distribution with fits. The fits shown here assume a linear residual background. A fit is shown for each of the four different fit regions. Fit region 3 (blue) is plotted on top. The dashed lines are the residual background.
+3b.) same as 3a, but with a quadratic residual background. Note that the blue fit in this panel is the default fit, from which the yield, mass, and width are extracted.
+3c.) same as 3a, but with a cubic residual background
+3d.) a legend describing the fit functions and a printout of information. The printed information should be green. If it is not, there is a problem.
+
+**Output ROOT file:
+
+The output ROOT file contains the following:
+1.) histogram hEventStat (see description below)
+2.) histogram hAEventsVsMulti (see description below)
+3.) histogram mass_unlike: the unlike-charge invariant-mass distribution
+4.) histogram mass_like: the like-charge combinatorial background
+5.) histogram mass_signal: background-subtracted invariant-mass distribution (mass_unlike - mass_like)
+6.) functions fit_pol*_fr*: the fit functions for mass_signal. pol1: linear residual background, pol2: quadratic residual background, pol3: cubic residual background. fr* indicates the fit region.
+7.) histogram fit_results: contains the results of the different fits (see below)
+
+***Histogram fit_results:
+
+This is an array of numbers (with uncertainties) stored in a TH2D. The x-axis is the name of the fit (e.g., "pol2_fr3") and the y-axis gives the name of the quantity stored. The quantities stored (numbered along the y-axis) are:
+
+1.) peak integral: integral of the phi peak in the fit function (parameter [0])
+2.) peak mass: parameter [1]
+3.) peak width: parameter [2]
+4.) peak resolution: parameter [3]: currently fixed to 1.1 MeV/c^2, but stored anyway
+5-8.) coefficients of the residual background polynomial (parameters [4] through [7])
+9.) fit chi^2
+10.) fit number of degrees of freedom
+11.) yield extracted from the fit: parameter [0] minus the yield below the kinematic cutoff and corrected for the width of the invariant-mass bins
+12.) yield/event (fit): the value above divided by the number of events and the width of the rapidity range (dy)
+13.) peak correction factor: the fraction of the peak that lies inside the range (1.01,1.03) GeV/c^2
+14.) yield (bin counting): the integral of the mass_signal histogram for (1.01,1.03) minus the integral of the residual background over the same interval
+15.) yield/event (bin counting) the value above divided by the number of events and the width of the rapidity range (dy)
+16.) yield/event (bin counting + tails) the value above further corrected by dividing by the peak correction factor (now accounting for the yield in the tails)
+
+The default fit is "pol2_fr3". Quantity 2 is reported as the mass, quantity 3 is reported as the width, and quantity 16 is reported as the yield. So, to find the default mass (which is printed out) do fit_results->GetBinContent(8,2). The default width is fit_results->GetBinContent(8,3) and the default yield is fit_results->GetBinContent(8,16).
+
+*More Information:
+
+In the analysis task, the kaons should be selected using track selection cuts and a 3-sigma cut (about the kaon mean) on TPC dE/dx. See the analysis task for the exact cuts.
+
+The output of the analysis task is expected to contain:
+1.) the histogram hEventStat, which gives the number of events that pass different selection criteria
+2.) the histogram hAEventsVsMulti, which gives the number of events in each multiplicity (for pp) or centrality (for Pb-Pb) bin
+3.) three THnSparse histograms with three dimensions and the suffixes "Unlike", "LikePP", and "LikeMM"
+
+For each THnSparse, axis 0 is the KK invariant mass, axis 1 is the transverse momentum, and axis 2 is the rapidity. This macro attempts to select the range 0.5<pT<1.5 GeV/c; if this is not possible an error message will be generated, as the expected yields were calculated assuming that pT range. This macro uses the full rapidity range on axis 2. The expected yields were calculated assuming a rapidity range of |y|<0.5. The macro corrects for the width of the rapidity range (dy), but using a different rapidity range may cause the yield to deviate from the expected value.
+
+*Expected Values
+
+**Expected Yields
+
+If the yield is less than 50% or more than 150% of the expected value, it is flagged as an orange problem. If the yield deviates from the expected value by more than a factor of 5, it is flagged as a red problem.
+
+The expected yields differ depending on the collision system.
+
+pp Collisions at 7 TeV (system="pp7"): The expected yield of 1.706e-3 is computed based on the published (ref. [1]) phi spectrum. The yield per event for 0.5<pT<1.5 GeV/c is multiplied by the efficiency and the branching ratio.
+
+pp Collisions at 2.76 TeV (system="pp276"): The expected yield of 1.416e-3 is computed by extrapolating from the pp 7 TeV yield, assuming that the phi yield scales as s^0.1 (energy dependence taken from ref. [2]).
+
+Pb-Pb Collisions at 2.76 TeV (system="PbPb276"): The expected yield of 0.245 is computed based on the nearly published (ref. [2]) phi spectrum.
+
+**Expected Mass
+
+The expected phi mass is the PDG value: 1.019455 GeV/c^2. A deviation from this value of more than 1 MeV/c^2 is flagged as an orange problem. A deviation from this value of more than 3 MeV/c^2 is flagged as a red problem.
+
+**Expected Width
+
+The expected phi width is the PDG value: 4.26 MeV/c^2. A deviation from this value of more than 1 MeV/c^2 is flagged as an orange problem. A deviation from this value of more than 3 MeV/c^2 is flagged as a red problem.
+
+*Troubleshooting
+
+If a parameter does not have the expected value, it will be flagged. In the output PDF file, a problematic parameter will be colored orange or red, with red indicating a larger deviation from the expected value.
+
+If the yield, mass, or width deviates from the expected value, you should check to see if the default fit is bad. Look at the blue fit in the upper right plot of the third canvas in the PDF file (panel 3b). This is the default fit. Does it describe the data? Does the residual background (dashed line) appear to be reasonable. If the default fit looks bad, look through the other fits and see if you can find a fit that looks OK. You can read out the yield, mass, and width from the fit_results histogram. Compare these values to the expected values. It might also make sense to compare the results of all fits to the expected value. Are the fit parameters always outside the acceptable range?
+
+If the phi yield deviates from the expected value, the following should be noted. The yield depends on
+1.) collision system and energy
+2.) the triggers used
+3.) multiplicity or centrality range
+4.) cuts used to select the decay daughters
+5.) pT and rapidity range for phi mesons
+
+The expected yields were calculated for
+1.) pp and Pb-Pb collisions at 2.76 TeV, pp collisions at 7 TeV
+2.) minimum-bias triggers
+3.) all multiplicities for pp, centrality 0-90% for Pb-Pb
+4.) the track selection cuts used in the original version of the analysis task
+5.) 0.5<pT<1.5 GeV/c and |y|<0.5 for phi mesons
+
+Changes to any of these criteria can cause the yield to change. If the criteria you are using do not exactly match the criteria for the expected yield, you will need to correct for those differences. The macro corrects for the rapidity range, so SMALL changes to the rapidity window will be accounted for to some extent.
+
+If the width deviates from the expected value, it is possible that the resolution is not 1.1 MeV/c^2 (the value to which it is currently fixed). Studying this issue is beyond the scope of these instructions, but you should be aware of the possibility.
+
+*References
+
+[1]: B. Abelev et al. (ALICE Collaboration), Eur. Phys. J. C 72, 2183 (2012)
+[2]: Paper in preparation: "K*(892)^0 and phi(1020) resonances in Pb-Pb collisions at 2.76 TeV", by the ALICE Collaboration, intended for publication in Phys. Rev. C (2014)
+[3]: Analysis Note: A. G. Knospe, "Yield of phi mesons at low pT in Pb-Pb collisions at 2.76 TeV (2010 data)", ALICE-ANA-2012-300, https://aliceinfo.cern.ch/Notes/node/42
+
+ -----------------------------------------------*/
--- /dev/null
+//////////////////////////////////////////////////\r
+//\r
+// This macro was written by Domenico Colella (domenico.colella@cern.ch).\r
+// 12 November 2013\r
+//\r
+// ------ Arguments\r
+// -- icasType = 0) Xi- 1) Xi+ 2) Omega- 3) Omega+\r
+// -- collidingsystem = 0) PbPb 1) pp\r
+// -- fileDir = "Input file directory"\r
+// -- filein = "Input file name"\r
+//\r
+// ------ QATask output content\r
+// The output produced by the QATask is a CFContainer with 4 steps and 21 variables.\r
+// The meaning of each variable within the container are listed here:\r
+// -- 0 = Max DCA Cascade Daughters pp: 2.0 PbPb: 0.3\r
+// -- 1 = Min DCA Bach To PV pp: 0.01 PbPb: 0.03\r
+// -- 2 = Min Cascade Cosine Of PA pp: 0.98 PbPb: 0.999\r
+// -- 3 = Min Cascade Radius Fid. Vol. pp: 0.2 PbPb: 0.9\r
+// -- 4 = Window Invariant Mass Lambda pp: 0.008 PbPb: 0.0008\r
+// -- 5 = Max DCA V0 Daughters pp: 1.5 PbPb: 1.0\r
+// -- 6 = Min V0 Cosine Of PA To PV pp: pT dep. PbPb: 0.98\r
+// -- 7 = Min V0 Radius Fid. Vol. pp: 0.2 PbPb: 0.9\r
+// -- 8 = Min DCA V0 To PV pp: 0.01 PbPb: 0.05\r
+// -- 9 = Min DCA Pos To PV pp: 0.05 PbPb: 0.1\r
+// -- 10 = Min DCA Neg To PV pp: 0.05 PbPb: 0.1\r
+// -- 11 = Invariant Mass distribution for Xi\r
+// -- 12 = Invariant Mass distribution for Omega\r
+// -- 13 = Transverse Momentum distribution\r
+// -- 14 = Rapidity distribution for Xi\r
+// -- 15 = Rapidity distribution for Omega\r
+// -- 16 = Proper length distribution for the cascade\r
+// -- 17 = Proper length distribution for the V0\r
+// -- 18 = Min V0 Cosine Of PA To Xi Vertex pp: pT dep. PbPb: pT dep.\r
+// -- 19 = Centrality\r
+// -- 20 = ESD track multiplicity\r
+// The last two variables are empty in the case proton-proton collisions.\r
+//\r
+// ------ Present Macro Check\r
+// With this macro are produced the plots of the distributions for the topological\r
+// variables used during the reconstruction of the cascades and defined in the two\r
+// classes AliCascadeVertexer.cxx and AliV0vertexer.cxx contained in /STEER/ESD/ \r
+// folder in Aliroot.\r
+//\r
+// -- First Canvas: DCA cascade daughters, Bachelor IP to PV, Cascade cosine of PA\r
+// Cascade radius of fiducial volume, Invariant mass Lambda,\r
+// DCA V0 daughters.\r
+// -- Second Canvas: V0 cosine of PA to PV, Min V0 Radius Fid. Vol., Min DCA V0 To PV\r
+// Min DCA Pos To PV, Min DCA Neg To PV, V0 cosine of PA to XiV\r
+//\r
+// In this plots, in correspondence to the min/max cut value adopted in the reconstruction\r
+// a line is drawn, to check if there is same problem in the cuts definition.\r
+//\r
+// Also, other specific distribution for the selected cascades are ploted as: the\r
+// invariant mass distribution, the transverse momentum distribution, the rapidity\r
+// distribution, proper length distribution for the cascade and the v0.\r
+//\r
+// -- Third Canvas: InvMass, Transverse momentum, Cascade proper length, V0 proper length\r
+//\r
+// In the end, only for thr PbPb collision the centrality distribution and the\r
+// track multiplicity distribution are sored.\r
+//\r
+// -- Fourth Canvas: Centrality, track multiplicity\r
+//\r
+//\r
+//////////////////////////////////////////////////////\r
+\r
+\r
+\r
+\r
+class AliCFContainer;\r
+\r
+void multistrangeQA(Int_t icasType = 0, // 0) Xi- 1) Xi+ 2) Omega- 3) Omega+\r
+ Int_t collidingsystem = 0, // 0) PbPb 1) pp\r
+ Char_t *fileDir = "./", // Input file directory\r
+ Char_t *filein = "AnalysisResults_AOD.root" // Input file name\r
+ ) {\r
+\r
+\r
+ /////////////\r
+ gStyle->SetOptStat(1110);\r
+ gStyle->SetOptStat(kFALSE);\r
+ gStyle->SetOptTitle(kFALSE);\r
+ gStyle->SetFrameLineWidth(2.5);\r
+ gStyle->SetCanvasColor(0);\r
+ gStyle->SetPadColor(0);\r
+ gStyle->SetHistLineWidth(2.5);\r
+ gStyle->SetLabelSize(0.05, "x");\r
+ gStyle->SetLabelSize(0.05, "y");\r
+ gStyle->SetTitleSize(0.05, "x");\r
+ gStyle->SetTitleSize(0.05, "y");\r
+ gStyle->SetTitleOffset(1.1, "x");\r
+ gStyle->SetPadBottomMargin(0.14);\r
+ gSystem->Load("libANALYSIS.so");\r
+ gSystem->Load("libANALYSISalice.so");\r
+ gSystem->Load("libCORRFW.so");\r
+\r
+ \r
+\r
+ TFile *f1 = new TFile(Form("%s/%s",fileDir,filein));\r
+ AliCFContainer *cf = (AliCFContainer*) (f1->Get("PWGLFStrangeness.outputCheckCascade/fCFContCascadeCuts"));\r
+ \r
+\r
+ //DEEFINE TEXT\r
+ TLatex* t1 = new TLatex(0.6,0.55,"#color[3]{OK!!}");\r
+ t1->SetTextSize(0.1);\r
+ t1->SetNDC();\r
+ TLatex* t2 = new TLatex(0.6,0.55,"#color[2]{NOT OK!!}");\r
+ t2->SetTextSize(0.1);\r
+ t2->SetNDC();\r
+ t2->SetTextColor(2);\r
+ TLatex* tcasc;\r
+ if (icasType == 0) tcasc = new TLatex(0.8,0.7,"#color[1]{#Xi^{-}}");\r
+ else if (icasType == 1) tcasc = new TLatex(0.8,0.7,"#color[1]{#Xi^{+}}");\r
+ else if (icasType == 2) tcasc = new TLatex(0.8,0.7,"#color[1]{#Omega^{-}}");\r
+ else if (icasType == 3) tcasc = new TLatex(0.8,0.7,"#color[1]{#Omega^{+}}");\r
+ tcasc->SetTextSize(0.2);\r
+ tcasc->SetNDC();\r
+ tcasc->SetTextColor(2);\r
+ TLatex* tpdgmass;\r
+ if (icasType == 0) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.321 GeV/c^{2}}");\r
+ else if (icasType == 1) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.321 GeV/c^{2}}");\r
+ else if (icasType == 2) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.672 GeV/c^{2}}");\r
+ else if (icasType == 3) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.672 GeV/c^{2}}");\r
+ tpdgmass->SetTextSize(0.07);\r
+ tpdgmass->SetNDC();\r
+ tpdgmass->SetTextColor(2);\r
+ \r
+ //DEFINE 1st CANVAS AND DRAW PLOTS\r
+ TCanvas *c1 = new TCanvas("c1","",1200,800);\r
+ c1->Divide(2,3); \r
+ //Pad 1: DCA cascade daughters\r
+ c1->cd(1);\r
+ gPad->SetLogy();\r
+ TH1D *hvar0 = cf->ShowProjection(0,icasType);\r
+ hvar0->Draw("histo");\r
+ Double_t x0;\r
+ if (collidingsystem == 0) x0 = 0.3;\r
+ else if (collidingsystem == 1) x0 = 2.0;\r
+ TLine *line0 = new TLine(x0,0.,x0,hvar0->GetBinContent(hvar0->GetMaximumBin()));\r
+ line0->SetLineColor(kRed);\r
+ line0->SetLineStyle(9);\r
+ line0->SetLineWidth(2.0);\r
+ line0->Draw("same");\r
+ Bool_t check_0 = checkOverTheLimit(hvar0, x0);\r
+ if (check_0) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ tcasc->Draw();\r
+ //Pad 2: Bachelor IP to PV\r
+ c1->cd(2);\r
+ gPad->SetLogy();\r
+ TH1D *hvar1 = cf->ShowProjection(1,icasType);\r
+ hvar1->GetXaxis()->SetRangeUser(0.,0.24);\r
+ hvar1->Draw("histo");\r
+ Double_t x1;\r
+ if (collidingsystem == 0) x1 = 0.03;\r
+ else if (collidingsystem == 1) x1 = 0.01;\r
+ TLine *line1 = new TLine(x1,0.,x1,hvar1->GetBinContent(hvar1->GetMaximumBin()));\r
+ line1->SetLineColor(kRed);\r
+ line1->SetLineStyle(9);\r
+ line1->SetLineWidth(2.0);\r
+ line1->Draw("same");\r
+ Bool_t check_1 = checkUnderTheLimit(hvar1, x1);\r
+ if (check_1) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ //Pad 3: Cascade cosine of Pointing Angle\r
+ c1->cd(3);\r
+ gPad->SetLogy();\r
+ TH1D *hvar2 = cf->ShowProjection(2,icasType);\r
+ Double_t max2 = hvar2->GetBinContent(hvar2->GetMaximumBin());\r
+ hvar2->GetYaxis()->SetRangeUser(0.01,max2*1.5);\r
+ hvar2->Draw("histo");\r
+ Double_t x2;\r
+ if (collidingsystem == 0) x2 = 0.999;\r
+ else if (collidingsystem == 1) x2 = 0.98;\r
+ TLine *line2 = new TLine(x2,0.,x2,hvar2->GetBinContent(hvar2->GetMaximumBin()));\r
+ line2->SetLineColor(kRed);\r
+ line2->SetLineStyle(9);\r
+ line2->SetLineWidth(2.0);\r
+ line2->Draw("same");\r
+ line1->Draw("same");\r
+ Bool_t check_2 = checkUnderTheLimit(hvar2, x2);\r
+ if (check_2) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ //Pad 4: Cascade radius of fiducial volume\r
+ c1->cd(4);\r
+ gPad->SetLogy();\r
+ TH1D *hvar3 = cf->ShowProjection(3,icasType);\r
+ hvar3->GetXaxis()->SetRangeUser(0.,3.8);\r
+ hvar3->Draw("histo");\r
+ Double_t x3;\r
+ if (collidingsystem == 0) x3 = 0.9;\r
+ else if (collidingsystem == 1) x3 = 0.2;\r
+ TLine *line3 = new TLine(x3,0.,x3,hvar3->GetBinContent(hvar3->GetMaximumBin()));\r
+ line3->SetLineColor(kRed);\r
+ line3->SetLineStyle(9);\r
+ line3->SetLineWidth(2.0);\r
+ line3->Draw("same");\r
+ Bool_t check_3 = checkUnderTheLimit(hvar3, x3);\r
+ if (check_3) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ //Pad 5: Invariant mass Lambda\r
+ c1->cd(5);\r
+ TH1D *hvar4 = cf->ShowProjection(4,icasType);\r
+ hvar4->Draw("histo");\r
+ Double_t x41 = 1.116 + 0.008;\r
+ TLine *line41 = new TLine(x41,0.,x41,hvar4->GetBinContent(hvar4->GetMaximumBin()));\r
+ line41->SetLineColor(kRed);\r
+ line41->SetLineStyle(9);\r
+ line41->SetLineWidth(2.0);\r
+ line41->Draw("same");\r
+ Double_t x42 = 1.115 - 0.008;\r
+ TLine *line42 = new TLine(x42,0.,x42,hvar4->GetBinContent(hvar4->GetMaximumBin()));\r
+ line42->SetLineColor(kRed);\r
+ line42->SetLineStyle(9);\r
+ line42->SetLineWidth(2.0);\r
+ line42->Draw("same");\r
+ Bool_t check_4_1 = checkUnderTheLimit(hvar3, x3);\r
+ Bool_t check_4_2 = checkOverTheLimit(hvar0, x0);\r
+ if (check_4_1 && check_4_2) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ //Pad 6: DCA V0 daughters\r
+ c1->cd(6);\r
+ gPad->SetLogy();\r
+ TH1D *hvar5 = cf->ShowProjection(5,icasType);\r
+ hvar5->Draw("histo");\r
+ Double_t x5;\r
+ if (collidingsystem == 0) x5 = 1.0;\r
+ else if (collidingsystem == 1) x5 = 1.5;\r
+ TLine *line5 = new TLine(x5,0.,x5,hvar5->GetBinContent(hvar5->GetMaximumBin()));\r
+ line5->SetLineColor(kRed);\r
+ line5->SetLineStyle(9);\r
+ line5->SetLineWidth(2.0);\r
+ line5->Draw("same");\r
+ Bool_t check_5 = checkOverTheLimit(hvar5, x5);\r
+ if (check_5) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ c1->SaveAs("MultistrangeQA.pdf(");\r
+ \r
+\r
+ //DEFINE 2st CANVAS AND DRAW PLOTS\r
+ TCanvas *c2 = new TCanvas("c2","",1200,800);\r
+ c2->Divide(2,3);\r
+ //Pad 1: V0 cosine of Pointing Angle to PV\r
+ c2->cd(1);\r
+ gPad->SetLogy();\r
+ TH1D *hvar6 = cf->ShowProjection(6,icasType);\r
+ Double_t max6 = hvar6->GetBinContent(hvar6->GetMaximumBin());\r
+ hvar6->GetYaxis()->SetRangeUser(0.01,max6*1.5);\r
+ hvar6->Draw("histo");\r
+ //Pad 2: Min V0 Radius Fid. Vol. \r
+ c2->cd(2);\r
+ gPad->SetLogy();\r
+ TH1D *hvar7 = cf->ShowProjection(7,icasType);\r
+ hvar7->GetXaxis()->SetRangeUser(0.,3.0);\r
+ hvar7->Draw("histo");\r
+ Double_t x7;\r
+ if (collidingsystem == 0) x7 = 0.9;\r
+ else if (collidingsystem == 1) x7 = 0.2;\r
+ TLine *line7 = new TLine(x7,0.,x7,hvar7->GetBinContent(hvar7->GetMaximumBin()));\r
+ line7->SetLineColor(kRed);\r
+ line7->SetLineStyle(9);\r
+ line7->SetLineWidth(2.0);\r
+ line7->Draw("same");\r
+ Bool_t check_7 = checkUnderTheLimit(hvar7, x7);\r
+ if (check_7) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ //Pad3: Min DCA V0 To PV\r
+ c2->cd(3);\r
+ gPad->SetLogy();\r
+ TH1D *hvar8 = cf->ShowProjection(8,icasType);\r
+ hvar8->GetXaxis()->SetRangeUser(0.,0.3);\r
+ hvar8->Draw("histo");\r
+ Double_t x8;\r
+ if (collidingsystem == 0) x8 = 0.05;\r
+ else if (collidingsystem == 1) x8 = 0.01;\r
+ TLine *line8 = new TLine(x8,0.,x8,hvar8->GetBinContent(hvar8->GetMaximumBin()));\r
+ line8->SetLineColor(kRed);\r
+ line8->SetLineStyle(9);\r
+ line8->SetLineWidth(2.0);\r
+ line8->Draw("same");\r
+ Bool_t check_8 = checkUnderTheLimit(hvar8, x8);\r
+ if (check_8) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ //Pad 4: Min DCA Pos To PV\r
+ c2->cd(4);\r
+ gPad->SetLogy();\r
+ TH1D *hvar9 = cf->ShowProjection(9,icasType);\r
+ hvar9->GetXaxis()->SetRangeUser(0.,0.2);\r
+ hvar9->Draw("histo");\r
+ Double_t x9;\r
+ if (collidingsystem == 0) x9 = 0.1;\r
+ else if (collidingsystem == 1) x9 = 0.05;\r
+ TLine *line9 = new TLine(x9,0.,x9,hvar9->GetBinContent(hvar9->GetMaximumBin()));\r
+ line9->SetLineColor(kRed);\r
+ line9->SetLineStyle(9);\r
+ line9->SetLineWidth(2.0);\r
+ line9->Draw("same");\r
+ Bool_t check_9 = checkUnderTheLimit(hvar9, x9);\r
+ if (check_9) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ //Pad 5: Min DCA Neg To PV\r
+ c2->cd(5);\r
+ gPad->SetLogy();\r
+ TH1D *hvar10 = cf->ShowProjection(10,icasType);\r
+ hvar10->GetXaxis()->SetRangeUser(0.,0.2);\r
+ hvar10->Draw("histo");\r
+ Double_t x10;\r
+ if (collidingsystem == 0) x10 = 0.1;\r
+ else if (collidingsystem == 1) x10 = 0.05;\r
+ TLine *line10 = new TLine(x10,0.,x10,hvar10->GetBinContent(hvar10->GetMaximumBin()));\r
+ line10->SetLineColor(kRed);\r
+ line10->SetLineStyle(9);\r
+ line10->SetLineWidth(2.0);\r
+ line10->Draw("same");\r
+ Bool_t check_10 = checkUnderTheLimit(hvar10, x10);\r
+ if (check_10) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }\r
+ else { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }\r
+ //Pad 6: V0 cosine of Pointing Angle to XiV\r
+ c2->cd(6);\r
+ gPad->SetLogy();\r
+ TH1D *hvar20 = cf->ShowProjection(18,icasType);\r
+ Double_t max20 = hvar20->GetBinContent(hvar20->GetMaximumBin());\r
+ hvar20->GetYaxis()->SetRangeUser(0.01,max20*1.5);\r
+ hvar20->Draw("histo");\r
+ c2->SaveAs("MultistrangeQA.pdf");\r
+\r
+ //DEFINE 3st CANVAS AND DRAW PLOTS\r
+ TCanvas *c3 = new TCanvas("c3","",1200,800);\r
+ c3->Divide(2,3);\r
+ //Pad 1: InvMass\r
+ c3->cd(1);\r
+ TH1D *hvar12 = cf->ShowProjection(11+icasType/2,icasType);\r
+ hvar12->Draw("histo");\r
+ tpdgmass->Draw(); \r
+ TLine *linemass;\r
+ if (icasType == 0) linemass = new TLine(1.32171,0.,1.32171,0.5*hvar12->GetBinContent(hvar12->GetMaximumBin()));\r
+ else if (icasType == 1) linemass = new TLine(1.32171,0.,1.32171,0.5*hvar12->GetBinContent(hvar12->GetMaximumBin()));\r
+ else if (icasType == 2) linemass = new TLine(1.67245,0.,1.67245,0.5*hvar12->GetBinContent(hvar12->GetMaximumBin()));\r
+ else if (icasType == 3) linemass = new TLine(1.67245,0.,1.67245,0.5*hvar12->GetBinContent(hvar12->GetMaximumBin()));\r
+ linemass->SetLineColor(kRed);\r
+ linemass->SetLineStyle(1);\r
+ linemass->SetLineWidth(2.0);\r
+ linemass->Draw("same");\r
+ //Pad 2: Transverse momentum\r
+ c3->cd(2);\r
+ TH1D *hvar13 = cf->ShowProjection(13,icasType);\r
+ hvar13->Draw("histo");\r
+ //Pad 3: Y\r
+ c3->cd(3);\r
+ TH1D *hvar14 = cf->ShowProjection(14+icasType/2,icasType);\r
+ hvar14->Draw("histo");\r
+ //Pad 4: Cascade proper length\r
+ c3->cd(4);\r
+ TH1D *hvar18;\r
+ hvar18 = cf->ShowProjection(16,icasType);\r
+ hvar18->GetXaxis()->SetRangeUser(0.,90.);\r
+ hvar18->Draw("histo");\r
+ //Pad 5: V0 proper length \r
+ c3->cd(5);\r
+ TH1D *hvar19;\r
+ hvar19 = cf->ShowProjection(17,icasType);\r
+ hvar19->GetXaxis()->SetRangeUser(0.,90.);\r
+ hvar19->Draw("histo");\r
+ //Pad 6\r
+ // empty \r
+ if (collidingsystem == 1) c3->SaveAs("MultistrangeQA.pdf)");\r
+ else if (collidingsystem == 0) c3->SaveAs("MultistrangeQA.pdf");\r
+\r
+ \r
+ //DEFINE 4st CANVAS AND DRAW PLOTS\r
+ TCanvas *c4 = new TCanvas("c4","",600,400);\r
+ c4->Divide(2,1);\r
+ //Pad1: invariant mass fit\r
+ c4->cd(1);\r
+ TH1D *hvar18 = cf->ShowProjection(11+icasType/2,icasType);\r
+ hvar18->Draw("histo");\r
+ // - SOME PARAMETER VALUE\r
+ Bool_t kfitgauss = kFALSE;\r
+ Bool_t kfitleft = kFALSE;\r
+ Bool_t kfitright = kFALSE;\r
+ Int_t ptbinNarrowY = 0;\r
+ if (icasType < 2) ptbinNarrowY = 10; // 6;\r
+ else ptbinNarrowY = 3; // 2;\r
+ // - SOME DEFINITIONS\r
+ Float_t lowlimmass;\r
+ Float_t uplimmass;\r
+ Float_t lowgausslim;\r
+ Float_t upgausslim;\r
+ if (icasType==0||icasType==1) {\r
+ lowlimmass=1.30;\r
+ uplimmass=1.34;\r
+ lowgausslim=1.312;\r
+ upgausslim=1.332;\r
+ } else {\r
+ lowlimmass=1.645;\r
+ uplimmass=1.70;\r
+ lowgausslim=1.668;\r
+ upgausslim=1.678;\r
+ }\r
+ TF1* fitinvmass = new TF1("fitinvmass","gaus(0)+pol2(3)",lowlimmass,uplimmass);\r
+ fitinvmass->SetParName(0, "cnstntG");\r
+ fitinvmass->SetParName(1, "meanG");\r
+ fitinvmass->SetParName(2, "sigmaG");\r
+ fitinvmass->SetParLimits(0,0.,500000.);\r
+ if (icasType==0||icasType==1) {\r
+ fitinvmass->SetParameter(1, 1.32171);\r
+ fitinvmass->SetParLimits(1, 1.31,1.33);\r
+ fitinvmass->SetParLimits(2,0.001,0.005);\r
+ } else {\r
+ fitinvmass->SetParameter(1, 1.67245);\r
+ fitinvmass->SetParLimits(1, 1.664,1.68);\r
+ fitinvmass->SetParLimits(2,0.0008,0.006);\r
+ }\r
+ hvar18->Fit("fitinvmass","rimeN");\r
+ fitinvmass->SetLineColor(kRed);\r
+ fitinvmass->Draw("same");\r
+ Float_t meanGauss = fitinvmass->GetParameter(1);\r
+ Float_t sigmaGauss = fitinvmass->GetParameter(2);\r
+ cout<<"Mean: "<<meanGauss<<endl;\r
+ cout<<"Sigma: "<<sigmaGauss<<endl;\r
+ //Pad2: Text\r
+ c4->cd(2);\r
+ Float_t refwidth = 0.002;\r
+ TPaveText *pave1 = new TPaveText(0.05,0.3,0.95,0.5);\r
+ pave1->SetFillColor(0);\r
+ pave1->SetTextSize(0.04);\r
+ pave1->SetTextAlign(12);\r
+ if (icasType < 2) pave1->AddText("PDG mass: 1.32171 GeV/c^{2}");\r
+ else pave1->AddText("PDG mass: 1.67245 GeV/c^{2}");\r
+ pave1->AddText(Form("#color[1]{Mass form Fit: %.5f #pm %.5f GeV/c^{2}}",meanGauss,sigmaGauss));\r
+ if (sigmaGauss > refwidth - 0.0003 && sigmaGauss < refwidth + 0.0003) pave1->AddText("#color[3]{OK!! The width is compatible with standard.}");\r
+ else pave1->AddText("#color[2]{NOT OK!! Problem.}");\r
+ pave1->Draw();\r
+ cout<<" "<<refwidth - 0.0003<<"<"<<sigmaGauss<<"<"<<refwidth + 0.0003<<endl;\r
+ \r
+ //DEFINE 5st CANVAS AND DRAW PLOTS\r
+ if (collidingsystem == 0) {\r
+ TCanvas *c5 = new TCanvas("c5","",1200,270);\r
+ c5->Divide(2,1);\r
+ //Pad 1: centrality\r
+ c5->cd(1);\r
+ TH1D *hvar16 = cf->ShowProjection(19,icasType);\r
+ hvar16->Draw("histo");\r
+ //Pad 2: track multiplicity\r
+ c5->cd(2);\r
+ TH1D *hvar17 = cf->ShowProjection(20,icasType);\r
+ hvar17->Draw("histo");\r
+ c5->SaveAs("MultistrangeQA.pdf)");\r
+ }\r
+ \r
+\r
+}\r
+\r
+\r
+\r
+\r
+Bool_t checkUnderTheLimit(TH1D *lHist, Double_t limit) {\r
+\r
+ Int_t binlimit = lHist->FindBin(limit);\r
+ Bool_t checkOk = kTRUE;\r
+ for (Int_t i = 1; i < binlimit; i++) {\r
+ Int_t content = 0;\r
+ content = lHist->GetBinContent(i);\r
+ if (content != 0) checkOk = kFALSE;\r
+ //cout<<"Content bin "<<i<<": "<<content<<endl;\r
+ }\r
+ return checkOk;\r
+\r
+}\r
+\r
+\r
+Bool_t checkOverTheLimit(TH1D *lHist, Double_t limit) {\r
+\r
+ Int_t binlimit = lHist->FindBin(limit);\r
+ Int_t lastbin = lHist->GetNbinsX();\r
+ Bool_t checkOk = kTRUE;\r
+ for (Int_t i = binlimit; i < lastbin+1; i++) {\r
+ Int_t content = 0;\r
+ content = lHist->GetBinContent(i);\r
+ if (content != 0) checkOk = kFALSE;\r
+ //cout<<"Content bin "<<i<<": "<<content<<endl;\r
+ }\r
+ return checkOk;\r
+\r
+}\r
+\r
+\r
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff