Macro compilble adapted to AliPHOSGetter
[u/mrichter/AliRoot.git] / PHOS / AnaPhotons_v1.C
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"
11 #include "TTree.h"
12 #include "TRandom.h"
13 #include "TObjectTable.h"
14 #include "AliRun.h"
15 #include "AliPHOSGetter.h"
16 #include "AliPHOSRecParticle.h"
17 #include "TLorentzVector.h"
18 #include "TGraphErrors.h"
19 #include "TF1.h"
20
21 TObjectTable * gObjectTable;
22 TRandom * gRandom;
23 AliRun * gAlice;
24
25
26 void 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
43   AliPHOSGetter * RecData = AliPHOSGetter::GetInstance(filename,"Gines")  ;
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
59   for(iEvent=0; iEvent<gAlice->TreeE()->GetEntries(); iEvent++)
60     //  for(iEvent=0; iEvent<1000; iEvent++)
61     {
62       //     if (iEvent==2) gObjectTable->Print();
63       //if (iEvent==15) gObjectTable->Print();
64       RecData->Event(iEvent);
65       printf(">>> Event %d \n",iEvent);
66       nRecParticle=RecData->NRecParticles();
67       average_multiplicity += ((Float_t) (nRecParticle) ) / ( (Float_t)gAlice->TreeE()->GetEntries() ) ;
68       // Construction de la masse invariante des pairs
69       if (nRecParticle > 1) 
70         {
71           for(iRecParticle1=0; iRecParticle1<nRecParticle; iRecParticle1++)
72             {
73               RecParticle1 = (AliPHOSRecParticle *)  RecData->RecParticle(iRecParticle1);
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                 {
84                   RecParticle2 = (AliPHOSRecParticle *)  RecData->RecParticle(iRecParticle2);
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);
105   Int_t Background = (Int_t)  (gAlice->TreeE()->GetEntries() * average_multiplicity * (average_multiplicity-1.)/2.) ;
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
152 void 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 }