PHOS/AliPHOSGammaJet.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16
  17 //_________________________________________________________________________
  18 // Class for the analysis of gamma-jet correlations
  19 //
  20 //
  21 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
  22 //////////////////////////////////////////////////////////////////////////////
  23
  24
  25 // --- ROOT system ---
  26 #include "TRandom.h"
  27
  28 #include "TLorentzVector.h"
  29 #include "TList.h"
  30 #include "TParticle.h"
  31 #include "AliPHOSGammaJet.h"
  32 #include "AliPHOSGetter.h"
  33
  34 ClassImp(AliPHOSGammaJet)
  35
  36 //____________________________________________________________________________
  37 AliPHOSGammaJet::AliPHOSGammaJet() {
  38   // ctor
  39 }
  40
  41 //____________________________________________________________________________
  42 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
  43   // ctor
  44   AliPHOSGetter::Instance(inputfilename) ;
  45 }
  46
  47 //____________________________________________________________________________
  48 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) {
  49 }
  50
  51 //____________________________________________________________________________
  52 AliPHOSGammaJet::~AliPHOSGammaJet() {
  53 }
  54
  55 //____________________________________________________________________________
  56 void AliPHOSGammaJet::Exec(Option_t * opt)
  57 {
  58   // does the job
  59
  60   if (! opt)
  61     return ;
  62
  63   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
  64   Int_t nEvent1 = gime->MaxEvent() ;
  65   Int_t iEvent = 0 ;
  66   for ( iEvent = 0 ; iEvent < nEvent1 ; iEvent++) {
  67     if (iEvent <= 100 || iEvent%10 == 0)
  68       Info("Exec", "Event %d", iEvent) ;
  69     TParticle * particle = 0 ;
  70     //-----------------Fill list with particles--------------------
  71     TLorentzVector pPi0, pGamma1, pGamma2 ;
  72     Double_t angle = 0 ;
  73     TList particleList ;
  74     Int_t n = -1;
  75     gime->Event(iEvent, "X") ;
  76     Int_t  nparticles = gime->NPrimaries() ;
  77     Int_t iParticle=0 ;
  78     for (iParticle=0 ; iParticle < nparticles ; iParticle++) {
  79       //Keep original partons
  80       particle = gime->Primary(iParticle) ;
  81       if((particle->GetStatusCode()== 21)){
  82         particleList.Add(particle);
  83         n++;
  84       }
  85       //Keep Stable particles within eta range
  86       Float_t etacut = 0. ;
  87       if((particle->GetStatusCode() == 1)&&
  88          (TMath::Abs(particle->Eta())<etacut)){
  89         // Keep particles different from Pi0
  90         if(particle->GetPdgCode() != 111){
  91           particleList.Add(particle);
  92           n++;
  93         }
  94         //Decay Pi0 and keep it with different status name
  95         //Keep decay photons
  96         if(particle->GetPdgCode() == 111) {
  97           //cout<<"Pi0 "<<endl;
  98           n += 3 ;
  99           particle->Momentum(pPi0);
 100           Pi0Decay(particle->GetMass(),pPi0,pGamma1,pGamma2,angle);
 101           TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
 102                                               pGamma1.Pz(),pGamma1.E(),0,0,0,0);
 103           TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(),
 104                                               pGamma2.Pz(),pGamma2.E(),0,0,0,0);
 105           photon1->SetWeight(1);
 106           photon2->SetWeight(2);
 107           particle->SetStatusCode(2);
 108           particleList.Add(particle);
 109           particleList.Add(photon1);
 110           particleList.Add(photon2);
 111           //photon1->Print();
 112           //photon2->Print();
 113         }//if pi0
 114       }//final particle etacut
 115     }//for (iParticle<nParticle)
 116     TLorentzVector gamma,charge,pi0, gammapair;
 117     Int_t idg = -1;
 118     GetGammaJet(particleList,gamma, idg);
 119     GetLeadingCharge(particleList,charge, idg);
 120     GetLeadingPi0(particleList,pi0);
 121     Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), idg) ;
 122     Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
 123     Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
 124     //    GetLeadingGammaPair(particleList, gammapair, idg,
 125     //                  0.04,0.2, 1.0,0.13,0.14);
 126     Info("Pi0Decay", "Pair: %f", gammapair.Energy()) ;
 127   }//loop: events
 128 }
 129
 130 //____________________________________________________________________________
 131 void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
 132                                TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
 133   // Perform isotropic decay pi0 -> 2 photons
 134   // p0 is pi0 4-momentum (input)
 135   // p1 and p2 are photon 4-momenta (output)
 136   //  cout<<"Boost vector"<<endl;
 137   TVector3 b = p0.BoostVector();
 138   //cout<<"Parameters"<<endl;
 139   //Double_t mPi0   = p0.M();
 140   Double_t phi    = TMath::TwoPi() * gRandom->Rndm();
 141   Double_t cosThe = 2 * gRandom->Rndm() - 1;
 142   Double_t cosPhi = TMath::Cos(phi);
 143   Double_t sinPhi = TMath::Sin(phi);
 144   Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
 145   Double_t ePi0   = mPi0/2.;
 146   //cout<<"ePi0 "<<ePi0<<endl;
 147   //cout<<"Components"<<endl;
 148   p1.SetPx(+ePi0*cosPhi*sinThe);
 149   p1.SetPy(+ePi0*sinPhi*sinThe);
 150   p1.SetPz(+ePi0*cosThe);
 151   p1.SetE(ePi0);
 152   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
 153   //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
 154   p2.SetPx(-ePi0*cosPhi*sinThe);
 155   p2.SetPy(-ePi0*sinPhi*sinThe);
 156   p2.SetPz(-ePi0*cosThe);
 157   p2.SetE(ePi0);
 158   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
 159   //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
 160   //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
 161   p1.Boost(b);
 162   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
 163   p2.Boost(b);
 164   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
 165   //cout<<"angle"<<endl;
 166   angle = p1.Angle(p2.Vect());
 167   //cout<<angle<<endl;
 168 }
 169 //____________________________________________________________________________
 170 void AliPHOSGammaJet::GetGammaJet(TList &particleList, TLorentzVector &gamma, Int_t & id)
 171 {
 172   // Get the lists of jet particles and gamma
 173   TParticle *particle = 0x0;
 174
 175   Int_t iPrimary=-1, id_motherg = -1;
 176
 177   TIter next(&particleList);
 178   while ( (particle = (TParticle*)next()) ) {
 179     iPrimary++;
 180     Int_t ksCode = particle->GetStatusCode();
 181     Int_t iMother= particle->GetMother(0);
 182
 183     if (ksCode == 21 && iMother == -1)
 184       if(particle->GetPdgCode()==22){
 185         id_motherg = iPrimary;
 186         //cout<<"idmother "<<id_motherg<<endl;
 187       }
 188     if(ksCode == 1){
 189
 190       if(id_motherg == particle->GetMother(0)){
 191         particle->Momentum(gamma);
 192         id = iPrimary;
 193         break;
 194       }
 195     }// kscode == 1
 196   }// while
 197 }
 198
 199 //____________________________________________________________________________
 200 void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id)
 201 {
 202   // Gets the leading particle from the list of charged particles
 203   TParticle *particle = 0x0;
 204
 205   Int_t iPrimary=-1;
 206   Double_t ptmax = 0, pti = 0;
 207   TIter next(&particleList);
 208   while ( (particle = (TParticle*)next()) ) {
 209     iPrimary++;
 210     Int_t ksCode = particle->GetStatusCode();
 211
 212     if((ksCode == 1)&&(id != iPrimary)
 213        &&(particle->GetPDG(0)->Charge()!=0)){
 214       pti = particle->Pt();
 215       if(pti> ptmax){
 216         ptmax = pti;
 217         particle->Momentum(charge);
 218       }//ptmax
 219     }// kscode == 1
 220   }// while
 221 }
 222
 223 //____________________________________________________________________________
 224 void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0)
 225 {
 226   // Gets the leading pi0 from the list of particles
 227   TParticle *particle = 0x0;
 228
 229   Int_t iPrimary=-1;
 230   Double_t ptmax = 0, pti = 0;
 231   TIter next(&particleList);
 232   while ( (particle = (TParticle*)next()) ) {
 233     iPrimary++;
 234     Int_t ksCode = particle->GetStatusCode();
 235
 236     if((ksCode == 2))
 237       {
 238         pti = particle->Pt();
 239         if(pti> ptmax){
 240           ptmax = pti;
 241           particle->Momentum(pi0);
 242         }//ptmax
 243       }// kscode == 1
 244   }// while
 245 }
 246
 247 //____________________________________________________________________________
 248 //  void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id,
 249 //                       Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2,
 250 //                       Double_t & invmasscut1,Double_t & invmasscut2)
 251 //  {
 252 //   TParticle *particle = 0x0;
 253
 254 //   Int_t  iPrimary=-1;
 255 //   Double_t El = 0, E12 = 0;
 256 //   TLorentzVector gamma_i,gamma_j;
 257 //   TIter next(&particleList);
 258 //   while ( (particle = (TParticle*)next()) ) {
 259 //     iPrimary++;
 260 //     Int_t ksCode = particle->GetStatusCode();
 261 //     Int_t ksPdg = particle->GetPdgCode();
 262 //     if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
 263 //       particle->Momentum(gamma_i);
 264 //       Int_t jPrimary=-1;
 265 //       TIter next2(&particleList);
 266 //       while ( (particle = (TParticle*)next2()) ) {
 267 //      jPrimary++;
 268 //      if(jPrimary>iPrimary){
 269 //        ksCode = particle->GetStatusCode();
 270 //        ksPdg = particle->GetPdgCode();
 271 //        if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
 272 //          particle->Momentum(gamma_j);
 273 //          if(gamma_j.Angle(gamma_i.Vect())<thetacut){
 274 //            Float_t invmass = (gamma_i+gamma_j).M();
 275 //            h_invmass->Fill(Eg,invmass);
 276 //            if((invmass>invmasscut1) && (invmass<invmasscut2)){
 277 //              E12 =  (gamma_i+gamma_j).Energy();
 278 //              if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){
 279 //                //cout<<E12<<" "<<E12/Eg<<endl;
 280 //                El = E12;
 281 //                id_i = iPrimary;
 282 //                id_j = jPrimary;
 283 //                gammapair = gamma_i+gamma_j;
 284 //              }//E12>El && (E12/Eg>0.2 && E12/Eg<1.)
 285 //            }//(invmass>0.125) && (invmass<0.145)
 286 //          }//gamma_j.Angle(gamma_i.Vect())<0.04
 287 //        }//(ksCode == 1)
 288 //      }
 289 //       }//while
 290 //       //         cout<<"jPrimary "<<jPrimary<<endl;
 291 //     }// if kscode 1
 292 //   }//while
 293 //   // cout<<"iPrimary "<<iPrimary<<endl;
 294 //  }