X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=CORRFW%2FAliCFParticleGenCuts.cxx;h=8b6336bc2121b6cf3ad84bb2bc83fea0226e37b3;hb=fc0882f371b96f1a712eaa5f3475a2ae575fb254;hp=8c5b69933fb7c06ded76dc0d32333740baa70910;hpb=b95f6a3670bc7dabdaee85ccf211a7e89811cc72;p=u%2Fmrichter%2FAliRoot.git diff --git a/CORRFW/AliCFParticleGenCuts.cxx b/CORRFW/AliCFParticleGenCuts.cxx index 8c5b69933fb..8b6336bc212 100644 --- a/CORRFW/AliCFParticleGenCuts.cxx +++ b/CORRFW/AliCFParticleGenCuts.cxx @@ -50,6 +50,8 @@ AliCFParticleGenCuts::AliCFParticleGenCuts() : fRequireIsPrimary(0), fRequireIsSecondary(0), fRequirePdgCode(0), + fRequireAbsolutePdg(0), + fProdVtxRange2D(0), fPdgCode(0), fProdVtxXMin (-1.e+09), fProdVtxYMin (-1.e+09), @@ -79,6 +81,9 @@ AliCFParticleGenCuts::AliCFParticleGenCuts() : for (int i=0; iClone(); + + for (int j=0; jClone(); } //______________________________ @@ -178,6 +193,8 @@ AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts fRequireIsPrimary=c.fRequireIsPrimary; fRequireIsSecondary=c.fRequireIsSecondary; fRequirePdgCode=c.fRequirePdgCode; + fRequireAbsolutePdg=c.fRequireAbsolutePdg; + fProdVtxRange2D=c.fProdVtxRange2D; fPdgCode=c.fPdgCode; fProdVtxXMin=c.fProdVtxXMin; fProdVtxYMin=c.fProdVtxYMin; @@ -205,6 +222,9 @@ AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts for (int i=0; iClone(); + + for (int j=0; jClone(); } return *this ; } @@ -248,6 +268,14 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart) Double32_t partVy=(Double32_t)mcPart->Yv(); Double32_t partVz=(Double32_t)mcPart->Zv(); + // calculate the production vertex ellipse + Double32_t prodVtxXYmin = 0.; + if (fProdVtxXMin!=0 && fProdVtxYMin!=0) + prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin); + Double32_t prodVtxXYmax = 0.; + if(fProdVtxXMax!=0 && fProdVtxYMax!=0) + prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax); + Double32_t decayVx=0.; Double32_t decayVy=0.; Double32_t decayVz=0.; @@ -301,7 +329,7 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart) // cut on PDG code if ( fRequirePdgCode ) { - if (IsA(mcPart,fPdgCode,kFALSE)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); + if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); } else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); @@ -310,6 +338,7 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart) Bool_t goodDecay = kTRUE ; Short_t nDaughters = mcPart->Particle()->GetNDaughters() ; if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ; + //now number of daughters is OK if (goodDecay) { // now check if decay channel is respected // first try @@ -325,9 +354,26 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart) if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} } } - } + if (!goodDecay && fRequireAbsolutePdg) { + //now tries inverting the sign of the daughters in case the anti-particle is also looked at + // third try + goodDecay = kTRUE ; + for (Int_t iDaughter = 0; iDaughterParticle(mcPart->Particle()->GetDaughter(iDaughter)) ; + if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} + } + if (!goodDecay) { + //fourth try inverting the order of the daughters + goodDecay = kTRUE ; + for (Int_t iDaughter = 0; iDaughterParticle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ; + if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} + } + } + } //end check anti-particle + } //end # daughters OK fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ; - } + } //end require decay channel else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel); @@ -336,10 +382,23 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart) if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); - if ( fCutValues->At(++iCutBit) > fProdVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE); - if ( fCutValues->At(++iCutBit) < fProdVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE); - if ( fCutValues->At(++iCutBit) > fProdVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE); - if ( fCutValues->At(++iCutBit) < fProdVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE); + + ++iCutBit; + if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin) + || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE); + + ++iCutBit; + if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax) + || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE); + + ++iCutBit; + if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin) + || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE); + + ++iCutBit; + if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax) + || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE); + if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE); if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE); if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE); @@ -373,6 +432,14 @@ void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart) Double32_t partVy=(Double32_t)mcPart->Yv(); Double32_t partVz=(Double32_t)mcPart->Zv(); + // calculate the production vertex ellipse + Double32_t prodVtxXYmin = 0.; + if (fProdVtxXMin!=0 && fProdVtxYMin!=0) + prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin); + Double32_t prodVtxXYmax = 0.; + if(fProdVtxXMax!=0 && fProdVtxYMax!=0) + prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax); + Double32_t decayVx=0.; Double32_t decayVy=0.; Double32_t decayVz=0.; @@ -429,7 +496,7 @@ void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart) // cut on PDG code if ( fRequirePdgCode ) { - if (IsA(mcPart,fPdgCode,kFALSE)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); + if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); } else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); @@ -466,10 +533,23 @@ void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart) if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); - if ( fCutValues->At(++iCutBit) > fProdVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE); - if ( fCutValues->At(++iCutBit) < fProdVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE); - if ( fCutValues->At(++iCutBit) > fProdVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE); - if ( fCutValues->At(++iCutBit) < fProdVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE); + + ++iCutBit; + if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin) + || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE); + + ++iCutBit; + if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax) + || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE); + + ++iCutBit; + if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin) + || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE); + + ++iCutBit; + if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax) + || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE); + if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE); if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE); if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE); @@ -496,6 +576,8 @@ void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts) for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++) fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber)); + fhProdVtxXY[afterCuts]->Fill(fCutValues->At(4),fCutValues->At(5)); + // fill cut statistics and cut correlation histograms with information from the bitmap if (afterCuts) return; @@ -524,6 +606,7 @@ void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) { qaList->Add(fhCutCorrelation); for (Int_t j=0; jAdd(fhProdVtxXY[j]); for(Int_t i=0; iAdd(fhQA[i][j]); } @@ -576,12 +659,12 @@ void AliCFParticleGenCuts::DefineHistograms() { fhQA[kCutCharge] [i] = new TH1F(Form("%s_charge%s" ,GetName(),str),"",2,0,2); fhQA[kCutPrimSec] [i] = new TH1F(Form("%s_primSec%s" ,GetName(),str),"",2,0,2); fhQA[kCutPDGCode] [i] = new TH1F(Form("%s_pdgCode%s" ,GetName(),str),"",2,0,2); - fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,0,10); - fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,0,10); - fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,0,10); - fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,0,10); - fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,0,10); - fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,0,10); + fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,-10,10); + fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,-10,10); + fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,-10,10); + fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,-10,10); + fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,-10,10); + fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,-10,10); fhQA[kCutDecVtxXMin] [i] = new TH1F(Form("%s_decVtxXMin%s" ,GetName(),str),"",100,0,10); fhQA[kCutDecVtxXMax] [i] = new TH1F(Form("%s_decVtxXMax%s" ,GetName(),str),"",100,0,10); fhQA[kCutDecVtxYMin] [i] = new TH1F(Form("%s_decVtxYMin%s" ,GetName(),str),"",100,0,10); @@ -593,6 +676,11 @@ void AliCFParticleGenCuts::DefineHistograms() { fhQA[kCutDecRxyMin] [i] = new TH1F(Form("%s_decRxyMin%s" ,GetName(),str),"",100,0,10); fhQA[kCutDecRxyMax] [i] = new TH1F(Form("%s_decRxyMax%s" ,GetName(),str),"",100,0,10); fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2); + fhProdVtxXY [i] = new TH2F(Form("%s_prodVtxXY%s" ,GetName(),str),"",100,0,10,100,0,10); + fhProdVtxXY [i] ->GetXaxis()->SetTitle("x_{production vertex}"); + fhProdVtxXY [i] ->GetYaxis()->SetTitle("y_{production vertex}"); + fhQA[kCutProdVtxXMax] [i] ->GetXaxis()->SetTitle("x_{production vertex}"); + fhQA[kCutProdVtxYMax] [i] ->GetXaxis()->SetTitle("y_{production vertex}"); } for(Int_t i=0; iSetLineColor(color); } @@ -644,7 +732,7 @@ Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) { Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) { // //Check on the pdg code of the MC particle. if abs=kTRUE then check on the - //absolute value. By default is set to kFALSE. + //absolute value. // TParticle* part = mcPart->Particle(); Int_t pdgCode = part->GetPdgCode(); @@ -660,7 +748,7 @@ Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) { Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) { // //Check on the pdg code of the MC particle. if abs=kTRUE then check on the - //absolute value. By default is set to kFALSE. + //absolute value. // Int_t pdgCode = mcPart->GetPdgCode(); @@ -672,7 +760,7 @@ Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs return kTRUE; } //______________________________ -void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) { +void AliCFParticleGenCuts::SetMCEventInfo(const TObject* mcEvent) { // // Sets pointer to MC event information (AliMCEvent) // @@ -687,7 +775,6 @@ void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) { AliError("argument must point to an AliMCEvent or an AliAODEvent !"); return ; } - - if (fIsAODMC) fMCInfo = dynamic_cast(mcEvent) ; - else fMCInfo = dynamic_cast (mcEvent) ; + + fMCInfo = (AliVEvent*)mcEvent ; }