Introducing new class : AliPHOSGammaJet

diff --git a/PHOS/AliPHOSGammaJet.cxx b/PHOS/AliPHOSGammaJet.cxx
new file mode 100644 (file)
index 0000000..26eb861
--- /dev/null
+++ b/PHOS/AliPHOSGammaJet.cxx
@@ -0,0 +1,292 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+/* $Id$ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma-jet correlations 
+//   
+// 
+//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN) 
+//////////////////////////////////////////////////////////////////////////////
+
+
+// --- ROOT system ---
+
+#include "TString.h"
+#include "TFile.h"
+#include "TLorentzVector.h"
+#include "TList.h"
+#include "TTree.h"
+#include "TParticle.h"
+#include "AliRun.h"
+#include "AliPHOSGammaJet.h" 
+#include "AliPHOSGetter.h" 
+
+ClassImp(AliPHOSGammaJet)
+
+//____________________________________________________________________________
+AliPHOSGammaJet::AliPHOSGammaJet() {
+  // ctor
+}
+
+//____________________________________________________________________________
+AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
+  // ctor 
+  AliPHOSGetter::Instance(inputfilename) ; 
+}
+
+//____________________________________________________________________________
+AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet * gj) {
+}
+
+//____________________________________________________________________________
+AliPHOSGammaJet::~AliPHOSGammaJet() {
+}
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::Exec(Option_t *option) 
+{
+  // does the job
+  
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  Int_t nEvent1 = gime->MaxEvent() ;   
+  Int_t iEvent = 0 ; 
+  for ( iEvent = 0 ; iEvent < nEvent1 ; iEvent++) {
+    if (iEvent <= 100 || iEvent%10 == 0)
+      Info("Exec", "Event %d", iEvent) ;     
+    TParticle * particle = 0 ;
+    //-----------------Fill list with particles--------------------
+    TLorentzVector pPi0, pGamma1, pGamma2 ;
+    Double_t angle = 0 ;
+    TList particleList ;
+    Int_t n = -1; 
+    gime->Event(iEvent, "X") ; 
+    Int_t  nparticles = gime->NPrimaries() ; 
+    Int_t iParticle=0 ;
+    for (iParticle=0 ; iParticle < nparticles ; iParticle++) {
+      //Keep original partons
+      particle = gime->Primary(iParticle) ; 
+      if((particle->GetStatusCode()== 21)){
+       particleList.Add(particle);
+       n++;
+      }
+      //Keep Stable particles within eta range
+      Float_t etacut = 0. ; 
+      if((particle->GetStatusCode() == 1)&&
+        (TMath::Abs(particle->Eta())<etacut)){
+       // Keep particles different from Pi0
+       if(particle->GetPdgCode() != 111){
+         particleList.Add(particle);
+         n++;
+       } 
+       //Decay Pi0 and keep it with different status name
+       //Keep decay photons
+       if(particle->GetPdgCode() == 111) {
+         //cout<<"Pi0 "<<endl;
+         n += 3 ; 
+         particle->Momentum(pPi0);
+         Pi0Decay(particle->GetMass(),pPi0,pGamma1,pGamma2,angle);     
+         TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
+                                             pGamma1.Pz(),pGamma1.E(),0,0,0,0);
+         TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(),
+                                             pGamma2.Pz(),pGamma2.E(),0,0,0,0);
+         photon1->SetWeight(1);
+         photon2->SetWeight(2);
+         particle->SetStatusCode(2);
+         particleList.Add(particle);
+         particleList.Add(photon1);
+         particleList.Add(photon2);
+         //photon1->Print();
+         //photon2->Print();  
+       }//if pi0
+      }//final particle etacut
+    }//for (iParticle<nParticle)
+    TLorentzVector gamma,charge,pi0, gammapair;
+    Int_t id_g = -1;
+    GetGammaJet(particleList,gamma, id_g);
+    GetLeadingCharge(particleList,charge, id_g);
+    GetLeadingPi0(particleList,pi0, id_g);
+    Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), id_g) ;
+    Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
+    Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
+    //    GetLeadingGammaPair(particleList, gammapair, id_g, 
+    //                 0.04,0.2, 1.0,0.13,0.14);
+    Info("Pi0Decay", "Pair: %f", gammapair.Energy()) ;
+  }//loop: events
+}    
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0, 
+                              TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
+  // Perform isotropic decay pi0 -> 2 photons
+  // p0 is pi0 4-momentum (input)
+  // p1 and p2 are photon 4-momenta (output)
+  //  cout<<"Boost vector"<<endl;
+  TVector3 b = p0.BoostVector();
+  //cout<<"Parameters"<<endl;
+  //Double_t mPi0   = p0.M();
+  Double_t phi    = TMath::TwoPi() * gRandom->Rndm();
+  Double_t cosThe = 2 * gRandom->Rndm() - 1;
+  Double_t cosPhi = TMath::Cos(phi);
+  Double_t sinPhi = TMath::Sin(phi);
+  Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
+  Double_t ePi0   = mPi0/2.;
+  //cout<<"ePi0 "<<ePi0<<endl;
+  //cout<<"Components"<<endl;
+  p1.SetPx(+ePi0*cosPhi*sinThe);
+  p1.SetPy(+ePi0*sinPhi*sinThe);
+  p1.SetPz(+ePi0*cosThe);
+  p1.SetE(ePi0);
+  //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
+  //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
+  p2.SetPx(-ePi0*cosPhi*sinThe);
+  p2.SetPy(-ePi0*sinPhi*sinThe);
+  p2.SetPz(-ePi0*cosThe);
+  p2.SetE(ePi0);
+  //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
+  //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
+  //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
+  p1.Boost(b);
+  //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
+  p2.Boost(b);
+  //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
+  //cout<<"angle"<<endl;
+  angle = p1.Angle(p2.Vect());
+  //cout<<angle<<endl;
+}
+//____________________________________________________________________________
+void AliPHOSGammaJet::GetGammaJet(TList &particleList, TLorentzVector &gamma, Int_t & id) 
+{
+  TParticle *particle = 0x0;
+  
+  Int_t iPrimary=-1, id_motherg = -1;
+  
+  TIter next(&particleList);
+  while ( (particle = (TParticle*)next()) ) {
+    iPrimary++;  
+    Int_t ksCode = particle->GetStatusCode();
+    Int_t iMother= particle->GetMother(0);
+    
+    if (ksCode == 21 && iMother == -1)
+      if(particle->GetPdgCode()==22){
+       id_motherg = iPrimary;
+       //cout<<"idmother "<<id_motherg<<endl;
+      }
+    if(ksCode == 1){
+      
+      if(id_motherg == particle->GetMother(0)){
+       particle->Momentum(gamma);
+       id = iPrimary;
+       break;
+      }
+    }// kscode == 1
+  }// while
+}
+ 
+//____________________________________________________________________________
+void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id) 
+{
+  TParticle *particle = 0x0;
+  
+  Int_t iPrimary=-1;
+  Double_t pt_max = 0, pt_i = 0;
+  TIter next(&particleList);
+  while ( (particle = (TParticle*)next()) ) {
+    iPrimary++;  
+    Int_t ksCode = particle->GetStatusCode();
+    
+    if((ksCode == 1)&&(id != iPrimary)
+       &&(particle->GetPDG(0)->Charge()!=0)){
+      pt_i = particle->Pt(); 
+      if(pt_i> pt_max){
+       pt_max = pt_i;
+       particle->Momentum(charge);
+      }//pt_max   
+    }// kscode == 1
+  }// while
+}
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0, Int_t & id) 
+{
+  TParticle *particle = 0x0;
+  
+  Int_t iPrimary=-1;
+  Double_t pt_max = 0, pt_i = 0;
+  TIter next(&particleList);
+  while ( (particle = (TParticle*)next()) ) {
+    iPrimary++;  
+    Int_t ksCode = particle->GetStatusCode();
+    
+    if((ksCode == 2))
+      //&&(id != iPrimary)&&(particle->GetPdgCode() == 111))
+      {
+       pt_i = particle->Pt(); 
+       if(pt_i> pt_max){
+         pt_max = pt_i;
+         particle->Momentum(pi0);
+       }//pt_max   
+      }// kscode == 1
+  }// while
+}
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id, 
+                        Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2,
+                        Double_t & invmasscut1,Double_t & invmasscut2) 
+{
+//   TParticle *particle = 0x0;
+  
+//   Int_t  iPrimary=-1;
+//   Double_t El = 0, E12 = 0;
+//   TLorentzVector gamma_i,gamma_j;
+//   TIter next(&particleList);
+//   while ( (particle = (TParticle*)next()) ) {
+//     iPrimary++;       
+//     Int_t ksCode = particle->GetStatusCode();
+//     Int_t ksPdg = particle->GetPdgCode();
+//     if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
+//       particle->Momentum(gamma_i);
+//       Int_t jPrimary=-1;
+//       TIter next2(&particleList);
+//       while ( (particle = (TParticle*)next2()) ) {
+//     jPrimary++;
+//     if(jPrimary>iPrimary){
+//       ksCode = particle->GetStatusCode();
+//       ksPdg = particle->GetPdgCode();
+//       if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
+//         particle->Momentum(gamma_j);
+//         if(gamma_j.Angle(gamma_i.Vect())<thetacut){
+//           Float_t invmass = (gamma_i+gamma_j).M();
+//           h_invmass->Fill(Eg,invmass);
+//           if((invmass>invmasscut1) && (invmass<invmasscut2)){
+//             E12 =  (gamma_i+gamma_j).Energy(); 
+//             if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){
+//               //cout<<E12<<" "<<E12/Eg<<endl;
+//               El = E12;
+//               id_i = iPrimary;
+//               id_j = jPrimary;
+//               gammapair = gamma_i+gamma_j;
+//             }//E12>El && (E12/Eg>0.2 && E12/Eg<1.)
+//           }//(invmass>0.125) && (invmass<0.145)
+//         }//gamma_j.Angle(gamma_i.Vect())<0.04
+//       }//(ksCode == 1)
+//     }
+//       }//while
+//       //        cout<<"jPrimary "<<jPrimary<<endl;
+//     }// if kscode 1
+//   }//while
+//   //        cout<<"iPrimary "<<iPrimary<<endl;
+}

diff --git a/PHOS/AliPHOSGammaJet.h b/PHOS/AliPHOSGammaJet.h
new file mode 100644 (file)
index 0000000..d0fd952
--- /dev/null
+++ b/PHOS/AliPHOSGammaJet.h
@@ -0,0 +1,38 @@
+#ifndef ALIPHOSGammaJet_H
+#define ALIPHOSGammaJet_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+/* $Id$ */
+
+//_________________________________________________________________________
+//  Class for the analysis of gamma-jet correlations     
+//                  
+//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
+
+// --- ROOT system ---
+#include "TTask.h"
+
+// --- AliRoot header files ---
+
+class AliPHOSGammaJet : public TTask {
+
+public: 
+
+  AliPHOSGammaJet() ; // default ctor
+  AliPHOSGammaJet(const TString inputfilename) ; //ctor 
+  AliPHOSGammaJet(const AliPHOSGammaJet * gj) ; // cpy ctor
+  ~AliPHOSGammaJet() ; // dtor
+  virtual void   Exec(Option_t *option); 
+  void GetGammaJet(TList & particleList, TLorentzVector & gamma, Int_t & id) ; 
+  void GetLeadingCharge(TList & particleList, TLorentzVector & charge, Int_t & id) ;
+  void GetLeadingPi0(TList & particleList, TLorentzVector & pi0, Int_t & id) ;
+  void GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id, 
+                          Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2,
+                          Double_t & invmasscut1,Double_t & invmasscut2);
+  void Pi0Decay(Double_t mPi0, TLorentzVector &p0, 
+               TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) ; 
+private: 
+
+  ClassDef(AliPHOSGammaJet,1)
+} ; 
+#endif //ALIPHOSGammaJet_H

