[u/mrichter/AliRoot.git] / PHOS / AnaPhotons_v1.C

//===============================================================
// In this Macro (which can (must) be compilated), you will find all the 
// analysis functions to build photon spectrum, invariant mass 
// spectrum of photon  pairs and combinatorial background calculations 
// in ALICE electromagnetic calorimeter
// Author: Gines MARTINEZ, Subatech, 15 june 2001
//============================================================== 
#include "TH2.h"
#include "TH1.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TTree.h"
#include "TRandom.h"
#include "TObjectTable.h"
#include "AliRun.h"
#include "AliPHOSGetter.h"
#include "AliPHOSRecParticle.h"
#include "TLorentzVector.h"
#include "TGraphErrors.h"
#include "TF1.h"

TObjectTable * gObjectTable;
TRandom      * gRandom;
AliRun       * gAlice;


Bool_t PhotonId(AliPHOSRecParticle * RecParticle, char *opt)
{
  Bool_t PhotonId = kTRUE ;
  if (strstr(opt,"narrow"))
    {
      if ( RecParticle->GetType() & 1) PhotonId= kFALSE ;
    }
  if (strstr(opt,"cpv"))
    {
      if (  (RecParticle->GetType() &  1<<2) >>2   ) PhotonId= kFALSE ;
    }
  if (strstr(opt,"time"))
    {
      
    }

  return PhotonId;
}

void AnaMinv(char * filename,char * opt)
{
  TH2F * h_Minv_lowpT       = new TH2F("h_Minv_lowpT","Minv vs pT low",500,0.0,1.0,40,0.,10.);
  TH2F * h_Minv_highpT      = new TH2F("h_Minv_highpT","Minv vs pT high",500,0.0,1.0,50,0.,100.);
  TH2F * h_Minv_lowpT_back  = new TH2F("h_Minv_lowpT_back","Minv vs pT low back",500,0.0,1.0,40,0.,10.);
  TH2F * h_Minv_highpT_back = new TH2F("h_Minv_highpT_back","Minv vs pT high back",500,0.0,1.0,50,0.,100.);


  TH1F * h_Pseudoeta = new TH1F("h_Pseudoeta","Pseudoeta photons",500,-1.0,1.0);
  TH1F * h_Pt        = new TH1F("h_Pt","Pt photons",400,0.,10.);
  TH2F * h_Peta_Pt   = new TH2F("h_Peta_Pt","Pseudo vs pT",40,0.,10.,50,-1.0,1.0);
  TH1F * h_Phi       = new TH1F("h_Phi","Phi photons",400,-4.,4.);
  TH2F * h_Peta_Phi  = new TH2F("h_Peta_Phi","Pseudo vs Phi",200,-4,4,200,-1.0,1.0);
  TH1F * h_Dispersion= new TH1F("h_Dispersion","Dispersion",50,0.,10.);
  TH1F * h_Type      = new TH1F("h_Type","Particle Type",9,-1.5,7.5);

  TH1F * h_DeltaR   = new TH1F("h_DeltaR","Delta R",400,0.,2.);
  TH1F * h_Asymmetry= new TH1F("h_Asymmetry","Asymmetry",400, -2., 2.);

  AliPHOSGetter * RecData = AliPHOSGetter::GetInstance(filename,"Gines")  ;
 
  AliPHOSRecParticle * RecParticle1;
  AliPHOSRecParticle * RecParticle2;

  Float_t RelativeRCut = 0.00001 ; 
  Float_t AsymmetryCut = 0.7 ;
  Float_t Asymmetry;
  Float_t Type;

  Int_t iEvent, iRecParticle1, iRecParticle2;
  Int_t nRecParticle;
  Float_t invariant_mass, invariant_mass_mixed;
  TClonesArray * RecParticles;
 
  Float_t average_multiplicity = 0.;

  for(iEvent=0; iEvent<gAlice->TreeE()->GetEntries(); iEvent++)
    {
      TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;

      RecData->Event(iEvent,"R");
      printf(">>> Event %d \n",iEvent);
      RecParticles = RecData->RecParticles();
      nRecParticle = RecParticles->GetEntries();
      printf(">>> >> Recparticle multilicity is %d \n",nRecParticle);
      average_multiplicity += ((Float_t) (nRecParticle) ) / ( (Float_t)gAlice->TreeE()->GetEntries() ) ;
      // Construction de la masse invariante des pairs
      if (nRecParticle > 1) 
	{
	  for(iRecParticle1=0; iRecParticle1<nRecParticle; iRecParticle1++)
	    {
	      RecParticle1 =  (AliPHOSRecParticle *) RecParticles->At(iRecParticle1);
	      RecParticle1->Momentum(P_photon1);
	      printf(">>> >> Photon momentum is %f \n", P_photon1.Pt() ); 
	      Type = RecParticle1->GetType();
	      h_Type->Fill(Type);
	      // Photon Id check
	      if ( PhotonId(RecParticle1, opt) )
		{
		  
		  h_Pseudoeta->Fill(P_photon1.PseudoRapidity());
		  h_Pt->Fill(P_photon1.Pt());
		  h_Phi->Fill(P_photon1.Phi());
		  h_Peta_Pt->Fill(P_photon1.Pt(), P_photon1.PseudoRapidity());
		  h_Peta_Phi->Fill(P_photon1.Phi(), P_photon1.PseudoRapidity() );
	    
		  for(iRecParticle2=iRecParticle1+1; iRecParticle2<nRecParticle; iRecParticle2++)
		    {
		      RecParticle2 = (AliPHOSRecParticle *) RecParticles->At(iRecParticle2);
		      RecParticle2->Momentum(P_photon2); 
		      Asymmetry = TMath::Abs((P_photon1.E()-P_photon2.E())/(P_photon1.E()+P_photon2.E()));
		      //Photon Id Check
		      if ( PhotonId(RecParticle2, opt) )
			{
			  if ( (P_photon1 != P_photon2) && 
			       (P_photon1.DeltaR(P_photon2) > RelativeRCut) &&
			       (Asymmetry < AsymmetryCut)                          )
			    {
			      h_DeltaR->Fill(P_photon1.DeltaR(P_photon2));
			      h_Asymmetry->Fill( Asymmetry );
			      
			      //   printf("A. p1 es %f \n",P_photon1->E());
			      invariant_mass = (P_photon1 + P_photon2).M();
			      // printf("B. p1 es %f \n",P_photon1->E());
			      h_Minv_lowpT->Fill(invariant_mass, (P_photon1 + P_photon2).Pt() );
			      h_Minv_highpT->Fill(invariant_mass,(P_photon1 + P_photon2).Pt() );
			    } 
			}
		    }
 		}
 	    }
	}
    }
  printf(">>> Average Multiplicity is %f \n",average_multiplicity);
  Int_t Background = (Int_t)  (gAlice->TreeE()->GetEntries() * average_multiplicity * (average_multiplicity-1.)/2.) ;
  printf(">>> Background is %d \n",Background);

  Double_t Pt_Mixed1, Pt_Mixed2;
  Double_t Y_Mixed1, Y_Mixed2;
  Double_t Phi_Mixed1, Phi_Mixed2;

  for(iEvent=0; iEvent<Background; iEvent++)
    {
      TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
      //printf(">>> Background Event %d \n",iEvent);
      Pt_Mixed1 =  h_Pt->GetRandom(); 
      Pt_Mixed2 =  h_Pt->GetRandom();
      h_Peta_Phi->GetRandom2(Phi_Mixed1, Y_Mixed1);
      h_Peta_Phi->GetRandom2(Phi_Mixed2, Y_Mixed2);
      P_photonMixed1.SetPtEtaPhiM( Pt_Mixed1, Y_Mixed1, Phi_Mixed1, 0.0);
      P_photonMixed2.SetPtEtaPhiM( Pt_Mixed2, Y_Mixed2, Phi_Mixed2, 0.0);
      Asymmetry = TMath::Abs((P_photonMixed1.E()-P_photonMixed2.E())/(P_photonMixed1.E()+P_photonMixed2.E()));
      
      if ( (P_photonMixed1.DeltaR(P_photonMixed2) > RelativeRCut) &&
           (Asymmetry < AsymmetryCut  )                               )
	{
	  invariant_mass_mixed = (P_photonMixed1 + P_photonMixed2).M();
	  h_Minv_lowpT_back->Fill(invariant_mass_mixed, (P_photonMixed1 + P_photonMixed2).Pt() );
	  h_Minv_highpT_back->Fill(invariant_mass_mixed,(P_photonMixed1 + P_photonMixed2).Pt() );
	}  
    }
  printf (">>> Background calculation finished ! \n");

  char outputname[80];
  sprintf(outputname,"%s.Minv%s",filename,opt);
  TFile output(outputname,"recreate");
  h_Minv_lowpT->Write();
  h_Minv_highpT->Write();
  h_Minv_lowpT_back->Write();  
  h_Minv_highpT_back->Write();
  h_Pseudoeta->Write();
  h_Pt->Write();
  h_Peta_Pt->Write();
  h_Phi->Write();
  h_Peta_Phi->Write();
  h_Dispersion->Write();
  h_Type->Write();
  h_Asymmetry->Write();
  h_DeltaR->Write();
  
  output.Close();
}


