#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliPIDResponse.h"
+#include "TRandom.h"
+
+using std::cout;
+using std::endl;
ClassImp(AliRDHFCuts)
fMinSPDMultiplicity(0),
fTriggerMask(AliVEvent::kAnyINT),
fUseOnlyOneTrigger(kFALSE),
-fTriggerClass("CINT1"),
fTrackCuts(0),
fnPtBins(1),
fnPtBinLimits(1),
fIsCandTrackSPDFirst(kFALSE),
fMaxPtCandTrackSPDFirst(0.),
fApplySPDDeadPbPb2011(kFALSE),
-fRemoveTrackletOutliers(kFALSE)
+fRemoveTrackletOutliers(kFALSE),
+fCutOnzVertexSPD(0),
+fKinkReject(kFALSE),
+ fUseTrackSelectionWithFilterBits(kTRUE),
+ fHistCentrDistr(0x0)
{
//
// Default Constructor
//
+ fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
}
//--------------------------------------------------------------------------
AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
fMinSPDMultiplicity(source.fMinSPDMultiplicity),
fTriggerMask(source.fTriggerMask),
fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
- fTriggerClass(source.fTriggerClass),
+ fTriggerClass(),
fTrackCuts(0),
fnPtBins(source.fnPtBins),
fnPtBinLimits(source.fnPtBinLimits),
fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst),
fApplySPDDeadPbPb2011(source.fApplySPDDeadPbPb2011),
- fRemoveTrackletOutliers(source.fRemoveTrackletOutliers)
+ fRemoveTrackletOutliers(source.fRemoveTrackletOutliers),
+ fCutOnzVertexSPD(source.fCutOnzVertexSPD),
+ fKinkReject(source.fKinkReject),
+ fUseTrackSelectionWithFilterBits(source.fUseTrackSelectionWithFilterBits),
+ fHistCentrDistr(0x0)
{
//
// Copy constructor
//
cout<<"Copy constructor"<<endl;
+ fTriggerClass[0] = source.fTriggerClass[0];
+ fTriggerClass[1] = source.fTriggerClass[1];
if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
if(source.fPidHF) SetPidHF(source.fPidHF);
+ if(source.fHistCentrDistr) fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
PrintAll();
}
fMinSPDMultiplicity=source.fMinSPDMultiplicity;
fTriggerMask=source.fTriggerMask;
fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
- fTriggerClass=source.fTriggerClass;
+ fTriggerClass[0]=source.fTriggerClass[0];
+ fTriggerClass[1]=source.fTriggerClass[1];
fnPtBins=source.fnPtBins;
fnPtBinLimits=source.fnPtBinLimits;
fnVars=source.fnVars;
fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
fApplySPDDeadPbPb2011=source.fApplySPDDeadPbPb2011;
fRemoveTrackletOutliers=source.fRemoveTrackletOutliers;
+ fCutOnzVertexSPD=source.fCutOnzVertexSPD;
+ fKinkReject=source.fKinkReject;
+ fUseTrackSelectionWithFilterBits=source.fUseTrackSelectionWithFilterBits;
+ if(fHistCentrDistr) delete fHistCentrDistr;
+ if(source.fHistCentrDistr)fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
delete fPidHF;
fPidHF=0;
}
+ if(fHistCentrDistr)delete fHistCentrDistr;
}
//---------------------------------------------------------------------------
Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
}else{
if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
return 2;
- }
+ }
+ // selection to flatten the centrality distribution
+ if(fHistCentrDistr){
+ if(!IsEventSelectedForCentrFlattening(centvalue))return 4;
+ }
}
}
return 0;
}
+
+
+//-------------------------------------------------
+void AliRDHFCuts::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
+ // set the histo for centrality flattening
+ // the centrality is flatten in the range minCentr,maxCentr
+ // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
+ // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
+ // negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat)
+ // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom
+
+ if(maxCentr<minCentr){
+ AliWarning("AliRDHFCuts::Wrong centralities values while setting the histogram for centrality flattening");
+ }
+
+ if(fHistCentrDistr)delete fHistCentrDistr;
+ fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
+ fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
+ Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
+ Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
+ fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
+ Double_t ref=0.,bincont=0.,binrefwidth=1.;
+ Int_t binref=0;
+ if(TMath::Abs(centrRef)<0.0001){
+ binref=fHistCentrDistr->GetMinimumBin();
+ binrefwidth=fHistCentrDistr->GetBinWidth(binref);
+ ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
+ }
+ else if(centrRef>0.){
+ binref=h->FindBin(centrRef);
+ if(binref<1||binref>h->GetNbinsX()){
+ AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
+ }
+ binrefwidth=fHistCentrDistr->GetBinWidth(binref);
+ ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
+ }
+ else{
+ if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
+ binref=fHistCentrDistr->GetMaximumBin();
+ binrefwidth=fHistCentrDistr->GetBinWidth(binref);
+ ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
+ }
+
+ for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
+ if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
+ bincont=h->GetBinContent(j);
+ fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
+ fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
+ }
+ else{
+ h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
+ }
+ }
+
+ fHistCentrDistr->SetBinContent(0,switchTRand);
+ return;
+
+}
+
+//-------------------------------------------------
+Bool_t AliRDHFCuts::IsEventSelectedForCentrFlattening(Float_t centvalue){
+ //
+ // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
+ // Can be faster if it was required that fHistCentrDistr covers
+ // exactly the desired centrality range (e.g. part of the lines below should be done during the
+ // setting of the histo) and TH1::SetMinimum called
+ //
+
+ if(!fHistCentrDistr) return kTRUE;
+ // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
+ // if(maxbin>fHistCentrDistr->GetNbinsX()){
+ // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
+ // }
+
+ Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
+ Double_t bincont=fHistCentrDistr->GetBinContent(bin);
+ Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
+
+ if(fHistCentrDistr->GetBinContent(0)<-0.9999){
+ if(gRandom->Uniform(1.)<bincont)return kTRUE;
+ return kFALSE;
+ }
+
+ if(centDigits*100.<bincont)return kTRUE;
+ return kFALSE;
+
+}
+//---------------------------------------------------------------------------
+void AliRDHFCuts::SetupPID(AliVEvent *event) {
+ // Set the PID response object in the AliAODPidHF
+ // in case of old PID sets the TPC dE/dx BB parameterization
+
+ if(fPidHF){
+ if(fPidHF->GetPidResponse()==0x0){
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
+ AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
+ fPidHF->SetPidResponse(pidResp);
+ }
+ if(fPidHF->GetUseCombined()) fPidHF->SetUpCombinedPID();
+ if(fPidHF->GetOldPid()) {
+
+ Bool_t isMC=kFALSE;
+ TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+ if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
+
+ // pp, from LHC10d onwards
+ if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
+ event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
+ // pp, 2011 low energy run
+ if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
+ fPidHF->SetppLowEn2011(kTRUE);
+ fPidHF->SetOnePad(kFALSE);
+ }
+ // PbPb LHC10h
+ if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
+ // MC
+ if(isMC) fPidHF->SetMC(kTRUE);
+ if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
+ fPidHF->SetMClowenpp2011(kTRUE);
+ fPidHF->SetBetheBloch();
+ }else{
+ // check that AliPIDResponse object was properly set in case of using OADB
+ if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
+ }
+ }
+}
//---------------------------------------------------------------------------
Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
//
//
//if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
+
+
fWhyRejection=0;
fEvRejectionBits=0;
Bool_t accept=kTRUE;
TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
- // settings for the TPC dE/dx BB parameterization
- if(fPidHF && fPidHF->GetOldPid()) {
- // pp, from LHC10d onwards
- if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
- event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
- // pp, 2011 low energy run
- if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
- fPidHF->SetppLowEn2011(kTRUE);
- fPidHF->SetOnePad(kFALSE);
- }
- // PbPb LHC10h
- if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
- // MC
- if(isMC) fPidHF->SetMC(kTRUE);
- if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
- fPidHF->SetMClowenpp2011(kTRUE);
- fPidHF->SetBetheBloch();
- }
- else if(fPidHF && !fPidHF->GetOldPid()) {
- if(fPidHF->GetPidResponse()==0x0){
- AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
- AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
- AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
- fPidHF->SetPidResponse(pidResp);
- }
- }
+ SetupPID(event);
// trigger class
TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
// don't do for MC and for PbPb 2010 data
if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
- if(!firedTriggerClasses.Contains(fTriggerClass.Data())) {
+ if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) &&
+ (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
fWhyRejection=5;
fEvRejectionBits+=1<<kNotSelTrigger;
accept=kFALSE;
}
}
+ // TEMPORARY FIX FOR GetEvent
+ Int_t nTracks=((AliAODEvent*)event)->GetNTracks();
+ for(Int_t itr=0; itr<nTracks; itr++){
+ AliAODTrack* tr=(AliAODTrack*)((AliAODEvent*)event)->GetTrack(itr);
+ tr->SetAODEvent((AliAODEvent*)event);
+ }
+
// TEMPORARY FIX FOR REFERENCES
// Fix references to daughter tracks
// if(fFixRefs) {
}
}
+ if(fCutOnzVertexSPD>0){
+ const AliVVertex *vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
+ if(!vSPD || (vSPD && vSPD->GetNContributors()<fMinVtxContr)){
+ accept=kFALSE;
+ fEvRejectionBits+=1<<kBadSPDVertex;
+ }else{
+ if(fCutOnzVertexSPD==1 && TMath::Abs(vSPD->GetZ())>12.) {
+ fEvRejectionBits+=1<<kZVtxSPDOutFid;
+ if(accept) fWhyRejection=6;
+ accept=kFALSE;
+ }
+ if(fCutOnzVertexSPD==2 && vertex){
+ if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) {
+ fEvRejectionBits+=1<<kZVtxSPDOutFid;
+ if(accept) fWhyRejection=6;
+ accept=kFALSE;
+ }
+ }
+ }
+ }
// pile-up rejection
if(fOptPileup==kRejectPileupEvent){
Int_t rejection=IsEventSelectedInCentrality(event);
if(rejection>1){
if(accept) fWhyRejection=rejection;
- fEvRejectionBits+=1<<kOutsideCentrality;
+ if(fWhyRejection==4)fEvRejectionBits+=1<<kCentralityFlattening;
+ else fEvRejectionBits+=1<<kOutsideCentrality;
accept=kFALSE;
}
+
}
// PbPb2011 outliers in tracklets vs. VZERO
- if(!isMC && event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
+ if(event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
if(fRemoveTrackletOutliers){
Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
- if(v0cent<cutval){
+ if(ntracklets<1000. && v0cent<cutval){
if(accept) fWhyRejection=2;
fEvRejectionBits+=1<<kOutsideCentrality;
accept=kFALSE;
esdTrack.RelateToVertex(primary,0.,3.);
if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
+
+
+ if(fKinkReject){
+ AliAODVertex *maybeKink=track->GetProdVertex();
+ if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE;
+ }
if(fOptPileup==kRejectTracksFromPileupVertex){
// to be implemented
esdTrack.GetXYZAt(7.6,0.,xyz2);
Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
if(phi1<0) phi1+=2*TMath::Pi();
- Int_t lad1=phi1/(2.*TMath::Pi()/20.);
+ Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
if(phi2<0) phi2+=2*TMath::Pi();
- Int_t lad2=phi2/(2.*TMath::Pi()/40.);
- Int_t mod1=(xyz1[2]+14)/7.;
- Int_t mod2=(xyz2[2]+14)/7.;
+ Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
+ Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
+ Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
Bool_t lay1ok=kFALSE;
if(mod1>=0 && mod1<4 && lad1<20){
lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
printf("Minimum vtx contr %d\n",fMinVtxContr);
printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
printf("Min SPD mult %d\n",fMinSPDMultiplicity);
- printf("Use PID %d\n",(Int_t)fUsePID);
+ printf("Use PID %d OldPid=%d\n",(Int_t)fUsePID,fPidHF ? fPidHF->GetOldPid() : -1);
printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
//--------------------------------------------------------------------------
void AliRDHFCuts::PrintTrigger() const{
+ // print the trigger selection
+
+ printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
cout<<" Trigger selection pattern: ";
if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
if(cent<0){
Int_t quality = centrality->GetQuality();
- if(quality<=1){
+ if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
}else{
Int_t runnum=aodEvent->GetRunNumber();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff