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 "TClonesArray.h"
  11 #include "TFile.h"
  12 #include "TTree.h"
  13 #include "TRandom.h"
  14 #include "TObjectTable.h"
  15 #include "AliRun.h"
  16 #include "AliPHOSGetter.h"
  17 #include "AliPHOSRecParticle.h"
  18 #include "TLorentzVector.h"
  19 #include "TGraphErrors.h"
  20 #include "TF1.h"
  21
  22 TObjectTable * gObjectTable;
  23 TRandom      * gRandom;
  24 AliRun       * gAlice;
  25
  26
  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)
  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);
  59   TH1F * h_Dispersion= new TH1F("h_Dispersion","Dispersion",50,0.,10.);
  60   TH1F * h_Type      = new TH1F("h_Type","Particle Type",9,-1.5,7.5);
  61
  62   TH1F * h_DeltaR   = new TH1F("h_DeltaR","Delta R",400,0.,2.);
  63   TH1F * h_Asymmetry= new TH1F("h_Asymmetry","Asymmetry",400, -2., 2.);
  64
  65   AliPHOSGetter * RecData = AliPHOSGetter::GetInstance(filename,"Gines")  ;
  66
  67   AliPHOSRecParticle * RecParticle1;
  68   AliPHOSRecParticle * RecParticle2;
  69
  70   Float_t RelativeRCut = 0.00001 ;
  71   Float_t AsymmetryCut = 0.7 ;
  72   Float_t Asymmetry;
  73   Float_t Type;
  74
  75   Int_t iEvent, iRecParticle1, iRecParticle2;
  76   Int_t nRecParticle;
  77   Float_t invariant_mass, invariant_mass_mixed;
  78   TClonesArray * RecParticles;
  79
  80   Float_t average_multiplicity = 0.;
  81
  82   for(iEvent=0; iEvent<gAlice->TreeE()->GetEntries(); iEvent++)
  83     {
  84       TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
  85
  86       RecData->Event(iEvent,"R");
  87       printf(">>> Event %d \n",iEvent);
  88       RecParticles = RecData->RecParticles();
  89       nRecParticle = RecParticles->GetEntries();
  90       printf(">>> >> Recparticle multilicity is %d \n",nRecParticle);
  91       average_multiplicity += ((Float_t) (nRecParticle) ) / ( (Float_t)gAlice->TreeE()->GetEntries() ) ;
  92       // Construction de la masse invariante des pairs
  93       if (nRecParticle > 1)
  94         {
  95           for(iRecParticle1=0; iRecParticle1<nRecParticle; iRecParticle1++)
  96             {
  97               RecParticle1 =  (AliPHOSRecParticle *) RecParticles->At(iRecParticle1);
  98               RecParticle1->Momentum(P_photon1);
  99               printf(">>> >> Photon momentum is %f \n", P_photon1.Pt() );
 100               Type = RecParticle1->GetType();
 101               h_Type->Fill(Type);
 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() );
 111
 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                     }
 135                 }
 136             }
 137         }
 138     }
 139   printf(">>> Average Multiplicity is %f \n",average_multiplicity);
 140   Int_t Background = (Int_t)  (gAlice->TreeE()->GetEntries() * average_multiplicity * (average_multiplicity-1.)/2.) ;
 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     {
 149       TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
 150       //printf(">>> Background Event %d \n",iEvent);
 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     }
 167   printf (">>> Background calculation finished ! \n");
 168
 169   char outputname[80];
 170   sprintf(outputname,"%s.Minv%s",filename,opt);
 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();
 181   h_Dispersion->Write();
 182   h_Type->Write();
 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
 225   Int_t Minv_1 = 46;
 226   Int_t Minv_2 = 76;
 227   if (strstr(particle,"eta"))
 228     {
 229       Minv_1 = 214;
 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));
 258
 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);
 271
 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.;
 286
 287     }
 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
 304
 305 }