void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, Option_t * opt)
{

  Int_t NumberOfPtBins = NumberPerPtBin;
  Float_t PtCalibration = 0.250;

  TFile * in = new TFile(filename);

  TH2F * h_Minv_pT = 0;  
  TH2F * h_Minv_pT_back = 0; 
  TH2F * frame = 0 ;

  if (strstr(opt,"low"))
    {
      h_Minv_pT      = (TH2F *) in->Get("h_Minv_lowpT"); ;
      h_Minv_pT_back = (TH2F *) in->Get("h_Minv_lowpT_back");
      PtCalibration  = 0.250;   
      frame = new TH2F("PtSpectrumlow","Pt Spectrum low",10, 0.,10.,10,0.1,10000);
    }  
  if (strstr(opt,"high"))
    {
      h_Minv_pT      = (TH2F *) in->Get("h_Minv_highpT"); ;
      h_Minv_pT_back = (TH2F *) in->Get("h_Minv_highpT_back");
      PtCalibration = 2.5;
      frame = new TH2F("PtSpectrumhigh","Pt Spectrum high",10, 0.,100.,10,0.1,10000);
    }

  if ( h_Minv_pT == 0 ) 
    {
      printf(">>> Bad Option! \n");
      return;
    }
  Int_t Norma_1 = 100; Float_t Norma_minv_1 = 0.2;
  Int_t Norma_2 = 200; Float_t Norma_minv_2 = 0.4;

  Int_t Minv_1 = 46;
  Int_t Minv_2 = 76;
  if (strstr(particle,"eta"))
    {
      Minv_1 = 214;
      Minv_2 = 314;
    }

  if (strstr(particle,"norma"))
    {
      Minv_1 = 100;
      Minv_2 = 200;
    }
  
  Int_t NHistos = 40/NumberOfPtBins;
  Int_t iHistos;

  TH1D * signal = 0;
  TH1D * background = 0;
  TH1D * ratio = 0;
  TH1D * difference = 0;

  Float_t Pt[NHistos];
  Float_t PtError[NHistos];
  Float_t Nmesons[NHistos];
  Float_t NmesonsError[NHistos];

  Float_t Ntota, Nback, Norma, NormaError, Renorma;

  for(iHistos=0; iHistos<NHistos; iHistos++)
    {
      signal     = h_Minv_pT->ProjectionX("signal",         NumberOfPtBins*iHistos+1,NumberOfPtBins*(iHistos+1));
      background = h_Minv_pT_back->ProjectionX("background",NumberOfPtBins*iHistos+1,NumberOfPtBins*(iHistos+1));
      
      ratio = new TH1D(*signal);
      ratio->Sumw2(); 
      ratio->Add(background,-1.0);
      ratio->Divide(background);
      difference = new TH1D(*signal);
      difference->Sumw2();
      ratio->Fit("pol0","","",Norma_minv_1,Norma_minv_2); 
      if (background->Integral(Norma_1,Norma_2) == 0)
	Renorma = 0.;
      else
	Renorma = signal->Integral(Norma_1,Norma_2)/background->Integral(Norma_1,Norma_2);
      difference->Add(background,(-1.)*Renorma); 
      

      Ntota = signal->Integral(Minv_1,Minv_2);
      Nback = background->Integral(Minv_1,Minv_2);
      Norma = ratio->GetFunction("pol0")->GetParameter(0);
      if (Renorma == 0.)
	NormaError = 0.;
      else
	NormaError =  ratio->GetFunction("pol0")->GetParError(0);
      printf("Ntotal %f Nback %f Norma %f and NormaError %f \n",Ntota, Nback, Norma, NormaError);
      printf("differencia is %f \n",difference->Integral(Minv_1,Minv_2));
      Nmesons[iHistos] = Ntota - Renorma * Nback;
      NmesonsError[iHistos] = TMath::Sqrt( Ntota + Nback*Renorma*Renorma + Nback*Nback*NormaError*NormaError );  
      Pt[iHistos] = (iHistos+0.5)*NumberOfPtBins*PtCalibration;
      PtError[iHistos] = NumberOfPtBins*PtCalibration/2.; 
      
    }

  char filenameout[80];
  sprintf(filenameout,"%s.PtSpectrum_%d_%s_%s",filename, NumberPerPtBin, particle, opt);
  TFile out(filenameout,"recreate");
  TGraphErrors * PtSpectrum = new TGraphErrors(NHistos, Pt, Nmesons, PtError, NmesonsError);
  PtSpectrum->SetName("PtSpectrum");
  PtSpectrum->Write();
  out.Close();

  frame->Draw();
  frame->SetStats(0);
  frame->SetXTitle("Neutral meson pT (GeV/c)");
  frame->SetYTitle("Number of neutral mesons per pT bin");
  PtSpectrum->SetMarkerStyle(27);
  PtSpectrum->Draw("P");


}
Commit	Line	Data
f09fd69e	1	//===============================================================
	2	// In this Macro (which can (must) be compilated), you will find all the
	3	// analysis functions to build photon spectrum, invariant mass
	4	// spectrum of photon pairs and combinatorial background calculations
	5	// in ALICE electromagnetic calorimeter
	6	// Author: Gines MARTINEZ, Subatech, 15 june 2001
	7	//==============================================================
	8	#include "TH2.h"
	9	#include "TH1.h"
a0531a78	10	#include "TClonesArray.h"
f09fd69e	11	#include "TFile.h"
079e348c	12	#include "TTree.h"
f09fd69e	13	#include "TRandom.h"
f09fd69e	14	#include "TObjectTable.h"
079e348c	15	#include "AliRun.h"
079e348c	16	#include "AliPHOSGetter.h"
f09fd69e	17	#include "AliPHOSRecParticle.h"
	18	#include "TLorentzVector.h"
	19	#include "TGraphErrors.h"
	20	#include "TF1.h"
	21
	22	TObjectTable * gObjectTable;
84ae9571	23	TRandom * gRandom;
84ae9571	24	AliRun * gAlice;
f09fd69e	25
f09fd69e	26
a0531a78	27	Bool_t PhotonId(AliPHOSRecParticle * RecParticle, char *opt)
	28	{
	29	Bool_t PhotonId = kTRUE ;
	30	if (strstr(opt,"narrow"))
	31	{
	32	if ( RecParticle->GetType() & 1) PhotonId= kFALSE ;
	33	}
	34	if (strstr(opt,"cpv"))
	35	{
	36	if ( (RecParticle->GetType() & 1<<2) >>2 ) PhotonId= kFALSE ;
	37	}
	38	if (strstr(opt,"time"))
	39	{
	40
	41	}
	42
	43	return PhotonId;
	44	}
	45
	46	void AnaMinv(char * filename,char * opt)
