X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG4%2FAliGammaMCReader.cxx;h=d606dacfb8d8ab83512f36ea177a71fc6f40327c;hb=89a458bd80c63f5645b3f5fe26d2dda36e9489df;hp=a593f2e2976a8b25856dfde9be79a33e76b55240;hpb=bdcfac304c0c7e81083cdd601f7e5bbbe0eb255a;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG4/AliGammaMCReader.cxx b/PWG4/AliGammaMCReader.cxx index a593f2e2976..d606dacfb8d 100644 --- a/PWG4/AliGammaMCReader.cxx +++ b/PWG4/AliGammaMCReader.cxx @@ -9,7 +9,7 @@ * 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 * + * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id$ */ @@ -17,8 +17,11 @@ /* History of cvs commits: * * $Log$ - * Revision 1.1.2.1 2007/07/26 10:32:09 schutz - * new analysis classes in the the new analysis framework + * Revision 1.5 2007/11/17 16:41:07 gustavo + * mc handler should not be set in the analysis task, but in the startup analysis macro (MG) + * + * Revision 1.4 2007/10/29 13:48:42 gustavo + * Corrected coding violations * * */ @@ -39,36 +42,33 @@ #include #include #include +#include +#include //---- ANALYSIS system ---- #include "AliGammaMCReader.h" #include "Riostream.h" #include "AliLog.h" +#include "AliStack.h" ClassImp(AliGammaMCReader) //____________________________________________________________________________ AliGammaMCReader::AliGammaMCReader() : - AliGammaReader(), - fEMCALIPDistance(0.), fPHOSIPDistance(0.), - fEMCALMinDistance(0.), fPHOSMinDistance(0.), - fDecayPi0(kFALSE) + AliGammaReader(), fDecayPi0(0), fCheckOverlapping(kFALSE), fNeutralParticlesArray(0x0) { //Ctor //Initialize parameters - fDataType = kMC; InitParameters(); - + fDataType = kMC; } //____________________________________________________________________________ AliGammaMCReader::AliGammaMCReader(const AliGammaMCReader & g) : - AliGammaReader(g), - fEMCALIPDistance(g.fEMCALIPDistance), fPHOSIPDistance(g.fPHOSIPDistance), - fEMCALMinDistance(g.fEMCALMinDistance), fPHOSMinDistance(g.fPHOSMinDistance), - fDecayPi0(g.fDecayPi0) + AliGammaReader(g), fDecayPi0(g.fDecayPi0), fCheckOverlapping(g.fCheckOverlapping), + fNeutralParticlesArray(g.fNeutralParticlesArray?new TArrayI(*g.fNeutralParticlesArray):0x0) { // cpy ctor } @@ -80,24 +80,123 @@ AliGammaMCReader & AliGammaMCReader::operator = (const AliGammaMCReader & source if(&source == this) return *this; - fEMCALIPDistance = source.fEMCALIPDistance; - fPHOSIPDistance = source.fPHOSIPDistance; - fEMCALMinDistance = source.fEMCALMinDistance; - fPHOSMinDistance = source.fPHOSMinDistance; - fDecayPi0 = source.fDecayPi0; - - + fDecayPi0 = source.fDecayPi0; + fCheckOverlapping = source.fCheckOverlapping; + delete fNeutralParticlesArray; + fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0; return *this; } +//_________________________________ +AliGammaMCReader::~AliGammaMCReader() { + //Dtor + + delete fNeutralParticlesArray ; + +} + +//___________________________________________________________________________ +void AliGammaMCReader::CaseDecayGamma(Int_t index, TParticle * particle, AliStack * sta, + TClonesArray * plEMCAL, Int_t &indexEMCAL, + TClonesArray * plPHOS, Int_t &indexPHOS){ + //In case pi0 are decayed by pythia, check if mother is pi0 and in such case look if + //there is overlapping. Send particle=pi0 if there is overlapping. + + TParticle * pmother =sta->Particle(particle->GetFirstMother()); + if(pmother->GetPdgCode() == 111 && pmother->GetNDaughters() == 2) {//Do not consider 3 particle decay case + Int_t idaug0 = pmother->GetDaughter(0); + Int_t idaug1 = pmother->GetDaughter(1); + TParticle * pdaug0 = sta -> Particle(idaug0); + TParticle * pdaug1 = sta -> Particle(idaug1); + + if((index == idaug0 && pdaug0->Pt() > fNeutralPtCut) || + (index == idaug1 && pdaug0->Pt() <= fNeutralPtCut)){//Check decay when first daughter arrives, do nothing with second. + + FillListWithDecayGammaOrPi0(pmother, pdaug0, pdaug1, plEMCAL, indexEMCAL, plPHOS, indexPHOS); + + } + }//mother is a pi0 with 2 daughters + else{ + + if(IsInPHOS(particle->Phi(),particle->Eta())) { + TParticle * pmother =sta->Particle(particle->GetFirstMother()); + SetPhotonStatus(particle,pmother); + new((*plPHOS)[indexPHOS++]) TParticle(*particle) ; + } + else if(IsInEMCAL(particle->Phi(),particle->Eta())){ + TParticle * pmother =sta->Particle(particle->GetFirstMother()); + SetPhotonStatus(particle,pmother); + new((*plEMCAL)[indexEMCAL++]) TParticle(*particle) ; + } + } +} + +//___________________________________________________________________________ + void AliGammaMCReader::CaseGeantDecay(TParticle * particle, AliStack * stack, + TClonesArray * plEMCAL, Int_t &indexEMCAL, + TClonesArray * plPHOS, Int_t &indexPHOS){ + //Find decay gamma from pi0, decayed by GEANT and put them in the list. + + Int_t ndaug = particle->GetNDaughters() ; + TParticle * d1 = new TParticle(); + TParticle * d2 = new TParticle(); + + if(ndaug > 0){//At least 1 daugther + d1 = stack->Particle(particle->GetDaughter(0)); + if (ndaug > 1 ) + d2 = stack->Particle(particle->GetDaughter(1)); + + FillListWithDecayGammaOrPi0(particle, d1, d2, plEMCAL, indexEMCAL, plPHOS, indexPHOS); + + }// ndaugh > 0 + } + +//___________________________________________________________________________ +void AliGammaMCReader::CasePi0Decay(TParticle * pPi0, + TClonesArray * plEMCAL, Int_t &indexEMCAL, + TClonesArray * plPHOS, Int_t &indexPHOS){ + + //Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma + + TLorentzVector lvPi0, lvGamma1, lvGamma2 ; + //Double_t angle = 0; + + lvPi0.SetPxPyPzE(pPi0->Px(),pPi0->Py(),pPi0->Pz(),pPi0->Energy()); + + //Decay + MakePi0Decay(lvPi0,lvGamma1,lvGamma2);//,angle); + + //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it + TParticle * pPhoton1 = new TParticle(22,1,0,0,0,0,lvGamma1.Px(),lvGamma1.Py(), + lvGamma1.Pz(),lvGamma1.E(),0,0,0,0); + TParticle * pPhoton2 = new TParticle(22,1,0,0,0,0,lvGamma2.Px(),lvGamma2.Py(), + lvGamma2.Pz(),lvGamma2.E(),0,0,0,0); + + FillListWithDecayGammaOrPi0(pPi0, pPhoton1, pPhoton2, plEMCAL, indexEMCAL, plPHOS, indexPHOS); + +} + +//_______________________________________________________________ +void AliGammaMCReader::InitParameters() +{ + + //Initialize the parameters of the analysis. + + fDecayPi0 = kGeantDecay; + fCheckOverlapping = kTRUE ; + + Int_t pdgarray[]={12,14,16};// skip neutrinos + fNeutralParticlesArray = new TArrayI(3, pdgarray); + +} //____________________________________________________________________________ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *, - TClonesArray * plCh, - TClonesArray * plEMCAL, + TClonesArray * plCh, + TClonesArray * plEMCAL, TClonesArray * plPHOS, - TClonesArray *plParton) + TClonesArray *plParton,TClonesArray *,TClonesArray *) { //Create list of particles from EMCAL, PHOS and CTS. AliStack * stack = (AliStack *) data ; @@ -110,34 +209,50 @@ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *, for (iParticle=0 ; iParticle < stack->GetNprimary() ; iParticle++) { TParticle * particle = stack->Particle(iParticle); - + //Keep partons if(particle->GetStatusCode() == 21 && iParticle>=2){//All partons, not nucleus new((*plParton)[indexParton++]) TParticle(*particle) ; } //Keep Stable particles - if((particle->GetStatusCode() == 0) && (particle->Pt() > 0)){ + if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){ charge = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode())->Charge(); //---------- Charged particles ---------------------- - if((charge != 0) && (particle->Pt() > fChargedPtCut)){ + if( fSwitchOnCTS && (charge != 0) && (particle->Pt() > fChargedPtCut)){ //Particles in CTS acceptance if(TMath::Abs(particle->Eta())Pt() > fNeutralPtCut && - TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos + else if((charge == 0) && particle->Pt() > fNeutralPtCut ){ //&& + //TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos if(particle->GetPdgCode()!=111){ - if(IsInPHOS(particle->Phi(),particle->Eta())) - new((*plPHOS)[indexPHOS++]) TParticle(*particle) ; - else if(IsInEMCAL(particle->Phi(),particle->Eta())) - new((*plEMCAL)[indexEMCAL++]) TParticle(*particle) ; + //In case that we force PYTHIA to decay pi0, and we want to check the overlapping of + // the decay gamma. + if(particle->GetPdgCode() == 22 && fDecayPi0==kDecayGamma){ + CaseDecayGamma(iParticle,particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS); //If pythia decays pi0 + } + else{ + //Take out some particles like neutrinos + if(!SkipNeutralParticles(particle->GetPdgCode()) && !particle->IsPrimary()){ // protection added (MG) + TParticle * pmother =stack->Particle(particle->GetFirstMother()); + if(IsInPHOS(particle->Phi(),particle->Eta())){ + if(particle->GetPdgCode()==22)SetPhotonStatus(particle,pmother); + new((*plPHOS)[indexPHOS++]) TParticle(*particle) ; + } + else if(IsInEMCAL(particle->Phi(),particle->Eta())){ + if(particle->GetPdgCode()==22) SetPhotonStatus(particle,pmother); + new((*plEMCAL)[indexEMCAL++]) TParticle(*particle) ; + } + }//skip neutral particles + }//Case kDecayGamma }//no pi0 else{ if(fDecayPi0 == kNoDecay){//keep the pi0 do not decay @@ -145,148 +260,112 @@ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *, new((*plPHOS)[indexPHOS++]) TParticle(*particle) ; else if(IsInEMCAL(particle->Phi(),particle->Eta())) new((*plEMCAL)[indexEMCAL++]) TParticle(*particle) ; - } + }//nodecay else if(fDecayPi0 == kDecay) - MakePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS); + CasePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS); else if(fDecayPi0 == kGeantDecay) - SetGeantDecay(particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS); + CaseGeantDecay(particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS); }//pi0 }//neutral particle }//stable particle + else if(particle->GetStatusCode() == 11 && (particle->GetPdgCode() == 111 || particle->GetPdgCode() == 221) && particle->Pt() > 5 ){//hardcoded pt threshold to avoid too many particles + if(IsInPHOS(particle->Phi(),particle->Eta())) + new((*plPHOS)[indexPHOS++]) TParticle(*particle) ; + else if(IsInEMCAL(particle->Phi(),particle->Eta())) + new((*plEMCAL)[indexEMCAL++]) TParticle(*particle) ; + } }//particle loop } -//___________________________________________________________________________ -Bool_t AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta){ - //Check if particle is in EMCAL acceptance - if(phi<0) - phi+=TMath::TwoPi(); - if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] && - TMath::Abs(eta) fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] && - TMath::Abs(eta)GetNDaughters() ; - if(ndaug<=2 && ndaug >0){//At least 1 daugther - TParticle * d1 = stack->Particle(particle->GetDaughter(0)); - if(d1->GetPdgCode()==22){ - if(IsInEMCAL(d1->Phi(),d1->Eta())) - new((*plEMCAL)[indexEMCAL++]) TParticle(*d1) ; - else if(IsInPHOS(d1->Phi(),d1->Eta())) - new((*plPHOS)[indexPHOS++]) TParticle(*d1) ; + //Check if decay gamma overlapp in calorimeter, in such case keep the pi0, if not keep both photons. + Bool_t overlap = kFALSE ; + TLorentzVector lv1 , lv2 ; + pdaug0->Momentum(lv1); + pdaug1->Momentum(lv2); + + if(fCheckOverlapping){ + Double_t angle = lv1.Angle(lv2.Vect()); + + if(IsInEMCAL(pPi0->Phi(), pPi0->Eta())) + { + if (angle < fEMCALMinAngle){ + pPi0->SetStatusCode(1); + new((*plEMCAL)[indexEMCAL++]) TParticle(*pPi0) ; + overlap = kTRUE; + } + }//in EMCAL? + + else if(IsInPHOS(pPi0->Phi(), pPi0->Eta())){ + if (angle < fPHOSMinAngle){ + pPi0->SetStatusCode(1); + new((*plPHOS)[indexPHOS++]) TParticle(*pPi0) ; + overlap = kTRUE; + } + }//in PHOS? + + }//fCheckOverlapping + + //Fill with gammas if not overlapp + if(!overlap){ + if(pdaug0->GetPdgCode() == 22 || TMath::Abs(pdaug0->GetPdgCode() ) == 11 ){ + pdaug0->SetStatusCode(kPi0DecayPhoton); + if(IsInEMCAL(pdaug0->Phi(),pdaug0->Eta()) && pdaug0->Pt() > fNeutralPtCut) + new((*plEMCAL)[indexEMCAL++]) TParticle(*pdaug0) ; + else if(IsInPHOS(pdaug0->Phi(),pdaug0->Eta()) && pdaug0->Pt() > fNeutralPtCut) + new((*plPHOS)[indexPHOS++]) TParticle(*pdaug0) ; } - if(ndaug>1){//second daugther if present - TParticle * d2 = stack->Particle(particle->GetDaughter(1)); - if(IsInEMCAL(d2->Phi(),d2->Eta())) - new((*plEMCAL)[indexEMCAL++]) TParticle(*d2) ; - else if(IsInPHOS(d2->Phi(),d2->Eta())) - new((*plPHOS)[indexPHOS++]) TParticle(*d2) ; + if(pdaug1->GetPdgCode() == 22 || TMath::Abs(pdaug1->GetPdgCode() ) == 11 ){ + pdaug1->SetStatusCode(kPi0DecayPhoton); + if(IsInEMCAL(pdaug1->Phi(),pdaug1->Eta()) && pdaug1->Pt() > fNeutralPtCut) + new((*plEMCAL)[indexEMCAL++]) TParticle(*pdaug1) ; + else if(IsInPHOS(pdaug1->Phi(),pdaug1->Eta()) && pdaug1->Pt() > fNeutralPtCut) + new((*plPHOS)[indexPHOS++]) TParticle(*pdaug1) ; } - } -} + }// overlap? + + } + + //___________________________________________________________________________ -void AliGammaMCReader::MakePi0Decay(TParticle * particle, - TClonesArray * plEMCAL, Int_t &indexEMCAL, - TClonesArray * plPHOS, Int_t &indexPHOS){ - - //Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma - - TLorentzVector pPi0, pGamma1, pGamma2 ; - Double_t angle = 0, cellDistance = 0.; - Bool_t checkPhoton = kTRUE; - - pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->Energy()); - - //Decay - Pi0Decay(pPi0,pGamma1,pGamma2,angle); - - //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it - if(IsInPHOS(particle->Phi(), particle->Eta())){ - cellDistance = angle*fPHOSIPDistance; - if (cellDistance < fPHOSMinDistance){ - new((*plPHOS)[indexPHOS++]) TParticle(*particle) ; - checkPhoton = kFALSE; - } - } - else if(IsInEMCAL(particle->Phi(), particle->Eta())){ - cellDistance = angle*fEMCALIPDistance; - if (cellDistance < fEMCALMinDistance) { - new((*plEMCAL)[indexEMCAL++]) TParticle(*particle) ; - checkPhoton = kFALSE; - } - } - else checkPhoton = kTRUE ; - - if (checkPhoton) { - //Gamma Not overlapped - TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(), - pGamma1.Pz(),pGamma1.E(),0,0,0,0); - if(photon1->Pt() > fNeutralPtCut) { - - if(IsInPHOS(photon1->Phi(), photon1->Eta())) - new((*plPHOS)[indexPHOS++]) TParticle(*photon1) ; - - else if(IsInEMCAL(photon1->Phi(), photon1->Eta())) - new((*plEMCAL)[indexEMCAL++]) TParticle(*photon1) ; - - }// photon 1 of pi0 in acceptance - - - TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(), - pGamma2.Pz(),pGamma2.E(),0,0,0,0); - - if(photon2->Pt() > fNeutralPtCut) { - - if(IsInPHOS(photon2->Phi(), photon2->Eta())) - new((*plPHOS)[indexPHOS++]) TParticle(*photon2) ; - - else if(IsInEMCAL(photon2->Phi(), photon2->Eta())) - new((*plEMCAL)[indexEMCAL++]) TParticle(*photon2) ; - - }// photon 2 of pi0 in acceptance - }//Not overlapped gamma +Bool_t AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta) const { + //Check if particle is in EMCAL acceptance + + if(fSwitchOnEMCAL){ + if(phi<0) + phi+=TMath::TwoPi(); + if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] && + TMath::Abs(eta) fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] && + TMath::Abs(eta) 2 photons // p0 is pi0 4-momentum (inut) // p1 and p2 are photon 4-momenta (output) @@ -321,7 +400,7 @@ void AliGammaMCReader::Pi0Decay(TLorentzVector &p0, TLorentzVector &p1, p2.Boost(b); //cout<<"p2: "<GetStatusCode() != 21 ){ + if(pmother->GetStatusCode() ==11){ + if(pmother->GetPdgCode()==111) pphoton->SetStatusCode(kPi0DecayPhoton);//decay gamma from pi0 + if(pmother->GetPdgCode()==221) pphoton->SetStatusCode(kEtaDecayPhoton);//decay gamma from eta + else pphoton->SetStatusCode(kOtherDecayPhoton);// decay gamma from other pphotons + } + else pphoton->SetStatusCode(kUnknown); + } + else if(pmother->GetPdgCode()==22) pphoton->SetStatusCode(kPromptPhoton);//Prompt photon + else pphoton->SetStatusCode(kFragmentPhoton); //Fragmentation photon + +} + +//________________________________________________________________ +Bool_t AliGammaMCReader::SkipNeutralParticles(Int_t pdg) const { + //Check if pdg is equal to one of the neutral particles list + //These particles will be skipped from analysis. + + for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++) + if(TMath::Abs(pdg) == fNeutralParticlesArray->At(i)) return kTRUE ; + + return kFALSE; }