#include <TList.h>
#include <TTree.h>
#include <TStopwatch.h>
+#include "TRandom.h"
#include "AliLog.h"
#include "AliEventplane.h"
fUseCentrality(kFALSE),
fCentralityType("QUALITY"),
fUseAOD049CentralityPatch(kFALSE),
+ fUseCentralityPatchPbPb2011(0),
fContinuousMix(kTRUE),
fNMix(0),
fMaxDiffMult(10),
fBigOutput(kFALSE),
fMixPrintRefresh(-1),
fMaxNDaughters(-1),
- fCheckP(kFALSE)
+ fCheckP(kFALSE),
+ fCheckFeedDown(kFALSE),
+ fOriginDselection(kFALSE),
+ fKeepDfromB(kFALSE),
+ fKeepDfromBOnly(kFALSE),
+ fRejectIfNoQuark(kFALSE)
{
//
// Dummy constructor ALWAYS needed for I/O.
fUseCentrality(kFALSE),
fCentralityType("QUALITY"),
fUseAOD049CentralityPatch(kFALSE),
+ fUseCentralityPatchPbPb2011(0),
fContinuousMix(kTRUE),
fNMix(0),
fMaxDiffMult(10),
fBigOutput(kFALSE),
fMixPrintRefresh(-1),
fMaxNDaughters(-1),
- fCheckP(kFALSE)
+ fCheckP(kFALSE),
+ fCheckFeedDown(kFALSE),
+ fOriginDselection(kFALSE),
+ fKeepDfromB(kFALSE),
+ fKeepDfromBOnly(kFALSE),
+ fRejectIfNoQuark(kFALSE)
{
//
// Default constructor.
fUseCentrality(copy.fUseCentrality),
fCentralityType(copy.fCentralityType),
fUseAOD049CentralityPatch(copy.fUseAOD049CentralityPatch),
+ fUseCentralityPatchPbPb2011(copy.fUseCentralityPatchPbPb2011),
fContinuousMix(copy.fContinuousMix),
fNMix(copy.fNMix),
fMaxDiffMult(copy.fMaxDiffMult),
fBigOutput(copy.fBigOutput),
fMixPrintRefresh(copy.fMixPrintRefresh),
fMaxNDaughters(copy.fMaxNDaughters),
- fCheckP(copy.fCheckP)
+ fCheckP(copy.fCheckP),
+ fCheckFeedDown(copy.fCheckFeedDown),
+ fOriginDselection(copy.fOriginDselection),
+ fKeepDfromB(copy.fOriginDselection),
+ fKeepDfromBOnly(copy.fKeepDfromBOnly),
+ fRejectIfNoQuark(copy.fRejectIfNoQuark)
{
//
// Copy constructor.
fUseCentrality = copy.fUseCentrality;
fCentralityType = copy.fCentralityType;
fUseAOD049CentralityPatch = copy.fUseAOD049CentralityPatch;
+ fUseCentralityPatchPbPb2011 = copy.fUseCentralityPatchPbPb2011;
fContinuousMix = copy.fContinuousMix;
fNMix = copy.fNMix;
fMaxDiffMult = copy.fMaxDiffMult;
fMixPrintRefresh = copy.fMixPrintRefresh;
fMaxNDaughters = copy.fMaxNDaughters;
fCheckP = copy.fCheckP;
+ fCheckFeedDown = copy.fCheckFeedDown;
+ fOriginDselection = copy.fOriginDselection;
+ fKeepDfromB = copy.fOriginDselection;
+ fKeepDfromBOnly = copy.fKeepDfromBOnly;
+ fRejectIfNoQuark = copy.fRejectIfNoQuark;
return (*this);
}
//
if (fUseCentrality) {
-
- if ((!fUseMC) && (!isESD) && (fUseAOD049CentralityPatch)) {
- return ApplyCentralityPatchAOD049();
- } else {
- AliCentrality *centrality = fInputEvent->GetCentrality();
+ if ((!fUseMC) && (fUseCentralityPatchPbPb2011)) {
+ return ApplyCentralityPatchPbPb2011();//
+ }
+ if ((!fUseMC) && (!isESD) && (fUseAOD049CentralityPatch)) {
+ return ApplyCentralityPatchAOD049();
+ } else {
+ AliCentrality *centrality = fInputEvent->GetCentrality();
if (!centrality) {
- AliError("Cannot compute centrality!");
- return -1.0;
+ AliError("Cannot compute centrality!");
+ return -1.0;
}
return centrality->GetCentralityPercentile(fCentralityType.Data());
- }
+ }
} else {
if (!fCentralityType.CompareTo("TRACKS"))
return fInputEvent->GetNumberOfTracks();
AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", type.Data()));
return -1.0;
}
-
- return 1E20;
}
//__________________________________________________________________________________________________
if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
(TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
(TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
+ if(fCheckFeedDown){
+ Int_t pdgGranma = 0;
+ Int_t mother = 0;
+ mother = part->GetMother();
+ Int_t istep = 0;
+ Int_t abspdgGranma =0;
+ Bool_t isFromB=kFALSE;
+ Bool_t isQuarkFound=kFALSE;
+ while (mother >=0 ){
+ istep++;
+ AliDebug(2,Form("mother at step %d = %d", istep, mother));
+ AliMCParticle* mcGranma = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(mother));
+ if (mcGranma){
+ pdgGranma = mcGranma->PdgCode();
+ AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
+ abspdgGranma = TMath::Abs(pdgGranma);
+ if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
+ isFromB=kTRUE;
+ }
+ if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
+ mother = mcGranma->GetMother();
+ }else{
+ AliError("Failed casting the mother particle!");
+ break;
+ }
+ }
+ if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
+ if(isFromB){
+ if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
+ }
+ else{
+ if (fKeepDfromBOnly) pdgGranma = -999;
+
+ if (pdgGranma == -99999){
+ AliDebug(2,"This particle does not have a quark in his genealogy\n");
+ continue;
+ }
+ if (pdgGranma == -9999){
+ AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
+ continue;
+ }
+
+ if (pdgGranma == -999){
+ AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
+ continue;
+ }
+
+ }
+ }
// assign momenta to computation object
miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
miniPair.FillRef(def->GetMotherMass());
if (!okMatch) continue;
if(fCheckP && (TMath::Abs(part->Px()-(daughter1->Px()+daughter2->Px()))/(TMath::Abs(part->Px())+1.e-13)) > 0.00001 &&
(TMath::Abs(part->Py()-(daughter1->Py()+daughter2->Py()))/(TMath::Abs(part->Py())+1.e-13)) > 0.00001 &&
- (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
-
+ (TMath::Abs(part->Pz()-(daughter1->Pz()+daughter2->Pz()))/(TMath::Abs(part->Pz())+1.e-13)) > 0.00001 ) continue;
+ if(fCheckFeedDown){
+ Int_t pdgGranma = 0;
+ Int_t mother = 0;
+ mother = part->GetMother();
+ Int_t istep = 0;
+ Int_t abspdgGranma =0;
+ Bool_t isFromB=kFALSE;
+ Bool_t isQuarkFound=kFALSE;
+ while (mother >=0 ){
+ istep++;
+ AliDebug(2,Form("mother at step %d = %d", istep, mother));
+ AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(list->At(mother));
+ if (mcGranma){
+ pdgGranma = mcGranma->GetPdgCode();
+ AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
+ abspdgGranma = TMath::Abs(pdgGranma);
+ if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
+ isFromB=kTRUE;
+ }
+ if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
+ mother = mcGranma->GetMother();
+ }else{
+ AliError("Failed casting the mother particle!");
+ break;
+ }
+ }
+ if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
+ if(isFromB){
+ if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
+ }
+ else{
+ if (fKeepDfromBOnly) pdgGranma = -999;
+ }
+
+ if (pdgGranma == -99999){
+ AliDebug(2,"This particle does not have a quark in his genealogy\n");
+ continue;
+ }
+ if (pdgGranma == -9999){
+ AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
+ continue;
+ }
+
+ if (pdgGranma == -999){
+ AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
+ continue;
+ }
+ }
// assign momenta to computation object
miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2);
miniPair.FillRef(def->GetMotherMass());
}
}
}
-
+//___________________________________________________________
+void AliRsnMiniAnalysisTask::SetDselection(UShort_t originDselection)
+{
+ // setting the way the D0 will be selected
+ // 0 --> only from c quarks
+ // 1 --> only from b quarks
+ // 2 --> from both c quarks and b quarks
+
+ fOriginDselection = originDselection;
+
+ if (fOriginDselection == 0) {
+ fKeepDfromB = kFALSE;
+ fKeepDfromBOnly = kFALSE;
+ }
+
+ if (fOriginDselection == 1) {
+ fKeepDfromB = kTRUE;
+ fKeepDfromBOnly = kTRUE;
+ }
+
+ if (fOriginDselection == 2) {
+ fKeepDfromB = kTRUE;
+ fKeepDfromBOnly = kFALSE;
+ }
+
+ return;
+}
//__________________________________________________________________________________________________
Bool_t AliRsnMiniAnalysisTask::EventsMatch(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2)
{
}
}
+//---------------------------------------------------------------------
+Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchPbPb2011(){
+ //This part rejects randomly events such that the centrality gets flat for LHC11h Pb-Pb data
+ //for 0-5% and 10-20% centrality bin
+
+ if (fCentralityType!="V0M") {
+ AliWarning("Wrong value (not centrality from V0).");
+ return -999.0;
+ }
+
+ AliCentrality *centrality = fInputEvent->GetCentrality();
+ if (!centrality) {
+ AliWarning("Cannot get centrality from AOD event.");
+ return -999.0;
+ }
+
+ Double_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
+ Double_t rnd_hc = -1., testf = 0.0, ff = 0, N1 = -1., N2 = -1.;
+
+ if(fUseCentralityPatchPbPb2011==510){
+ N1 = 1.9404e+06;
+ N2 = 1.56435e+06; //N2 is the reference
+ ff = 5.04167e+06 - 1.49885e+06*cent + 2.35998e+05*cent*cent -1.22873e+04*cent*cent*cent;
+ } else {
+ if(fUseCentralityPatchPbPb2011==1020){
+ N2 = 2.0e+05; //N2 is the reference
+ N1 = 3.7e+05;
+ ff = -1.73979e+06 - 3.05316e+06*cent + 1.05517e+06*cent*cent - 133205*cent*cent*cent + 8187.45*cent*cent*cent*cent - 247.875*cent*cent*cent*cent*cent + 2.9676*cent*cent*cent*cent*cent*cent;
+ } else {
+ AliError(Form("Patch for the requested centrality (%i) is not available", fUseCentralityPatchPbPb2011));
+ return -999.0;
+ }
+ }
+ testf = ( N2 + (N1-ff) ) / N1;
+ rnd_hc = gRandom->Rndm();
+ //AliDebugClass(1, Form("Flat Centrality %d", fUseCentralityPatchPbPb2011));
+
+ if (rnd_hc < 0 || rnd_hc > 1 )
+ {
+ AliWarning("Wrong Random number generated");
+ return -999.0;
+ }
+
+ if (rnd_hc < testf)
+ return cent;
+ else
+ return -999.0;
+}
//---------------------------------------------------------------------
Double_t AliRsnMiniAnalysisTask::ApplyCentralityPatchAOD049()
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff