#include "TH2F.h"
#include "TBits.h"
#include "TList.h"
+#include "TArrayF.h"
+#include "TDecayChannel.h"
ClassImp(AliCFParticleGenCuts)
fDecayLengthMax(1.e+09),
fDecayRxyMin(-1),
fDecayRxyMax(1.e+09),
+ fDecayChannel(0x0),
fhCutStatistics(0x0),
fhCutCorrelation(0x0),
+ fCutValues(new TArrayF(kNCuts)),
fBitmap(new TBits(0))
{
//
fDecayLengthMax(1.e+09),
fDecayRxyMin(-1.),
fDecayRxyMax(1.e+09),
+ fDecayChannel(0x0),
fhCutStatistics(0x0),
fhCutCorrelation(0x0),
+ fCutValues(new TArrayF(kNCuts)),
fBitmap(new TBits(0))
{
//
fDecayLengthMax(c.fDecayLengthMin),
fDecayRxyMin(c.fDecayLengthMin),
fDecayRxyMax(c.fDecayLengthMin),
+ fDecayChannel(c.fDecayChannel),
fhCutStatistics(new TH1F(*c.fhCutStatistics)),
fhCutCorrelation(new TH2F(*c.fhCutCorrelation)),
+ fCutValues(new TArrayF(*c.fCutValues)),
fBitmap(new TBits(*c.fBitmap))
{
//
//copy ctor
//
- if (c.fhCutStatistics) fhCutStatistics = new TH1F(*c.fhCutStatistics) ;
- if (c.fhCutCorrelation) fhCutCorrelation = new TH2F(*c.fhCutCorrelation) ;
-
for (int i=0; i<kNCuts; i++)
for (int j=0; j<kNStepQA; j++)
fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
fDecayLengthMax=c.fDecayLengthMax;
fDecayRxyMin=c.fDecayRxyMin;
fDecayRxyMax=c.fDecayRxyMax;
+ fDecayChannel=c.fDecayChannel;
+ fCutValues=new TArrayF(*c.fCutValues);
fBitmap=new TBits(*c.fBitmap);
if (fhCutStatistics) fhCutStatistics =new TH1F(*c.fhCutStatistics) ;
if (fIsQAOn) FillHistograms(obj,1);
return kTRUE;
}
-
+
//__________________________________________________________________________________
void AliCFParticleGenCuts::SelectionBitMap(TObject* obj)
{
// test if the track passes the single cuts
// and store the information in a bitmap
//
-
- for (UInt_t i=0; i<kNCuts; i++) fBitmap->SetBitNumber(i,kFALSE);
+
+ for (UInt_t i=0; i<kNCuts; i++) {
+ fBitmap->SetBitNumber(i,kFALSE);
+ fCutValues->SetAt((Double32_t)0,i) ;
+ }
if (!obj) return ;
TString className(obj->ClassName());
AliError("argument must point to an AliMCParticle !");
return ;
}
-
- AliMCParticle* mcPart = (AliMCParticle*) obj ;
+
+ AliMCParticle* mcPart = dynamic_cast<AliMCParticle*>(obj) ;
TParticle* part = mcPart->Particle();
AliStack* stack = fMCInfo->Stack();
+
+ // fill the cut array
Double32_t partVx=(Double32_t)part->Vx();
Double32_t partVy=(Double32_t)part->Vy();
Double32_t partVz=(Double32_t)part->Vz();
decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
}
-
- Int_t iCutBit = 0;
-
+ fCutValues->SetAt(partVx ,kCutProdVtxXMin) ;
+ fCutValues->SetAt(partVy ,kCutProdVtxYMin) ;
+ fCutValues->SetAt(partVz ,kCutProdVtxZMin) ;
+ fCutValues->SetAt(partVx ,kCutProdVtxXMax) ;
+ fCutValues->SetAt(partVy ,kCutProdVtxYMax) ;
+ fCutValues->SetAt(partVz ,kCutProdVtxZMax) ;
+ fCutValues->SetAt(decayVx ,kCutDecVtxXMin) ;
+ fCutValues->SetAt(decayVy ,kCutDecVtxYMin) ;
+ fCutValues->SetAt(decayVz ,kCutDecVtxZMin) ;
+ fCutValues->SetAt(decayVx ,kCutDecVtxXMax) ;
+ fCutValues->SetAt(decayVy ,kCutDecVtxYMax) ;
+ fCutValues->SetAt(decayVz ,kCutDecVtxZMax) ;
+ fCutValues->SetAt(decayL ,kCutDecLgthMin) ;
+ fCutValues->SetAt(decayL ,kCutDecLgthMax) ;
+ fCutValues->SetAt(decayRxy,kCutDecRxyMin) ;
+ fCutValues->SetAt(decayRxy,kCutDecRxyMax) ;
+
// cut on charge
if ( fRequireIsCharged || fRequireIsNeutral ) {
- if(fRequireIsCharged && IsCharged(mcPart)) fBitmap->SetBitNumber(iCutBit,kTRUE);
- if(fRequireIsNeutral && !IsCharged(mcPart)) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
+ if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
}
- else fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
+ else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
// cut on primary/secondary
if ( fRequireIsPrimary || fRequireIsSecondary) {
- if (fRequireIsPrimary && IsPrimary(mcPart,stack)) fBitmap->SetBitNumber(iCutBit,kTRUE);
- if (fRequireIsSecondary && !IsPrimary(mcPart,stack)) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if (fRequireIsPrimary && IsPrimary(mcPart,stack)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
+ if (fRequireIsSecondary && !IsPrimary(mcPart,stack)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
}
- else fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
-
+ else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
+
// cut on PDG code
- if ( fRequirePdgCode ){
- if (IsA(mcPart,fPdgCode,kTRUE)) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fRequirePdgCode ) {
+ if (IsA(mcPart,fPdgCode,kFALSE)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
+ }
+ else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
+
+ // cut on decay channel
+ if ( fDecayChannel ) {
+ Bool_t goodDecay = kTRUE ;
+ Short_t nDaughters = mcPart->Particle()->GetNDaughters() ;
+ if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
+ if (goodDecay) {
+ // now check if decay channel is respected
+ // first try
+ for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
+ TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
+ if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+ }
+ if (!goodDecay) {
+ //second try inverting the order of the daughters
+ goodDecay = kTRUE ;
+ for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
+ TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
+ if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+ }
+ }
+ }
+ fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
}
- else fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
-
- // production vertex cuts
- if ( partVx > fProdVtxXMin ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( partVx < fProdVtxXMax ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( partVy > fProdVtxYMin ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( partVy < fProdVtxYMax ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( partVz > fProdVtxZMin ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( partVz < fProdVtxZMax ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
+ else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
- // decay vertex cuts
- if ( decayVx > fDecayVtxXMin ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( decayVx < fDecayVtxXMax ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( decayVy > fDecayVtxYMin ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( decayVy < fDecayVtxYMax ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( decayVz > fDecayVtxZMin ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( decayVz < fDecayVtxZMax ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
-
- // decay length cuts
- if ( decayL > fDecayLengthMin ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( decayL < fDecayLengthMax ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- // transverse decay length cuts
- if ( decayRxy > fDecayRxyMin ) fBitmap->SetBitNumber(iCutBit,kTRUE);
- iCutBit++;
- if ( decayRxy < fDecayRxyMax ) fBitmap->SetBitNumber(iCutBit,kTRUE);
-
+ // now array of cut is build, fill the bitmap consequently
+ Int_t iCutBit = -1;
+ 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);
+ 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);
+ if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
+ if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
}
+
//__________________________________________________________________________________
-void AliCFParticleGenCuts::FillHistograms(TObject* obj, Bool_t afterCuts)
+void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
{
//
// fill the QA histograms
//
for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++)
- fhQA[iCutNumber][afterCuts]->Fill(fBitmap->TestBitNumber(iCutNumber));
+ fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
// fill cut statistics and cut correlation histograms with information from the bitmap
if (afterCuts) return;
// Number of single cuts in this class
UInt_t ncuts = fBitmap->GetNbits();
for(UInt_t bit=0; bit<ncuts;bit++) {
- if (!fBitmap->TestBitNumber(bit)) {
+ if (!fBitmap->TestBitNumber(bit)) {
fhCutStatistics->Fill(bit+1);
for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
if (!fBitmap->TestBitNumber(bit2))
fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin") ; k++;
fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax") ; k++;
+ fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
for (int i=0; i<kNStepQA; i++) {
if (i==0) sprintf(str," ");
else sprintf(str,"_cut");
- fhQA[kCutCharge] [i] = new TH1F(Form("%s_charge%s" ,GetName(),str),"",2,-0.5,1.5);
- fhQA[kCutPrimSec] [i] = new TH1F(Form("%s_primSec%s" ,GetName(),str),"",2,-0.5,1.5);
- fhQA[kCutPDGCode] [i] = new TH1F(Form("%s_pdgCode%s" ,GetName(),str),"",10000,-5000,5000);
- 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[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);
- fhQA[kCutDecVtxYMax] [i] = new TH1F(Form("%s_decVtxYMax%s" ,GetName(),str),"",100,0,10);
- fhQA[kCutDecVtxZMin] [i] = new TH1F(Form("%s_decVtxZMin%s" ,GetName(),str),"",100,0,10);
- fhQA[kCutDecVtxZMax] [i] = new TH1F(Form("%s_decVtxZMax%s" ,GetName(),str),"",100,0,10);
- fhQA[kCutDecLgthMin] [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
- fhQA[kCutDecLgthMax] [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
- 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[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[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);
+ fhQA[kCutDecVtxYMax] [i] = new TH1F(Form("%s_decVtxYMax%s" ,GetName(),str),"",100,0,10);
+ fhQA[kCutDecVtxZMin] [i] = new TH1F(Form("%s_decVtxZMin%s" ,GetName(),str),"",100,0,10);
+ fhQA[kCutDecVtxZMax] [i] = new TH1F(Form("%s_decVtxZMax%s" ,GetName(),str),"",100,0,10);
+ fhQA[kCutDecLgthMin] [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
+ fhQA[kCutDecLgthMax] [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
+ 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);
}
for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
}
//
TParticle* part = mcPart->Particle();
Int_t pdgCode = part->GetPdgCode();
- if (abs) pdgCode = TMath::Abs(pdgCode);
+ if (abs) {
+ pdgCode = TMath::Abs(pdgCode);
+ pdg = TMath::Abs(pdg);
+ }
if (pdgCode != pdg ) return kFALSE;
return kTRUE;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff