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