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
  27 #include "TString.h"
  28 #include "TFile.h"
  29 #include "TLorentzVector.h"
  30 #include "TList.h"
  31 #include "TTree.h"
  32 #include "TParticle.h"
  33 #include "AliRun.h"
  34 #include "AliPHOSGammaJet.h"
  35 #include "AliPHOSGetter.h"
  36
  37 ClassImp(AliPHOSGammaJet)
  38
  39 //____________________________________________________________________________
  40 AliPHOSGammaJet::AliPHOSGammaJet() {
  41   // ctor
  42 }
  43
  44 //____________________________________________________________________________
  45 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
  46   // ctor
  47   AliPHOSGetter::Instance(inputfilename) ;
  48 }
  49
  50 //____________________________________________________________________________
  51 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) {
  52 }
  53
  54 //____________________________________________________________________________
  55 AliPHOSGammaJet::~AliPHOSGammaJet() {
  56 }
  57
  58 //____________________________________________________________________________
  59 void AliPHOSGammaJet::Exec(Option_t *)
  60 {
  61   // does the job
  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 id_g = -1;
 118     GetGammaJet(particleList,gamma, id_g);
 119     GetLeadingCharge(particleList,charge, id_g);
 120     GetLeadingPi0(particleList,pi0);
 121     Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), id_g) ;
 122     Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
 123     Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
 124     //    GetLeadingGammaPair(particleList, gammapair, id_g,
 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   TParticle *particle = 0x0;
 173
 174   Int_t iPrimary=-1, id_motherg = -1;
 175
 176   TIter next(&particleList);
 177   while ( (particle = (TParticle*)next()) ) {
 178     iPrimary++;
 179     Int_t ksCode = particle->GetStatusCode();
 180     Int_t iMother= particle->GetMother(0);
 181
 182     if (ksCode == 21 && iMother == -1)
 183       if(particle->GetPdgCode()==22){
 184         id_motherg = iPrimary;
 185         //cout<<"idmother "<<id_motherg<<endl;
 186       }
 187     if(ksCode == 1){
 188
 189       if(id_motherg == particle->GetMother(0)){
 190         particle->Momentum(gamma);
 191         id = iPrimary;
 192         break;
 193       }
 194     }// kscode == 1
 195   }// while
 196 }
 197
 198 //____________________________________________________________________________
 199 void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id)
 200 {
 201   TParticle *particle = 0x0;
 202
 203   Int_t iPrimary=-1;
 204   Double_t pt_max = 0, pt_i = 0;
 205   TIter next(&particleList);
 206   while ( (particle = (TParticle*)next()) ) {
 207     iPrimary++;
 208     Int_t ksCode = particle->GetStatusCode();
 209
 210     if((ksCode == 1)&&(id != iPrimary)
 211        &&(particle->GetPDG(0)->Charge()!=0)){
 212       pt_i = particle->Pt();
 213       if(pt_i> pt_max){
 214         pt_max = pt_i;
 215         particle->Momentum(charge);
 216       }//pt_max
 217     }// kscode == 1
 218   }// while
 219 }
 220
 221 //____________________________________________________________________________
 222 void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0)
 223 {
 224   TParticle *particle = 0x0;
 225
 226   Int_t iPrimary=-1;
 227   Double_t pt_max = 0, pt_i = 0;
 228   TIter next(&particleList);
 229   while ( (particle = (TParticle*)next()) ) {
 230     iPrimary++;
 231     Int_t ksCode = particle->GetStatusCode();
 232
 233     if((ksCode == 2))
 234       {
 235         pt_i = particle->Pt();
 236         if(pt_i> pt_max){
 237           pt_max = pt_i;
 238           particle->Momentum(pi0);
 239         }//pt_max
 240       }// kscode == 1
 241   }// while
 242 }
 243
 244 //____________________________________________________________________________
 245 //  void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id,
 246 //                       Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2,
 247 //                       Double_t & invmasscut1,Double_t & invmasscut2)
 248 //  {
 249 //   TParticle *particle = 0x0;
 250
 251 //   Int_t  iPrimary=-1;
 252 //   Double_t El = 0, E12 = 0;
 253 //   TLorentzVector gamma_i,gamma_j;
 254 //   TIter next(&particleList);
 255 //   while ( (particle = (TParticle*)next()) ) {
 256 //     iPrimary++;
 257 //     Int_t ksCode = particle->GetStatusCode();
 258 //     Int_t ksPdg = particle->GetPdgCode();
 259 //     if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
 260 //       particle->Momentum(gamma_i);
 261 //       Int_t jPrimary=-1;
 262 //       TIter next2(&particleList);
 263 //       while ( (particle = (TParticle*)next2()) ) {
 264 //      jPrimary++;
 265 //      if(jPrimary>iPrimary){
 266 //        ksCode = particle->GetStatusCode();
 267 //        ksPdg = particle->GetPdgCode();
 268 //        if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
 269 //          particle->Momentum(gamma_j);
 270 //          if(gamma_j.Angle(gamma_i.Vect())<thetacut){
 271 //            Float_t invmass = (gamma_i+gamma_j).M();
 272 //            h_invmass->Fill(Eg,invmass);
 273 //            if((invmass>invmasscut1) && (invmass<invmasscut2)){
 274 //              E12 =  (gamma_i+gamma_j).Energy();
 275 //              if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){
 276 //                //cout<<E12<<" "<<E12/Eg<<endl;
 277 //                El = E12;
 278 //                id_i = iPrimary;
 279 //                id_j = jPrimary;
 280 //                gammapair = gamma_i+gamma_j;
 281 //              }//E12>El && (E12/Eg>0.2 && E12/Eg<1.)
 282 //            }//(invmass>0.125) && (invmass<0.145)
 283 //          }//gamma_j.Angle(gamma_i.Vect())<0.04
 284 //        }//(ksCode == 1)
 285 //      }
 286 //       }//while
 287 //       //         cout<<"jPrimary "<<jPrimary<<endl;
 288 //     }// if kscode 1
 289 //   }//while
 290 //   // cout<<"iPrimary "<<iPrimary<<endl;
 291 //  }