f09fd69e	47	{
	48	TH2F * h_Minv_lowpT = new TH2F("h_Minv_lowpT","Minv vs pT low",500,0.0,1.0,40,0.,10.);
	49	TH2F * h_Minv_highpT = new TH2F("h_Minv_highpT","Minv vs pT high",500,0.0,1.0,50,0.,100.);
	50	TH2F * h_Minv_lowpT_back = new TH2F("h_Minv_lowpT_back","Minv vs pT low back",500,0.0,1.0,40,0.,10.);
	51	TH2F * h_Minv_highpT_back = new TH2F("h_Minv_highpT_back","Minv vs pT high back",500,0.0,1.0,50,0.,100.);
	52
	53
	54	TH1F * h_Pseudoeta = new TH1F("h_Pseudoeta","Pseudoeta photons",500,-1.0,1.0);
	55	TH1F * h_Pt = new TH1F("h_Pt","Pt photons",400,0.,10.);
	56	TH2F * h_Peta_Pt = new TH2F("h_Peta_Pt","Pseudo vs pT",40,0.,10.,50,-1.0,1.0);
	57	TH1F * h_Phi = new TH1F("h_Phi","Phi photons",400,-4.,4.);
	58	TH2F * h_Peta_Phi = new TH2F("h_Peta_Phi","Pseudo vs Phi",200,-4,4,200,-1.0,1.0);
84ae9571	59	TH1F * h_Dispersion= new TH1F("h_Dispersion","Dispersion",50,0.,10.);
a0531a78	60	TH1F * h_Type = new TH1F("h_Type","Particle Type",9,-1.5,7.5);
f09fd69e	61
84ae9571	62	TH1F * h_DeltaR = new TH1F("h_DeltaR","Delta R",400,0.,2.);
f09fd69e	63	TH1F * h_Asymmetry= new TH1F("h_Asymmetry","Asymmetry",400, -2., 2.);
f09fd69e	64
079e348c	65	AliPHOSGetter * RecData = AliPHOSGetter::GetInstance(filename,"Gines") ;
f09fd69e	66
	67	AliPHOSRecParticle * RecParticle1;
	68	AliPHOSRecParticle * RecParticle2;
f09fd69e	69
84ae9571	70	Float_t RelativeRCut = 0.00001 ;
f09fd69e	71	Float_t AsymmetryCut = 0.7 ;
f09fd69e	72	Float_t Asymmetry;
84ae9571	73	Float_t Type;
f09fd69e	74
	75	Int_t iEvent, iRecParticle1, iRecParticle2;
	76	Int_t nRecParticle;
	77	Float_t invariant_mass, invariant_mass_mixed;
a0531a78	78	TClonesArray * RecParticles;
84ae9571	79
f09fd69e	80	Float_t average_multiplicity = 0.;
f09fd69e	81
079e348c	82	for(iEvent=0; iEvent<gAlice->TreeE()->GetEntries(); iEvent++)
f09fd69e	83	{
84ae9571	84	TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
84ae9571	85
a0531a78	86	RecData->Event(iEvent,"R");
f09fd69e	87	printf(">>> Event %d \n",iEvent);
a0531a78	88	RecParticles = RecData->RecParticles();
	89	nRecParticle = RecParticles->GetEntries();
	90	printf(">>> >> Recparticle multilicity is %d \n",nRecParticle);
079e348c	91	average_multiplicity += ((Float_t) (nRecParticle) ) / ( (Float_t)gAlice->TreeE()->GetEntries() ) ;
f09fd69e	92	// Construction de la masse invariante des pairs
	93	if (nRecParticle > 1)
	94	{
	95	for(iRecParticle1=0; iRecParticle1<nRecParticle; iRecParticle1++)
	96	{
a0531a78	97	RecParticle1 = (AliPHOSRecParticle *) RecParticles->At(iRecParticle1);
f09fd69e	98	RecParticle1->Momentum(P_photon1);
a0531a78	99	printf(">>> >> Photon momentum is %f \n", P_photon1.Pt() );
84ae9571	100	Type = RecParticle1->GetType();
84ae9571	101	h_Type->Fill(Type);
a0531a78	102	// Photon Id check
	103	if ( PhotonId(RecParticle1, opt) )
	104	{
	105
	106	h_Pseudoeta->Fill(P_photon1.PseudoRapidity());
	107	h_Pt->Fill(P_photon1.Pt());
	108	h_Phi->Fill(P_photon1.Phi());
	109	h_Peta_Pt->Fill(P_photon1.Pt(), P_photon1.PseudoRapidity());
	110	h_Peta_Phi->Fill(P_photon1.Phi(), P_photon1.PseudoRapidity() );
84ae9571	111
a0531a78	112	for(iRecParticle2=iRecParticle1+1; iRecParticle2<nRecParticle; iRecParticle2++)
	113	{
	114	RecParticle2 = (AliPHOSRecParticle *) RecParticles->At(iRecParticle2);
	115	RecParticle2->Momentum(P_photon2);
	116	Asymmetry = TMath::Abs((P_photon1.E()-P_photon2.E())/(P_photon1.E()+P_photon2.E()));
	117	//Photon Id Check
	118	if ( PhotonId(RecParticle2, opt) )
	119	{
	120	if ( (P_photon1 != P_photon2) &&
	121	(P_photon1.DeltaR(P_photon2) > RelativeRCut) &&
	122	(Asymmetry < AsymmetryCut) )
	123	{
	124	h_DeltaR->Fill(P_photon1.DeltaR(P_photon2));
	125	h_Asymmetry->Fill( Asymmetry );
	126
	127	// printf("A. p1 es %f \n",P_photon1->E());
	128	invariant_mass = (P_photon1 + P_photon2).M();
	129	// printf("B. p1 es %f \n",P_photon1->E());
	130	h_Minv_lowpT->Fill(invariant_mass, (P_photon1 + P_photon2).Pt() );
	131	h_Minv_highpT->Fill(invariant_mass,(P_photon1 + P_photon2).Pt() );
	132	}
	133	}
	134	}
f09fd69e	135	}
	136	}
	137	}
	138	}
	139	printf(">>> Average Multiplicity is %f \n",average_multiplicity);
079e348c	140	Int_t Background = (Int_t) (gAlice->TreeE()->GetEntries() * average_multiplicity * (average_multiplicity-1.)/2.) ;
f09fd69e	141	printf(">>> Background is %d \n",Background);
	142
	143	Double_t Pt_Mixed1, Pt_Mixed2;
	144	Double_t Y_Mixed1, Y_Mixed2;
	145	Double_t Phi_Mixed1, Phi_Mixed2;
	146
	147	for(iEvent=0; iEvent<Background; iEvent++)
	148	{
84ae9571	149	TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
a0531a78	150	//printf(">>> Background Event %d \n",iEvent);
f09fd69e	151	Pt_Mixed1 = h_Pt->GetRandom();
	152	Pt_Mixed2 = h_Pt->GetRandom();
	153	h_Peta_Phi->GetRandom2(Phi_Mixed1, Y_Mixed1);
	154	h_Peta_Phi->GetRandom2(Phi_Mixed2, Y_Mixed2);
	155	P_photonMixed1.SetPtEtaPhiM( Pt_Mixed1, Y_Mixed1, Phi_Mixed1, 0.0);
	156	P_photonMixed2.SetPtEtaPhiM( Pt_Mixed2, Y_Mixed2, Phi_Mixed2, 0.0);
	157	Asymmetry = TMath::Abs((P_photonMixed1.E()-P_photonMixed2.E())/(P_photonMixed1.E()+P_photonMixed2.E()));
	158
	159	if ( (P_photonMixed1.DeltaR(P_photonMixed2) > RelativeRCut) &&
	160	(Asymmetry < AsymmetryCut ) )
	161	{
	162	invariant_mass_mixed = (P_photonMixed1 + P_photonMixed2).M();
	163	h_Minv_lowpT_back->Fill(invariant_mass_mixed, (P_photonMixed1 + P_photonMixed2).Pt() );
	164	h_Minv_highpT_back->Fill(invariant_mass_mixed,(P_photonMixed1 + P_photonMixed2).Pt() );
	165	}
	166	}
a0531a78	167	printf (">>> Background calculation finished ! \n");
f09fd69e	168
f09fd69e	169	char outputname[80];
a0531a78	170	sprintf(outputname,"%s.Minv%s",filename,opt);
f09fd69e	171	TFile output(outputname,"recreate");
	172	h_Minv_lowpT->Write();
	173	h_Minv_highpT->Write();
	174	h_Minv_lowpT_back->Write();
	175	h_Minv_highpT_back->Write();
	176	h_Pseudoeta->Write();
	177	h_Pt->Write();
	178	h_Peta_Pt->Write();
	179	h_Phi->Write();
	180	h_Peta_Phi->Write();
84ae9571	181	h_Dispersion->Write();
84ae9571	182	h_Type->Write();
f09fd69e	183	h_Asymmetry->Write();
	184	h_DeltaR->Write();
	185
	186	output.Close();
	187	}
	188
	189
	190	void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, Option_t * opt)
	191	{
	192
	193	Int_t NumberOfPtBins = NumberPerPtBin;
	194	Float_t PtCalibration = 0.250;
	195
	196	TFile * in = new TFile(filename);
	197
	198	TH2F * h_Minv_pT = 0;
	199	TH2F * h_Minv_pT_back = 0;
	200	TH2F * frame = 0 ;
	201
	202	if (strstr(opt,"low"))
	203	{
	204	h_Minv_pT = (TH2F *) in->Get("h_Minv_lowpT"); ;
	205	h_Minv_pT_back = (TH2F *) in->Get("h_Minv_lowpT_back");
	206	PtCalibration = 0.250;
	207	frame = new TH2F("PtSpectrumlow","Pt Spectrum low",10, 0.,10.,10,0.1,10000);
	208	}
	209	if (strstr(opt,"high"))
	210	{
	211	h_Minv_pT = (TH2F *) in->Get("h_Minv_highpT"); ;
	212	h_Minv_pT_back = (TH2F *) in->Get("h_Minv_highpT_back");
	213	PtCalibration = 2.5;
	214	frame = new TH2F("PtSpectrumhigh","Pt Spectrum high",10, 0.,100.,10,0.1,10000);
	215	}
	216
	217	if ( h_Minv_pT == 0 )
	218	{
	219	printf(">>> Bad Option! \n");
	220	return;
	221	}
	222	Int_t Norma_1 = 100; Float_t Norma_minv_1 = 0.2;
	223	Int_t Norma_2 = 200; Float_t Norma_minv_2 = 0.4;
	224
84ae9571	225	Int_t Minv_1 = 46;
f09fd69e	226	Int_t Minv_2 = 76;
	227	if (strstr(particle,"eta"))
	228	{
84ae9571	229	Minv_1 = 214;
f09fd69e	230	Minv_2 = 314;
	231	}
	232
	233	if (strstr(particle,"norma"))
	234	{
	235	Minv_1 = 100;
	236	Minv_2 = 200;
	237	}
	238
	239	Int_t NHistos = 40/NumberOfPtBins;
	240	Int_t iHistos;
	241
	242	TH1D * signal = 0;
	243	TH1D * background = 0;
	244	TH1D * ratio = 0;
	245	TH1D * difference = 0;
	246
	247	Float_t Pt[NHistos];
	248	Float_t PtError[NHistos];
	249	Float_t Nmesons[NHistos];
	250	Float_t NmesonsError[NHistos];
	251
	252	Float_t Ntota, Nback, Norma, NormaError, Renorma;
	253
	254	for(iHistos=0; iHistos<NHistos; iHistos++)
	255	{
	256	signal = h_Minv_pT->ProjectionX("signal", NumberOfPtBinsiHistos+1,NumberOfPtBins(iHistos+1));
	257	background = h_Minv_pT_back->ProjectionX("background",NumberOfPtBinsiHistos+1,NumberOfPtBins(iHistos+1));
84ae9571	258
f09fd69e	259	ratio = new TH1D(*signal);
	260	ratio->Sumw2();
	261	ratio->Add(background,-1.0);
	262	ratio->Divide(background);
	263	difference = new TH1D(*signal);
	264	difference->Sumw2();
	265	ratio->Fit("pol0","","",Norma_minv_1,Norma_minv_2);
	266	if (background->Integral(Norma_1,Norma_2) == 0)
	267	Renorma = 0.;
	268	else
	269	Renorma = signal->Integral(Norma_1,Norma_2)/background->Integral(Norma_1,Norma_2);
	270	difference->Add(background,(-1.)*Renorma);
84ae9571	271
f09fd69e	272
	273	Ntota = signal->Integral(Minv_1,Minv_2);
	274	Nback = background->Integral(Minv_1,Minv_2);
	275	Norma = ratio->GetFunction("pol0")->GetParameter(0);
	276	if (Renorma == 0.)
	277	NormaError = 0.;
	278	else
	279	NormaError = ratio->GetFunction("pol0")->GetParError(0);
	280	printf("Ntotal %f Nback %f Norma %f and NormaError %f \n",Ntota, Nback, Norma, NormaError);
	281	printf("differencia is %f \n",difference->Integral(Minv_1,Minv_2));
	282	Nmesons[iHistos] = Ntota - Renorma * Nback;
	283	NmesonsError[iHistos] = TMath::Sqrt( Ntota + NbackRenormaRenorma + NbackNbackNormaError*NormaError );
	284	Pt[iHistos] = (iHistos+0.5)NumberOfPtBinsPtCalibration;
	285	PtError[iHistos] = NumberOfPtBins*PtCalibration/2.;
84ae9571	286
f09fd69e	287	}
f09fd69e	288
	289	char filenameout[80];
	290	sprintf(filenameout,"%s.PtSpectrum_%d_%s_%s",filename, NumberPerPtBin, particle, opt);
	291	TFile out(filenameout,"recreate");
	292	TGraphErrors * PtSpectrum = new TGraphErrors(NHistos, Pt, Nmesons, PtError, NmesonsError);
	293	PtSpectrum->SetName("PtSpectrum");
	294	PtSpectrum->Write();
	295	out.Close();
	296
	297	frame->Draw();
	298	frame->SetStats(0);
	299	frame->SetXTitle("Neutral meson pT (GeV/c)");
	300	frame->SetYTitle("Number of neutral mesons per pT bin");
	301	PtSpectrum->SetMarkerStyle(27);
	302	PtSpectrum->Draw("P");
	303
84ae9571	304
f09fd69e	305	}