]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AnaPhotons_v1.C
Typo corrected
[u/mrichter/AliRoot.git] / PHOS / AnaPhotons_v1.C
index 9b07afd7bdc9e9055d5a78d8581cd53ee8e875c0..1f82c39ea7dc602416cdf130295178e066ae9e22 100644 (file)
@@ -7,20 +7,43 @@
 //============================================================== 
 #include "TH2.h"
 #include "TH1.h"
+#include "TClonesArray.h"
 #include "TFile.h"
+#include "TTree.h"
 #include "TRandom.h"
 #include "TObjectTable.h"
-#include "AliPHOSIndexToObject.h"
+#include "AliRun.h"
+#include "AliPHOSGetter.h"
 #include "AliPHOSRecParticle.h"
 #include "TLorentzVector.h"
 #include "TGraphErrors.h"
 #include "TF1.h"
 
 TObjectTable * gObjectTable;
-TRandom * gRandom;
+TRandom      * gRandom;
+AliRun       * gAlice;
 
 
-void AnaMinv(char * filename)
+Bool_t PhotonId(AliPHOSRecParticle * RecParticle, char *opt)
+{
+  Bool_t PhotonId = kTRUE ;
+  if (strstr(opt,"narrow"))
+    {
+      if ( RecParticle->GetType() & 1) PhotonId= kFALSE ;
+    }
+  if (strstr(opt,"cpv"))
+    {
+      if (  (RecParticle->GetType() &  1<<2) >>2   ) PhotonId= kFALSE ;
+    }
+  if (strstr(opt,"time"))
+    {
+      
+    }
+
+  return PhotonId;
+}
+
+void AnaMinv(char * filename,char * opt)
 {
   TH2F * h_Minv_lowpT       = new TH2F("h_Minv_lowpT","Minv vs pT low",500,0.0,1.0,40,0.,10.);
   TH2F * h_Minv_highpT      = new TH2F("h_Minv_highpT","Minv vs pT high",500,0.0,1.0,50,0.,100.);
@@ -33,73 +56,88 @@ void AnaMinv(char * filename)
   TH2F * h_Peta_Pt   = new TH2F("h_Peta_Pt","Pseudo vs pT",40,0.,10.,50,-1.0,1.0);
   TH1F * h_Phi       = new TH1F("h_Phi","Phi photons",400,-4.,4.);
   TH2F * h_Peta_Phi  = new TH2F("h_Peta_Phi","Pseudo vs Phi",200,-4,4,200,-1.0,1.0);
+  TH1F * h_Dispersion= new TH1F("h_Dispersion","Dispersion",50,0.,10.);
+  TH1F * h_Type      = new TH1F("h_Type","Particle Type",9,-1.5,7.5);
 
-  TH1F * h_DeltaR   = new TH1F("h_DeltaR","Delta R",400,0,4);
+  TH1F * h_DeltaR   = new TH1F("h_DeltaR","Delta R",400,0.,2.);
   TH1F * h_Asymmetry= new TH1F("h_Asymmetry","Asymmetry",400, -2., 2.);
 
-  AliPHOSIndexToObject * RecData = AliPHOSIndexToObject::GetInstance(filename)  ;
+  AliPHOSGetter * RecData = AliPHOSGetter::GetInstance(filename,"Gines")  ;
  
   AliPHOSRecParticle * RecParticle1;
   AliPHOSRecParticle * RecParticle2;
 
-  Float_t RelativeRCut = 0.0001 ; 
+  Float_t RelativeRCut = 0.00001 ; 
   Float_t AsymmetryCut = 0.7 ;
   Float_t Asymmetry;
+  Float_t Type;
 
   Int_t iEvent, iRecParticle1, iRecParticle2;
   Int_t nRecParticle;
   Float_t invariant_mass, invariant_mass_mixed;
-  TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
+  TClonesArray * RecParticles;
   Float_t average_multiplicity = 0.;
 
-  for(iEvent=0; iEvent<RecData->GetMaxEvent(); iEvent++)
-    //  for(iEvent=0; iEvent<1000; iEvent++)
+  for(iEvent=0; iEvent<gAlice->TreeE()->GetEntries(); iEvent++)
     {
-      //     if (iEvent==2) gObjectTable->Print();
-      //if (iEvent==15) gObjectTable->Print();
-      RecData->GetEvent(iEvent);
+      TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
+
+      RecData->Event(iEvent,"R");
       printf(">>> Event %d \n",iEvent);
-      nRecParticle=RecData->GimeNRecParticles();
-      average_multiplicity += ((Float_t) (nRecParticle) ) / ( (Float_t)RecData->GetMaxEvent() ) ;
+      RecParticles = RecData->RecParticles();
+      nRecParticle = RecParticles->GetEntries();
+      printf(">>> >> Recparticle multilicity is %d \n",nRecParticle);
+      average_multiplicity += ((Float_t) (nRecParticle) ) / ( (Float_t)gAlice->TreeE()->GetEntries() ) ;
       // Construction de la masse invariante des pairs
       if (nRecParticle > 1) 
        {
          for(iRecParticle1=0; iRecParticle1<nRecParticle; iRecParticle1++)
            {
-             RecParticle1 = (AliPHOSRecParticle *)  RecData->GimeRecParticle(iRecParticle1);
+             RecParticle1 =  (AliPHOSRecParticle *) RecParticles->At(iRecParticle1);
              RecParticle1->Momentum(P_photon1);
-
-             h_Pseudoeta->Fill(P_photon1.PseudoRapidity());
-             h_Pt->Fill(P_photon1.Pt());
-             h_Phi->Fill(P_photon1.Phi());
-             h_Peta_Pt->Fill(P_photon1.Pt(), P_photon1.PseudoRapidity());
-             h_Peta_Phi->Fill(P_photon1.Phi(), P_photon1.PseudoRapidity() );
-        
-             for(iRecParticle2=iRecParticle1+1; iRecParticle2<nRecParticle; iRecParticle2++)
-               {
-                 RecParticle2 = (AliPHOSRecParticle *)  RecData->GimeRecParticle(iRecParticle2);
-                 RecParticle2->Momentum(P_photon2); 
-                 Asymmetry = TMath::Abs((P_photon1.E()-P_photon2.E())/(P_photon1.E()+P_photon2.E()));
-                 if ( (P_photon1 != P_photon2) && 
-                      (P_photon1.DeltaR(P_photon2) > RelativeRCut) &&
-                      (Asymmetry < AsymmetryCut)                          )
-                   {
-                     h_DeltaR->Fill(P_photon1.DeltaR(P_photon2));
-                     h_Asymmetry->Fill( Asymmetry );
-
-                     //   printf("A. p1 es %f \n",P_photon1->E());
-                     invariant_mass = (P_photon1 + P_photon2).M();
-                     // printf("B. p1 es %f \n",P_photon1->E());
-                     h_Minv_lowpT->Fill(invariant_mass, (P_photon1 + P_photon2).Pt() );
-                     h_Minv_highpT->Fill(invariant_mass,(P_photon1 + P_photon2).Pt() );
-                   }  
+             printf(">>> >> Photon momentum is %f \n", P_photon1.Pt() ); 
+             Type = RecParticle1->GetType();
+             h_Type->Fill(Type);
+             // Photon Id check
+             if ( PhotonId(RecParticle1, opt) )
+               {
+                 
+                 h_Pseudoeta->Fill(P_photon1.PseudoRapidity());
+                 h_Pt->Fill(P_photon1.Pt());
+                 h_Phi->Fill(P_photon1.Phi());
+                 h_Peta_Pt->Fill(P_photon1.Pt(), P_photon1.PseudoRapidity());
+                 h_Peta_Phi->Fill(P_photon1.Phi(), P_photon1.PseudoRapidity() );
+           
+                 for(iRecParticle2=iRecParticle1+1; iRecParticle2<nRecParticle; iRecParticle2++)
+                   {
+                     RecParticle2 = (AliPHOSRecParticle *) RecParticles->At(iRecParticle2);
+                     RecParticle2->Momentum(P_photon2); 
+                     Asymmetry = TMath::Abs((P_photon1.E()-P_photon2.E())/(P_photon1.E()+P_photon2.E()));
+                     //Photon Id Check
+                     if ( PhotonId(RecParticle2, opt) )
+                       {
+                         if ( (P_photon1 != P_photon2) && 
+                              (P_photon1.DeltaR(P_photon2) > RelativeRCut) &&
+                              (Asymmetry < AsymmetryCut)                          )
+                           {
+                             h_DeltaR->Fill(P_photon1.DeltaR(P_photon2));
+                             h_Asymmetry->Fill( Asymmetry );
+                             
+                             //   printf("A. p1 es %f \n",P_photon1->E());
+                             invariant_mass = (P_photon1 + P_photon2).M();
+                             // printf("B. p1 es %f \n",P_photon1->E());
+                             h_Minv_lowpT->Fill(invariant_mass, (P_photon1 + P_photon2).Pt() );
+                             h_Minv_highpT->Fill(invariant_mass,(P_photon1 + P_photon2).Pt() );
+                           } 
+                       }
+                   }
                }
            }
        }
     }
   printf(">>> Average Multiplicity is %f \n",average_multiplicity);
-  Int_t Background = (Int_t)  (RecData->GetMaxEvent() * average_multiplicity * (average_multiplicity-1.)/2.) ;
+  Int_t Background = (Int_t)  (gAlice->TreeE()->GetEntries() * average_multiplicity * (average_multiplicity-1.)/2.) ;
   printf(">>> Background is %d \n",Background);
 
   Double_t Pt_Mixed1, Pt_Mixed2;
@@ -108,7 +146,8 @@ void AnaMinv(char * filename)
 
   for(iEvent=0; iEvent<Background; iEvent++)
     {
-      //      printf(">>> Background Event %d \n",iEvent);
+      TLorentzVector P_photon1, P_photon2, P_photonMixed1, P_photonMixed2 ;
+      //printf(">>> Background Event %d \n",iEvent);
       Pt_Mixed1 =  h_Pt->GetRandom(); 
       Pt_Mixed2 =  h_Pt->GetRandom();
       h_Peta_Phi->GetRandom2(Phi_Mixed1, Y_Mixed1);
@@ -125,10 +164,10 @@ void AnaMinv(char * filename)
          h_Minv_highpT_back->Fill(invariant_mass_mixed,(P_photonMixed1 + P_photonMixed2).Pt() );
        }  
     }
-  
+  printf (">>> Background calculation finished ! \n");
 
   char outputname[80];
-  sprintf(outputname,"%s.Minv",filename);
+  sprintf(outputname,"%s.Minv%s",filename,opt);
   TFile output(outputname,"recreate");
   h_Minv_lowpT->Write();
   h_Minv_highpT->Write();
@@ -139,6 +178,8 @@ void AnaMinv(char * filename)
   h_Peta_Pt->Write();
   h_Phi->Write();
   h_Peta_Phi->Write();
+  h_Dispersion->Write();
+  h_Type->Write();
   h_Asymmetry->Write();
   h_DeltaR->Write();
   
@@ -181,11 +222,11 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O
   Int_t Norma_1 = 100; Float_t Norma_minv_1 = 0.2;
   Int_t Norma_2 = 200; Float_t Norma_minv_2 = 0.4;
 
-  Int_t Minv_1 = 56;
+  Int_t Minv_1 = 46;
   Int_t Minv_2 = 76;
   if (strstr(particle,"eta"))
     {
-      Minv_1 = 234;
+      Minv_1 = 214;
       Minv_2 = 314;
     }
 
@@ -214,8 +255,7 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O
     {
       signal     = h_Minv_pT->ProjectionX("signal",         NumberOfPtBins*iHistos+1,NumberOfPtBins*(iHistos+1));
       background = h_Minv_pT_back->ProjectionX("background",NumberOfPtBins*iHistos+1,NumberOfPtBins*(iHistos+1));
-      //signal->Rebin();
-      //background->Rebin();
+      
       ratio = new TH1D(*signal);
       ratio->Sumw2(); 
       ratio->Add(background,-1.0);
@@ -228,10 +268,7 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O
       else
        Renorma = signal->Integral(Norma_1,Norma_2)/background->Integral(Norma_1,Norma_2);
       difference->Add(background,(-1.)*Renorma); 
-       
-      //ratio->Draw();
-      //      background->Draw("same");
-      //     difference->Draw();
+      
 
       Ntota = signal->Integral(Minv_1,Minv_2);
       Nback = background->Integral(Minv_1,Minv_2);
@@ -246,11 +283,8 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O
       NmesonsError[iHistos] = TMath::Sqrt( Ntota + Nback*Renorma*Renorma + Nback*Nback*NormaError*NormaError );  
       Pt[iHistos] = (iHistos+0.5)*NumberOfPtBins*PtCalibration;
       PtError[iHistos] = NumberOfPtBins*PtCalibration/2.; 
-      //   ratio->Delete("");
-      //difference->Delete("");
+      
     }
-  // in->Close();
-
 
   char filenameout[80];
   sprintf(filenameout,"%s.PtSpectrum_%d_%s_%s",filename, NumberPerPtBin, particle, opt);
@@ -267,4 +301,5 @@ void AnaPtSpectrum(char * filename, Int_t NumberPerPtBin, Option_t * particle, O
   PtSpectrum->SetMarkerStyle(27);
   PtSpectrum->Draw("P");
 
+
 }