]>
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" |
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 | ||
22 | TObjectTable * gObjectTable; | |
84ae9571 | 23 | TRandom * gRandom; |
24 | AliRun * gAlice; | |
f09fd69e | 25 | |
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.); |
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 | ||
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", 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 | } |