[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / QAPlotsBoth.C

Float_t QAPlotsBoth( AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,
	      AliSpectraBothEventCuts* ecuts_data, AliSpectraBothEventCuts* ecuts_mc,
	      AliSpectraBothTrackCuts* tcuts_data, AliSpectraBothTrackCuts* tcuts_mc,
	      TList * flistqa,TList * flistcanvas,Bool_t fullicorr=kTRUE)
{
TString pidmethods[3]={"TPC","TOF","TPCTOF"};	
	Double_t neventsdata =  ecutsdata->NumberOfPhysSelEvents();
	Double_t neventsmc =  ecutsmc->NumberOfPhysSelEvents();
	
	
	for(Int_t ipart=0;ipart<3;ipart++)
	{
			
			for(Int_t imethod=0;imethod<3;imethod++)
			{
				if(!hman_data->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
					continue;
				 TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
				// nsig_data->RebinX(20);			 
				// nsig_data->RebinY(4);
				// nsig_data->Sumw2();
				if(!nsig_data)
					continue;

				if(!hman_mc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
					continue;
				 TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
				 //nsig_mc->RebinX(20);			 
				// nsig_mc->RebinY(4);
				// nsig_mc->Sumw2();
				if(!nsig_mc)
					continue;
				 

				Int_t ibin=1;
				Float_t binsize=nsig_mc->GetXaxis()->GetBinWidth(1);

				TH1F* maxposdata=(TH1F*)nsig_data->ProjectionX(Form("%s%sdatamaxpos",Particle[ipart].Data(),pidmethods[imethod].Data()),-1,-1);
				maxposdata->Reset();
				maxposdata->SetTitle(";p_{T} (GeV/c);max in (-2,2)");
				TH1F* maxposmc=(TH1F*)nsig_data->ProjectionX(Form("%s%smcmaxpos",Particle[ipart].Data(),pidmethods[imethod].Data()),-1,-1);
				maxposmc->Reset();
				maxposmc->SetTitle(";p_{T} (GeV/c);max in (-2,2)");

	
				 while (ibin*binsize<3.0)
				 {
					// TCanvas* c=new TCanvas(Form("canvas%s%s%d",Particle[ipart].Data(),pidmethods[imethod].Data(),ibin),Form("canvas%s%s%d",Particle[ipart].Data(),pidmethods[imethod].Data(),ibin),700,500);
					
					 TH1F *nsig_data_Proj1=(TH1F*)nsig_data->ProjectionY(Form("%s%sdata[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
					 TH1F *nsig_mc_Proj1=(TH1F*)nsig_mc->ProjectionY(Form("%s%smc[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_mc->GetXaxis()->GetBinLowEdge(ibin),nsig_mc->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
					 nsig_data_Proj1->GetXaxis()->SetRangeUser(-3,3);
					 nsig_data_Proj1->SetLineColor(kRed);
					 if(nsig_data_Proj1->Integral()<1&&nsig_mc_Proj1->Integral()<1)
					 {	
						ibin++;	
						continue;
					 } 
									//	nsig_data_Proj1->Sumw2();
				//	nsig_mc_Proj1->Sumw2();
					  //nsig_data_Proj1->GetXaxis()->SetRangeUser(-5,5);
					 
					 //c->cd()->SetLogy();
					 //nsig_data_Proj1->Draw();
					 //nsig_mc_Proj1->Draw("same");
					nsig_data_Proj1->GetXaxis()->SetRange(nsig_data_Proj1->GetXaxis()->FindBin(-2.0),nsig_data_Proj1->GetXaxis()->FindBin(2.0));
					if(nsig_data_Proj1->GetMaximumBin()<=nsig_data_Proj1->GetXaxis()->FindBin(2.0)&&nsig_data_Proj1->GetMaximumBin()>=nsig_data_Proj1->GetXaxis()->FindBin(-2.0))
					{
						maxposdata->SetBinContent(ibin,nsig_data_Proj1->GetXaxis()->GetBinCenter(nsig_data_Proj1->GetMaximumBin()));	
						maxposdata->SetBinError(ibin,nsig_data_Proj1->GetXaxis()->GetBinWidth(nsig_data_Proj1->GetMaximumBin())/2.0);	
	 				}
				//cout<<Form("%s%sdatamaxpos",Particle[ipart].Data(),pidmethods[imethod].Data())<<" "<<nsig_data_Proj1->GetMaximumBin()<<" "<<nsig_data_Proj1->GetXaxis()->FindBin(2.0)<<" "<<nsig_data_Proj1->GetXaxis()->FindBin(-2.0)<<" "<<nsig_data_Proj1->GetXaxis()->GetBinCenter(nsig_data_Proj1->GetMaximumBin())<<" "<<ibin<<endl;

					nsig_data_Proj1->GetXaxis()->SetRange(0,nsig_data_Proj1->GetXaxis()->GetNbins());
									
					nsig_mc_Proj1->GetXaxis()->SetRange(nsig_mc_Proj1->GetXaxis()->FindBin(-2.0),nsig_mc_Proj1->GetXaxis()->FindBin(2.0));
					if(nsig_mc_Proj1->GetMaximumBin()<=nsig_mc_Proj1->GetXaxis()->FindBin(2.0)&&nsig_mc_Proj1->GetMaximumBin()>=nsig_mc_Proj1->GetXaxis()->FindBin(-2.0))
					{
						maxposmc->SetBinContent(ibin,nsig_mc_Proj1->GetXaxis()->GetBinCenter(nsig_mc_Proj1->GetMaximumBin()));	
						maxposmc->SetBinError(ibin,nsig_mc_Proj1->GetXaxis()->GetBinWidth(nsig_mc_Proj1->GetMaximumBin())/2.0);	
	 				}
				//		cout<<Form("%s%smcmaxpos",Particle[ipart].Data(),pidmethods[imethod].Data())<<" "<<nsig_mc_Proj1->GetMaximumBin()<<" "<<nsig_mc_Proj1->GetXaxis()->FindBin(2.0)<<" "<<nsig_mc_Proj1->GetXaxis()->FindBin(-2.0)<<" "<<ibin<<" "<<nsig_mc_Proj1->GetXaxis()->GetBinCenter(nsig_mc_Proj1->GetMaximumBin())<<endl;

					nsig_mc_Proj1->GetXaxis()->SetRange(0,nsig_mc_Proj1->GetXaxis()->GetNbins());

					if(neventsdata>0.0&&neventsmc>0.0)
					{
						nsig_data_Proj1->Sumw2();
						nsig_mc_Proj1->Sumw2();
						nsig_data_Proj1->Scale(1.0/neventsdata);
					 	nsig_mc_Proj1->Scale(1.0/neventsmc);
					}
					
					flistqa->Add(nsig_data_Proj1);
					flistqa->Add(nsig_mc_Proj1);
					ibin++;
				 }
				flistqa->Add(maxposmc);
				flistqa->Add(maxposdata);
			}
	}
	TH1F* fHistoVtxAftSeldata=(TH1F*)ecuts_data->GetHistoVtxAftSel();
	TH1F* fHistoVtxAftSelmc=(TH1F*)ecuts_mc->GetHistoVtxAftSel();
	Int_t binstartx=1;	
         while(fHistoVtxAftSeldata->GetBinContent(binstartx)<1&&binstartx<fHistoVtxAftSeldata->GetXaxis()->GetNbins())
                binstartx++;
        binstartx++;
	flistcanvas->Add(plot_on_canvas("vertex",fHistoVtxAftSeldata,fHistoVtxAftSelmc));
	/*

	TF1* fdata=new TF1("dataveretxfit","gausn");
	TF1* fmc=new TF1("mcveretxfit","gausn");
	//we strat fit a second not empty bin
	Float_t minfit=fHistoVtxAftSeldata->GetXaxis()->GetBinCenter(binstartx);
	cout<<"fit starts "<<minfit<<endl;
	fHistoVtxAftSeldata->Fit("dataveretxfit","0","",minfit,-1.0*minfit);
	fHistoVtxAftSelmc->Fit("mcveretxfit","0","",minfit,-1.0*minfit);
	Float_t datavertexratio=fHistoVtxAftSeldata->Integral(-1,-1,"width")/fdata->GetParameter(0);
	Float_t mcvertexratio=fHistoVtxAftSelmc->Integral(-1,-1,"width")/fmc->GetParameter(0);
	*/
	//Event cut histo 
	TH1I* histodata=ecuts_data->GetHistoCuts();
	TH1I* histomc=ecuts_mc->GetHistoCuts();
	
	Int_t events_data=histodata->GetBinContent(3);
	Int_t events_mc=histomc->GetBinContent(3);

	if(events_data==0&&events_mc==0)
		return 0;


	Float_t datavertexratio=((Float_t)(events_data))/((Float_t)histodata->GetBinContent(4));
         Float_t mcvertexratio=((Float_t)(events_mc))/((Float_t)histomc->GetBinContent(4));
	 TH1F* fHistoEtaAftSeldata=(TH1F*)ecuts_data->GetHistoEtaAftSel();
	 TH1F* fHistoEtaAftSelmc=(TH1F*)ecuts_mc->GetHistoEtaAftSel();
	flistcanvas->Add(plot_on_canvas("ETA",fHistoEtaAftSeldata,fHistoEtaAftSelmc));


	 TH1F* fITSclustershistdata=(TH1F*)tcuts_data->GetHistoNclustersITS();
 	  TH1F* fITSclustershistmc=(TH1F*)tcuts_mc->GetHistoNclustersITS();

	flistcanvas->Add(plot_on_canvas("NITS",fITSclustershistdata,fITSclustershistmc));
	cout<<" data "<<datavertexratio<<" mc "<<mcvertexratio<<endl;
	
	TH2F* hmul=(TH2F*)hman_mc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul");	
	hmul->Sumw2();	
	TCanvas* cbc=new TCanvas("broken chunks","broken chunks",1200,600);
	cbc->Divide(2,1);
	cbc->cd(1);
	hmul->Draw();
	cbc->cd(2);
	TH1F* nonzero=(TH1F*)hmul->ProjectionX("nonzerotracks",2,-1);
	nonzero->SetMarkerColor(kRed);
	nonzero->SetMarkerStyle(21);
	TH1F* binzero=(TH1F*)hmul->ProjectionX("binzerotracks",1,1);
	binzero->SetMarkerColor(kBlack);
	binzero->SetMarkerStyle(22);
	binzero->Sumw2();
	nonzero->Sumw2();
	binzero->Divide(nonzero);
	TF1* badchunk=new TF1("badchunkfit","pol0",10,40);
	binzero->Fit("badchunkfit","R");
	Float_t badchunksfraction=badchunk->GetParameter(0);
	cout<<"Bad chunks "<<badchunksfraction<<endl;	
	binzero->Draw("E1");
	flistcanvas->Add(cbc);
	if(TMath::Abs(hmul->GetEntries()/events_mc-1.0)>0.001)
		cout<<"MC merging problem"<<endl;
	if(TMath::Abs(hman_data->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()/events_data-1.0)>0.001)
		cout<<"Data merging problem"<<endl;
	
	if(fullicorr)
		return (1.0-badchunksfraction)*mcvertexratio/datavertexratio;
	else 
		return (1.0-badchunksfraction)*mcvertexratio;
}
TCanvas* plot_on_canvas(TString name, TH1* h1,TH1* h2)
{
	TCanvas* cvrt=new TCanvas(name.Data(),name.Data(),600,600);
	cvrt->cd();
	h1->SetLineColor(kRed);
	h2->SetLineColor(kBlue);
	h1->Sumw2();
	h2->Sumw2();
	TLegend *lvtr=new TLegend(0.2,0.2,0.5,0.3,"","NDC");
	lvtr->SetLineColor(kWhite);
	lvtr->AddEntry(h1,"data","l");
	lvtr->AddEntry(h2,"MC","l");
	if(h1->GetBinContent(h1->GetXaxis()->FindBin(0.0))>0.0&&h2->GetBinContent(h2->GetXaxis()->FindBin(0.0))>0.0)
	{
		h1->Scale(1.0/h1->GetBinContent(h1->GetXaxis()->FindBin(0.0)));
		h2->Scale(1.0/h2->GetBinContent(h2->GetXaxis()->FindBin(0.0)));
        }
	h1->DrawCopy("L");
	h2->DrawCopy("Lsame");
	lvtr->Draw();
	return cvrt;
}
Commit	Line	Data
88d02ae4	1	Float_t QAPlotsBoth( AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,
	2	AliSpectraBothEventCuts* ecuts_data, AliSpectraBothEventCuts* ecuts_mc,
	3	AliSpectraBothTrackCuts* tcuts_data, AliSpectraBothTrackCuts* tcuts_mc,
3a71f081	4	TList * flistqa,TList * flistcanvas,Bool_t fullicorr=kTRUE)
88d02ae4	5	{
	6	TString pidmethods[3]={"TPC","TOF","TPCTOF"};
	7	Double_t neventsdata = ecutsdata->NumberOfPhysSelEvents();
	8	Double_t neventsmc = ecutsmc->NumberOfPhysSelEvents();
	9
eb9ebb0a	10
eb9ebb0a	11
88d02ae4	12	for(Int_t ipart=0;ipart<3;ipart++)
88d02ae4	13	{
eb9ebb0a	14
88d02ae4	15	for(Int_t imethod=0;imethod<3;imethod++)
88d02ae4	16	{
5c3e7f80	17	if(!hman_data->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
5c3e7f80	18	continue;
88d02ae4	19	TH2F nsig_data = (TH2F)((TH2F*)hman_data->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
	20	// nsig_data->RebinX(20);
	21	// nsig_data->RebinY(4);
	22	// nsig_data->Sumw2();
be25efef	23	if(!nsig_data)
be25efef	24	continue;
88d02ae4	25
5c3e7f80	26	if(!hman_mc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
5c3e7f80	27	continue;
88d02ae4	28	TH2F nsig_mc = (TH2F)((TH2F*)hman_mc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
	29	//nsig_mc->RebinX(20);
	30	// nsig_mc->RebinY(4);
	31	// nsig_mc->Sumw2();
be25efef	32	if(!nsig_mc)
be25efef	33	continue;
88d02ae4	34
be25efef	35
88d02ae4	36	Int_t ibin=1;
88d02ae4	37	Float_t binsize=nsig_mc->GetXaxis()->GetBinWidth(1);
eb9ebb0a	38
cd28560b	39	TH1F* maxposdata=(TH1F*)nsig_data->ProjectionX(Form("%s%sdatamaxpos",Particle[ipart].Data(),pidmethods[imethod].Data()),-1,-1);
eb9ebb0a	40	maxposdata->Reset();
eb9ebb0a	41	maxposdata->SetTitle(";p_{T} (GeV/c);max in (-2,2)");
cd28560b	42	TH1F* maxposmc=(TH1F*)nsig_data->ProjectionX(Form("%s%smcmaxpos",Particle[ipart].Data(),pidmethods[imethod].Data()),-1,-1);
eb9ebb0a	43	maxposmc->Reset();
	44	maxposmc->SetTitle(";p_{T} (GeV/c);max in (-2,2)");
	45
	46
88d02ae4	47	while (ibin*binsize<3.0)
	48	{
	49	// TCanvas* c=new TCanvas(Form("canvas%s%s%d",Particle[ipart].Data(),pidmethods[imethod].Data(),ibin),Form("canvas%s%s%d",Particle[ipart].Data(),pidmethods[imethod].Data(),ibin),700,500);
	50
cd28560b	51	TH1F nsig_data_Proj1=(TH1F)nsig_data->ProjectionY(Form("%s%sdata[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
cd28560b	52	TH1F nsig_mc_Proj1=(TH1F)nsig_mc->ProjectionY(Form("%s%smc[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_mc->GetXaxis()->GetBinLowEdge(ibin),nsig_mc->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
88d02ae4	53	nsig_data_Proj1->GetXaxis()->SetRangeUser(-3,3);
	54	nsig_data_Proj1->SetLineColor(kRed);
	55	if(nsig_data_Proj1->Integral()<1&&nsig_mc_Proj1->Integral()<1)
	56	{
	57	ibin++;
	58	continue;
	59	}
eb9ebb0a	60	// nsig_data_Proj1->Sumw2();
88d02ae4	61	// nsig_mc_Proj1->Sumw2();
	62	//nsig_data_Proj1->GetXaxis()->SetRangeUser(-5,5);
	63
	64	//c->cd()->SetLogy();
	65	//nsig_data_Proj1->Draw();
	66	//nsig_mc_Proj1->Draw("same");
eb9ebb0a	67	nsig_data_Proj1->GetXaxis()->SetRange(nsig_data_Proj1->GetXaxis()->FindBin(-2.0),nsig_data_Proj1->GetXaxis()->FindBin(2.0));
	68	if(nsig_data_Proj1->GetMaximumBin()<=nsig_data_Proj1->GetXaxis()->FindBin(2.0)&&nsig_data_Proj1->GetMaximumBin()>=nsig_data_Proj1->GetXaxis()->FindBin(-2.0))
	69	{
	70	maxposdata->SetBinContent(ibin,nsig_data_Proj1->GetXaxis()->GetBinCenter(nsig_data_Proj1->GetMaximumBin()));
	71	maxposdata->SetBinError(ibin,nsig_data_Proj1->GetXaxis()->GetBinWidth(nsig_data_Proj1->GetMaximumBin())/2.0);
	72	}
a4ca64e3	73	//cout<<Form("%s%sdatamaxpos",Particle[ipart].Data(),pidmethods[imethod].Data())<<" "<<nsig_data_Proj1->GetMaximumBin()<<" "<<nsig_data_Proj1->GetXaxis()->FindBin(2.0)<<" "<<nsig_data_Proj1->GetXaxis()->FindBin(-2.0)<<" "<<nsig_data_Proj1->GetXaxis()->GetBinCenter(nsig_data_Proj1->GetMaximumBin())<<" "<<ibin<<endl;
eb9ebb0a	74
	75	nsig_data_Proj1->GetXaxis()->SetRange(0,nsig_data_Proj1->GetXaxis()->GetNbins());
	76
	77	nsig_mc_Proj1->GetXaxis()->SetRange(nsig_mc_Proj1->GetXaxis()->FindBin(-2.0),nsig_mc_Proj1->GetXaxis()->FindBin(2.0));
	78	if(nsig_mc_Proj1->GetMaximumBin()<=nsig_mc_Proj1->GetXaxis()->FindBin(2.0)&&nsig_mc_Proj1->GetMaximumBin()>=nsig_mc_Proj1->GetXaxis()->FindBin(-2.0))
	79	{
	80	maxposmc->SetBinContent(ibin,nsig_mc_Proj1->GetXaxis()->GetBinCenter(nsig_mc_Proj1->GetMaximumBin()));
	81	maxposmc->SetBinError(ibin,nsig_mc_Proj1->GetXaxis()->GetBinWidth(nsig_mc_Proj1->GetMaximumBin())/2.0);
	82	}
a4ca64e3	83	// cout<<Form("%s%smcmaxpos",Particle[ipart].Data(),pidmethods[imethod].Data())<<" "<<nsig_mc_Proj1->GetMaximumBin()<<" "<<nsig_mc_Proj1->GetXaxis()->FindBin(2.0)<<" "<<nsig_mc_Proj1->GetXaxis()->FindBin(-2.0)<<" "<<ibin<<" "<<nsig_mc_Proj1->GetXaxis()->GetBinCenter(nsig_mc_Proj1->GetMaximumBin())<<endl;
eb9ebb0a	84
	85	nsig_mc_Proj1->GetXaxis()->SetRange(0,nsig_mc_Proj1->GetXaxis()->GetNbins());
	86
1e6bff53	87	if(neventsdata>0.0&&neventsmc>0.0)
1e6bff53	88	{
50c532c4	89	nsig_data_Proj1->Sumw2();
50c532c4	90	nsig_mc_Proj1->Sumw2();
1e6bff53	91	nsig_data_Proj1->Scale(1.0/neventsdata);
	92	nsig_mc_Proj1->Scale(1.0/neventsmc);
	93	}
eb9ebb0a	94
33d374d5	95	flistqa->Add(nsig_data_Proj1);
33d374d5	96	flistqa->Add(nsig_mc_Proj1);
88d02ae4	97	ibin++;
88d02ae4	98	}
eb9ebb0a	99	flistqa->Add(maxposmc);
eb9ebb0a	100	flistqa->Add(maxposdata);
88d02ae4	101	}
	102	}
	103	TH1F* fHistoVtxAftSeldata=(TH1F*)ecuts_data->GetHistoVtxAftSel();
	104	TH1F* fHistoVtxAftSelmc=(TH1F*)ecuts_mc->GetHistoVtxAftSel();
547e52a0	105	Int_t binstartx=1;
	106	while(fHistoVtxAftSeldata->GetBinContent(binstartx)<1&&binstartx<fHistoVtxAftSeldata->GetXaxis()->GetNbins())
	107	binstartx++;
	108	binstartx++;
33d374d5	109	flistcanvas->Add(plot_on_canvas("vertex",fHistoVtxAftSeldata,fHistoVtxAftSelmc));
3a71f081	110	/*
3a71f081	111
88d02ae4	112	TF1* fdata=new TF1("dataveretxfit","gausn");
88d02ae4	113	TF1* fmc=new TF1("mcveretxfit","gausn");
547e52a0	114	//we strat fit a second not empty bin
	115	Float_t minfit=fHistoVtxAftSeldata->GetXaxis()->GetBinCenter(binstartx);
	116	cout<<"fit starts "<<minfit<<endl;
	117	fHistoVtxAftSeldata->Fit("dataveretxfit","0","",minfit,-1.0*minfit);
	118	fHistoVtxAftSelmc->Fit("mcveretxfit","0","",minfit,-1.0*minfit);
88d02ae4	119	Float_t datavertexratio=fHistoVtxAftSeldata->Integral(-1,-1,"width")/fdata->GetParameter(0);
88d02ae4	120	Float_t mcvertexratio=fHistoVtxAftSelmc->Integral(-1,-1,"width")/fmc->GetParameter(0);
3a71f081	121	*/
	122	//Event cut histo
	123	TH1I* histodata=ecuts_data->GetHistoCuts();
	124	TH1I* histomc=ecuts_mc->GetHistoCuts();
88d02ae4	125
3a71f081	126	Int_t events_data=histodata->GetBinContent(3);
	127	Int_t events_mc=histomc->GetBinContent(3);
	128
	129	if(events_data==0&&events_mc==0)
	130	return 0;
88d02ae4	131
3a71f081	132
	133	Float_t datavertexratio=((Float_t)(events_data))/((Float_t)histodata->GetBinContent(4));
	134	Float_t mcvertexratio=((Float_t)(events_mc))/((Float_t)histomc->GetBinContent(4));
88d02ae4	135	TH1F* fHistoEtaAftSeldata=(TH1F*)ecuts_data->GetHistoEtaAftSel();
88d02ae4	136	TH1F* fHistoEtaAftSelmc=(TH1F*)ecuts_mc->GetHistoEtaAftSel();
33d374d5	137	flistcanvas->Add(plot_on_canvas("ETA",fHistoEtaAftSeldata,fHistoEtaAftSelmc));
88d02ae4	138
	139
	140	TH1F* fITSclustershistdata=(TH1F*)tcuts_data->GetHistoNclustersITS();
	141	TH1F* fITSclustershistmc=(TH1F*)tcuts_mc->GetHistoNclustersITS();
	142
33d374d5	143	flistcanvas->Add(plot_on_canvas("NITS",fITSclustershistdata,fITSclustershistmc));
88d02ae4	144	cout<<" data "<<datavertexratio<<" mc "<<mcvertexratio<<endl;
	145
	146	TH2F* hmul=(TH2F*)hman_mc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul");
	147	hmul->Sumw2();
	148	TCanvas* cbc=new TCanvas("broken chunks","broken chunks",1200,600);
	149	cbc->Divide(2,1);
	150	cbc->cd(1);
	151	hmul->Draw();
	152	cbc->cd(2);
	153	TH1F* nonzero=(TH1F*)hmul->ProjectionX("nonzerotracks",2,-1);
	154	nonzero->SetMarkerColor(kRed);
	155	nonzero->SetMarkerStyle(21);
	156	TH1F* binzero=(TH1F*)hmul->ProjectionX("binzerotracks",1,1);
	157	binzero->SetMarkerColor(kBlack);
	158	binzero->SetMarkerStyle(22);
	159	binzero->Sumw2();
	160	nonzero->Sumw2();
	161	binzero->Divide(nonzero);
	162	TF1* badchunk=new TF1("badchunkfit","pol0",10,40);
	163	binzero->Fit("badchunkfit","R");
	164	Float_t badchunksfraction=badchunk->GetParameter(0);
547e52a0	165	cout<<"Bad chunks "<<badchunksfraction<<endl;
88d02ae4	166	binzero->Draw("E1");
33d374d5	167	flistcanvas->Add(cbc);
3a71f081	168	if(TMath::Abs(hmul->GetEntries()/events_mc-1.0)>0.001)
	169	cout<<"MC merging problem"<<endl;
	170	if(TMath::Abs(hman_data->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()/events_data-1.0)>0.001)
	171	cout<<"Data merging problem"<<endl;
88d02ae4	172
3a71f081	173	if(fullicorr)
	174	return (1.0-badchunksfraction)*mcvertexratio/datavertexratio;
	175	else
	176	return (1.0-badchunksfraction)*mcvertexratio;
88d02ae4	177	}
	178	TCanvas* plot_on_canvas(TString name, TH1* h1,TH1* h2)
	179	{
	180	TCanvas* cvrt=new TCanvas(name.Data(),name.Data(),600,600);
	181	cvrt->cd();
	182	h1->SetLineColor(kRed);
	183	h2->SetLineColor(kBlue);
	184	h1->Sumw2();
	185	h2->Sumw2();
	186	TLegend *lvtr=new TLegend(0.2,0.2,0.5,0.3,"","NDC");
	187	lvtr->SetLineColor(kWhite);
	188	lvtr->AddEntry(h1,"data","l");
	189	lvtr->AddEntry(h2,"MC","l");
1e6bff53	190	if(h1->GetBinContent(h1->GetXaxis()->FindBin(0.0))>0.0&&h2->GetBinContent(h2->GetXaxis()->FindBin(0.0))>0.0)
	191	{
	192	h1->Scale(1.0/h1->GetBinContent(h1->GetXaxis()->FindBin(0.0)));
	193	h2->Scale(1.0/h2->GetBinContent(h2->GetXaxis()->FindBin(0.0)));
	194	}
88d02ae4	195	h1->DrawCopy("L");
	196	h2->DrawCopy("Lsame");
	197	lvtr->Draw();
	198	return cvrt;
	199	}
	200