Inheritance from TObject. Automatic streamers.
[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"
a0531a78 10#include "TClonesArray.h"
f09fd69e 11#include "TFile.h"
079e348c 12#include "TTree.h"
f09fd69e 13#include "TRandom.h"
14#include "TObjectTable.h"
079e348c 15#include "AliRun.h"
16#include "AliPHOSGetter.h"
f09fd69e 17#include "AliPHOSRecParticle.h"
18#include "TLorentzVector.h"
19#include "TGraphErrors.h"
20#include "TF1.h"
21
22TObjectTable * gObjectTable;
84ae9571 23TRandom * gRandom;
24AliRun * gAlice;
f09fd69e 25
26
a0531a78 27Bool_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
46void 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.);
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 ;
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.;
81
079e348c 82 for(iEvent=0; iEvent<gAlice->TreeE()->GetEntries(); iEvent++)
f09fd69e 83 {
84ae9571 84 TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
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();
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
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();
182 h_Type->Write();
f09fd69e 183 h_Asymmetry->Write();
184 h_DeltaR->Write();
185
186 output.Close();
187}
188
189
190void 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", NumberOfPtBins*iHistos+1,NumberOfPtBins*(iHistos+1));
257 background = h_Minv_pT_back->ProjectionX("background",NumberOfPtBins*iHistos+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 + Nback*Renorma*Renorma + Nback*Nback*NormaError*NormaError );
284 Pt[iHistos] = (iHistos+0.5)*NumberOfPtBins*PtCalibration;
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}