X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAnaPhotons_v1.C;h=1f82c39ea7dc602416cdf130295178e066ae9e22;hb=c2140715d75ba09019bc2373b20c7e2e077062c0;hp=50b9d35f5c94f32c65775d85d34a914818e30014;hpb=079e348c27ec78e2953e8eab94a6561b6c428b35;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AnaPhotons_v1.C b/PHOS/AnaPhotons_v1.C index 50b9d35f5c9..1f82c39ea7d 100644 --- a/PHOS/AnaPhotons_v1.C +++ b/PHOS/AnaPhotons_v1.C @@ -7,6 +7,7 @@ //============================================================== #include "TH2.h" #include "TH1.h" +#include "TClonesArray.h" #include "TFile.h" #include "TTree.h" #include "TRandom.h" @@ -19,11 +20,30 @@ #include "TF1.h" TObjectTable * gObjectTable; -TRandom * gRandom; -AliRun * gAlice; +TRandom * gRandom; +AliRun * gAlice; -void AnaMinv(char * filename) +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.); @@ -36,67 +56,82 @@ void AnaMinv(char * filename) 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,4); + 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.0001 ; + 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; - TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ; + TClonesArray * RecParticles; + Float_t average_multiplicity = 0.; for(iEvent=0; iEventTreeE()->GetEntries(); iEvent++) - // for(iEvent=0; iEvent<1000; iEvent++) { - // if (iEvent==2) gObjectTable->Print(); - //if (iEvent==15) gObjectTable->Print(); - RecData->Event(iEvent); + TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ; + + RecData->Event(iEvent,"R"); printf(">>> Event %d \n",iEvent); - nRecParticle=RecData->NRecParticles(); + 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; iRecParticle1RecParticle(iRecParticle1); + RecParticle1 = (AliPHOSRecParticle *) RecParticles->At(iRecParticle1); RecParticle1->Momentum(P_photon1); - - 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; iRecParticle2RecParticle(iRecParticle2); - RecParticle2->Momentum(P_photon2); - Asymmetry = TMath::Abs((P_photon1.E()-P_photon2.E())/(P_photon1.E()+P_photon2.E())); - 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(">>> >> 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; iRecParticle2At(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() ); + } + } + } } } } @@ -111,7 +146,8 @@ void AnaMinv(char * filename) for(iEvent=0; iEvent>> Background Event %d \n",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); @@ -128,10 +164,10 @@ void AnaMinv(char * filename) 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",filename); + sprintf(outputname,"%s.Minv%s",filename,opt); TFile output(outputname,"recreate"); h_Minv_lowpT->Write(); h_Minv_highpT->Write(); @@ -142,6 +178,8 @@ void AnaMinv(char * filename) h_Peta_Pt->Write(); h_Phi->Write(); h_Peta_Phi->Write(); + h_Dispersion->Write(); + h_Type->Write(); h_Asymmetry->Write(); h_DeltaR->Write(); @@ -184,11 +222,11 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O 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 = 56; + Int_t Minv_1 = 46; Int_t Minv_2 = 76; if (strstr(particle,"eta")) { - Minv_1 = 234; + Minv_1 = 214; Minv_2 = 314; } @@ -217,8 +255,7 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O { 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)); - //signal->Rebin(); - //background->Rebin(); + ratio = new TH1D(*signal); ratio->Sumw2(); ratio->Add(background,-1.0); @@ -231,10 +268,7 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O else Renorma = signal->Integral(Norma_1,Norma_2)/background->Integral(Norma_1,Norma_2); difference->Add(background,(-1.)*Renorma); - - //ratio->Draw(); - // background->Draw("same"); - // difference->Draw(); + Ntota = signal->Integral(Minv_1,Minv_2); Nback = background->Integral(Minv_1,Minv_2); @@ -249,11 +283,8 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O NmesonsError[iHistos] = TMath::Sqrt( Ntota + Nback*Renorma*Renorma + Nback*Nback*NormaError*NormaError ); Pt[iHistos] = (iHistos+0.5)*NumberOfPtBins*PtCalibration; PtError[iHistos] = NumberOfPtBins*PtCalibration/2.; - // ratio->Delete(""); - //difference->Delete(""); + } - // in->Close(); - char filenameout[80]; sprintf(filenameout,"%s.PtSpectrum_%d_%s_%s",filename, NumberPerPtBin, particle, opt); @@ -270,4 +301,5 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O PtSpectrum->SetMarkerStyle(27); PtSpectrum->Draw("P"); + }