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