diff --git a/PHOS/PHOSLinkDef.h b/PHOS/PHOSLinkDef.h
index 9837c867b6d727065a10e91261c72bbb72bd59b4..79b9a9bb51e172adfc14d097b78b915df655cc74 100644 (file)
--- a/PHOS/PHOSLinkDef.h
+++ b/PHOS/PHOSLinkDef.h
@@ -63,5 +63,6 @@
 #pragma link C++ class AliPHOSFastGlobalReconstruction+;
 #pragma link C++ class AliPHOSJet+;
 #pragma link C++ class AliPHOSJetFinder+;
+#pragma link C++ class AliPHOSGammaJet+;
 
 #endif

diff --git a/PHOS/libPHOS.pkg b/PHOS/libPHOS.pkg
index b877c4d67d17b4645472b90f221d7c8c163bd370..f7f3e80c3f9aadeb54322522083d037c0d368ffb 100644 (file)
--- a/PHOS/libPHOS.pkg
+++ b/PHOS/libPHOS.pkg
@@ -23,7 +23,8 @@ SRCS          =    AliPHOS.cxx AliPHOSv0.cxx AliPHOSv1.cxx AliPHOSv2.cxx \
                  AliPHOSRecManager.cxx AliPHOSRecCpvManager.cxx AliPHOSRecEmcManager.cxx \
                  AliPHOSClusterizerv2.cxx AliPHOSPIDv1.cxx AliPHOSJet.cxx AliPHOSJetFinder.cxx\
                  AliPHOSLoader.cxx AliPHOSRaw2Digits.cxx AliPHOSBeamTestEvent.cxx AliPHOSCalibrationDB.cxx AliPHOSCalibrationData.cxx\
-                 AliPHOSCalibrator.cxx AliPHOSCalibrManager.cxx AliPHOSConTableDB.cxx AliPHOSFastGlobalReconstruction.cxx AliPHOSGetter.cxx 
+                 AliPHOSCalibrator.cxx AliPHOSCalibrManager.cxx AliPHOSConTableDB.cxx AliPHOSFastGlobalReconstruction.cxx AliPHOSGetter.cxx \
+                AliPHOSGammaJet.cxx    
 
 HDRS:= $(SRCS:.cxx=.h)