X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=PWGHF%2FvertexingHF%2FAliAnalysisTaskSELc2V0bachelorTMVA.cxx;h=6ba1ace0a9997aba1a5c0bdd2f4986f3a9827ac1;hb=041f8e44a1a06a5ee06451abb7786409bffe4ce2;hp=da29ba2d28c46c030f2306ec495281754c646409;hpb=0979b1fe1fe85e9090c3581a182d6052ecf072f3;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx b/PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx old mode 100755 new mode 100644 index da29ba2d28c..6ba1ace0a99 --- a/PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx +++ b/PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx @@ -187,7 +187,10 @@ AliAnalysisTaskSE(), fCutKFDeviationFromVtxV0(0.), fCurrentEvent(-1), fBField(0), - fKeepingOnlyPYTHIABkg(kFALSE) + fKeepingOnlyPYTHIABkg(kFALSE), + fHistoMCLcK0SpGen(0x0), + fHistoMCLcK0SpGenAcc(0x0), + fHistoMCLcK0SpGenLimAcc(0x0) { // // Default ctor @@ -307,7 +310,10 @@ AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Cha fCutKFDeviationFromVtxV0(0.), fCurrentEvent(-1), fBField(0), - fKeepingOnlyPYTHIABkg(kFALSE) + fKeepingOnlyPYTHIABkg(kFALSE), + fHistoMCLcK0SpGen(0x0), + fHistoMCLcK0SpGenAcc(0x0), + fHistoMCLcK0SpGenLimAcc(0x0) { // @@ -407,7 +413,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*) // a query. It always runs on the client, it can be used to present // the results graphically or save the results to file. - //AliInfo("Terminate",""); + AliInfo("Terminate"); AliAnalysisTaskSE::Terminate(); fOutput = dynamic_cast (GetOutputData(1)); @@ -415,6 +421,27 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*) AliError("fOutput not available"); return; } + + + //AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries())); + //AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries())); + //AliDebug(2, Form("At Reco level, %lld Lc --> K0S + p were found", fVariablesTreeSgn->GetEntries())); + if(fHistoMCLcK0SpGen) { + AliInfo(Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries())); + } else { + AliInfo("fHistoMCLcK0SpGen not available"); + } + if(fHistoMCLcK0SpGenAcc) { + AliInfo(Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries())); + } else { + AliInfo("fHistoMCLcK0SpGenAcc not available"); + } + if(fVariablesTreeSgn) { + AliInfo(Form("At Reco level, %lld Lc --> K0S + p were found", fVariablesTreeSgn->GetEntries())); + } else { + AliInfo("fVariablesTreeSgn not available"); + } + fOutputKF = dynamic_cast (GetOutputData(6)); if (!fOutputKF) { AliError("fOutputKF not available"); @@ -605,6 +632,12 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() { //fOutput->Add(fVariablesTreeSgn); //fOutput->Add(fVariablesTreeBkg); + const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.}; + + fHistoMCLcK0SpGen = new TH1F("fHistoMCLcK0SpGen", "fHistoMCLcK0SpGen", 14, ptbins); + fHistoMCLcK0SpGenAcc = new TH1F("fHistoMCLcK0SpGenAcc", "fHistoMCLcK0SpGenAcc", 14, ptbins); + fHistoMCLcK0SpGenLimAcc = new TH1F("fHistoMCLcK0SpGenLimAcc", "fHistoMCLcK0SpGenLimAcc", 14, ptbins); + fOutput->Add(fHistoEvents); fOutput->Add(fHistoLc); fOutput->Add(fHistoLcOnTheFly); @@ -614,6 +647,9 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() { fOutput->Add(fHistoCodesBkg); fOutput->Add(fHistoLcpKpiBeforeCuts); fOutput->Add(fHistoBackground); + fOutput->Add(fHistoMCLcK0SpGen); + fOutput->Add(fHistoMCLcK0SpGenAcc); + fOutput->Add(fHistoMCLcK0SpGenLimAcc); PostData(1, fOutput); PostData(4, fVariablesTreeSgn); @@ -824,7 +860,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *) } fCurrentEvent++; - Printf("Processing event = %d", fCurrentEvent); + AliDebug(2, Form("Processing event = %d", fCurrentEvent)); AliAODEvent* aodEvent = dynamic_cast(fInputEvent); TClonesArray *arrayLctopKos=0; @@ -856,32 +892,21 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *) fAnalCuts->SetTriggerClass(""); fAnalCuts->SetTriggerMask(AliVEvent::kINT7); } - - fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo); - - // AOD primary vertex - fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex(); - if (!fVtx1) return; - if (fVtx1->GetNContributors()<1) return; - - fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent); - - if ( !fIsEventSelected ) { - fHistoEvents->Fill(0); - return; // don't take into account not selected events + + Int_t runnumber = aodEvent->GetRunNumber(); + if (aodEvent->GetTriggerMask() == 0 && (runnumber >= 195344 && runnumber <= 195677)){ + AliDebug(3,"Event rejected because of null trigger mask"); + return; } - fHistoEvents->Fill(1); - - // Setting magnetic field for KF vertexing - fBField = aodEvent->GetMagneticField(); - AliKFParticle::SetField(fBField); - + + fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo); + // mc analysis TClonesArray *mcArray = 0; AliAODMCHeader *mcHeader=0; if (fUseMCInfo) { - // MC array need for maching + // MC array need for matching mcArray = dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); if (!mcArray) { AliError("Could not find Monte-Carlo in AOD"); @@ -893,7 +918,34 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *) AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n"); return; } + + Double_t zMCVertex = mcHeader->GetVtxZ(); + if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()){ + AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fAnalCuts->GetMaxVtxZ(), fAnalCuts->GetMaxVtxZ())); + AliInfo(Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fAnalCuts->GetMaxVtxZ(), fAnalCuts->GetMaxVtxZ())); + return; + } + + //Printf("Filling MC histo"); + FillMCHisto(mcArray); + } + + // AOD primary vertex + fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex(); + if (!fVtx1) return; + if (fVtx1->GetNContributors()<1) return; + + fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent); + + if ( !fIsEventSelected ) { + fHistoEvents->Fill(0); + return; // don't take into account not selected events } + fHistoEvents->Fill(1); + + // Setting magnetic field for KF vertexing + fBField = aodEvent->GetMagneticField(); + AliKFParticle::SetField(fBField); Int_t nSelectedAnal = 0; if (fIsK0sAnalysis) { @@ -909,6 +961,112 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *) PostData(5, fVariablesTreeBkg); PostData(6, fOutputKF); +} +//------------------------------------------------------------------------------- +void AliAnalysisTaskSELc2V0bachelorTMVA::FillMCHisto(TClonesArray *mcArray){ + + // method to fill MC histo: how many Lc --> K0S + p are there at MC level + for (Int_t iPart=0; iPartGetEntriesFast(); iPart++) { + AliAODMCParticle* mcPart = dynamic_cast(mcArray->At(iPart)); + if (!mcPart){ + AliError("Failed casting particle from MC array!, Skipping particle"); + continue; + } + Int_t pdg = mcPart->GetPdgCode(); + if (TMath::Abs(pdg) != 4122){ + AliDebug(2, Form("MC particle %d is not a Lc: its pdg code is %d", iPart, pdg)); + continue; + } + AliDebug(2, Form("Step 0 ok: MC particle %d is a Lc: its pdg code is %d", iPart, pdg)); + Int_t labeldaugh0 = mcPart->GetDaughter(0); + Int_t labeldaugh1 = mcPart->GetDaughter(1); + if (labeldaugh0 <= 0 || labeldaugh1 <= 0){ + AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!")); + continue; + } + else if (labeldaugh1 - labeldaugh0 == 1){ + AliDebug(2, Form("Step 1 ok: The MC particle has correct daughters!!")); + AliAODMCParticle* daugh0 = dynamic_cast(mcArray->At(labeldaugh0)); + AliAODMCParticle* daugh1 = dynamic_cast(mcArray->At(labeldaugh1)); + if(!daugh0 || !daugh1){ + AliDebug(2,"Particle daughters not properly retrieved!"); + return; + } + Int_t pdgCodeDaugh0 = TMath::Abs(daugh0->GetPdgCode()); + Int_t pdgCodeDaugh1 = TMath::Abs(daugh1->GetPdgCode()); + AliAODMCParticle* bachelorMC = daugh0; + AliAODMCParticle* v0MC = daugh1; + AliDebug(2, Form("pdgCodeDaugh0 = %d, pdgCodeDaugh1 = %d", pdgCodeDaugh0, pdgCodeDaugh1)); + if ((pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) || (pdgCodeDaugh0 == 2212 && pdgCodeDaugh1 == 311)){ + // we are in the case of Lc --> K0 + p; now we have to check if the K0 decays in K0S, and if this goes in pi+pi- + /// first, we set the bachelor and the v0: above we assumed first proton and second V0, but we could have to change it: + if (pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) { + bachelorMC = daugh1; + v0MC = daugh0; + } + AliDebug(2, Form("Number of Daughters of v0 = %d", v0MC->GetNDaughters())); + if (v0MC->GetNDaughters() != 1) { + AliDebug(2, "The K0 does not decay in 1 body only! Impossible... Continuing..."); + continue; + } + else { // So far: Lc --> K0 + p, K0 with 1 daughter + AliDebug(2, "Step 2 ok: The K0 does decay in 1 body only! "); + Int_t labelK0daugh = v0MC->GetDaughter(0); + AliAODMCParticle* partK0S = dynamic_cast(mcArray->At(labelK0daugh)); + if(!partK0S){ + AliError("Error while casting particle! returning a NULL array"); + continue; + } + else { // So far: Lc --> K0 + p, K0 with 1 daughter that we can access + if (partK0S->GetNDaughters() != 2 || TMath::Abs(partK0S->GetPdgCode() != 310)){ + AliDebug(2, "The K0 daughter is not a K0S or does not decay in 2 bodies"); + continue; + } + else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies + AliDebug(2, "Step 3 ok: The K0 daughter is a K0S and does decay in 2 bodies"); + Int_t labelK0Sdaugh0 = partK0S->GetDaughter(0); + Int_t labelK0Sdaugh1 = partK0S->GetDaughter(1); + AliAODMCParticle* daughK0S0 = dynamic_cast(mcArray->At(labelK0Sdaugh0)); + AliAODMCParticle* daughK0S1 = dynamic_cast(mcArray->At(labelK0Sdaugh1)); + if (!daughK0S0 || ! daughK0S1){ + AliDebug(2, "Could not access K0S daughters, continuing..."); + continue; + } + else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies that we can access + AliDebug(2, "Step 4 ok: Could access K0S daughters, continuing..."); + Int_t pdgK0Sdaugh0 = daughK0S0->GetPdgCode(); + Int_t pdgK0Sdaugh1 = daughK0S1->GetPdgCode(); + if (TMath::Abs(pdgK0Sdaugh0) != 211 || TMath::Abs(pdgK0Sdaugh1) != 211){ + AliDebug(2, "The K0S does not decay in pi+pi-, continuing"); + //AliInfo("The K0S does not decay in pi+pi-, continuing"); + } + else { // Full chain: Lc --> K0 + p, K0 --> K0S, K0S --> pi+pi- + if (fAnalCuts->IsInFiducialAcceptance(mcPart->Pt(), mcPart->Y())) { + AliDebug(2, Form("----> Filling histo with pt = %f", mcPart->Pt())); + if(TMath::Abs(mcPart->Y()) < 0.5) fHistoMCLcK0SpGenLimAcc->Fill(mcPart->Pt()); + //AliInfo(Form("\nparticle = %d, Filling MC Gen histo\n", iPart)); + fHistoMCLcK0SpGen->Fill(mcPart->Pt()); + if(!(TMath::Abs(bachelorMC->Eta()) > 0.9 || bachelorMC->Pt() < 0.1 || + TMath::Abs(daughK0S0->Eta()) > 0.9 || daughK0S0->Pt() < 0.1 || + TMath::Abs(daughK0S1->Eta()) > 0.9 || daughK0S1->Pt() < 0.1)) { + fHistoMCLcK0SpGenAcc->Fill(mcPart->Pt()); + } + } + else { + AliDebug(2, "not in fiducial acceptance! Skipping"); + continue; + } + } + } + } + } + } + } + } + } // closing loop over mcArray + + return; + } //------------------------------------------------------------------------------- @@ -976,7 +1134,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *a // find associated MC particle for Lc -> p+K0 and K0S->pi+pi fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE); if (fmcLabelLc>=0) { - AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades)); + AliDebug(2, Form("----> cascade number %d (total cascade number = %d) is a Lc!", iLctopK0s, nCascades)); AliAODMCParticle *partLc = dynamic_cast(mcArray->At(fmcLabelLc)); if(partLc){ @@ -1023,7 +1181,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *a AliAODTrack * v0Pos = dynamic_cast(lcK0spr->Getv0PositiveTrack()); AliAODTrack * v0Neg = dynamic_cast(lcK0spr->Getv0NegativeTrack()); - if (!v0Neg || !v0Neg) { + if (!v0Neg || !v0Pos) { AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s)); continue; } @@ -1119,7 +1277,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *a } } - FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray); + FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray, iLctopK0s); } @@ -1131,7 +1289,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF Int_t isLc, Int_t &nSelectedAnal, AliRDHFCutsLctoV0 *cutsAnal, - TClonesArray *mcArray){ + TClonesArray *mcArray, Int_t iLctopK0s){ // // Fill histos for Lc -> K0S+proton // @@ -1189,8 +1347,8 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass if (isInV0window == 0) { - AliDebug(2, "No: The candidate has NOT passed the mass cuts!"); - if (isLc) Printf("SIGNAL candidate rejected"); + AliDebug(2, "No: The candidate has NOT passed the V0 window cuts!"); + if (isLc) Printf("SIGNAL candidate rejected: V0 window cuts"); return; } else AliDebug(2, "Yes: The candidate has passed the mass cuts!"); @@ -1199,12 +1357,13 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF if (!isInCascadeWindow) { AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!"); - if (isLc) Printf("SIGNAL candidate rejected"); + if (isLc) Printf("SIGNAL candidate rejected: cascade window cuts"); return; } else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!"); Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts + AliDebug(2, Form("recoAnalysisCuts = %d", cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate) & (AliRDHFCutsLctoV0::kLcToK0Spr))); if (!isCandidateSelectedCuts){ AliDebug(2, "No: Analysis cuts kCandidate level NOT passed"); if (isLc) Printf("SIGNAL candidate rejected"); @@ -1226,13 +1385,13 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF); AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed)); Double_t probProton = -1.; - Double_t probPion = -1.; - Double_t probKaon = -1.; + // Double_t probPion = -1.; + // Double_t probKaon = -1.; if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) { AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton])); probProton = probTPCTOF[AliPID::kProton]; - probPion = probTPCTOF[AliPID::kPion]; - probKaon = probTPCTOF[AliPID::kKaon]; + // probPion = probTPCTOF[AliPID::kPion]; + // probKaon = probTPCTOF[AliPID::kKaon]; } else { // if you don't have both TOF and TPC, try only TPC fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC); @@ -1241,8 +1400,8 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF AliDebug(2,Form(" detUsed (TPC case) = %d", detUsed)); if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) { probProton = probTPCTOF[AliPID::kProton]; - probPion = probTPCTOF[AliPID::kPion]; - probKaon = probTPCTOF[AliPID::kKaon]; + // probPion = probTPCTOF[AliPID::kPion]; + // probKaon = probTPCTOF[AliPID::kKaon]; AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton])); } else { @@ -1543,7 +1702,7 @@ void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF if (fUseMCInfo) { if (isLc){ - AliDebug(2, "Filling Sgn"); + AliDebug(2, Form("Reco particle %d --> Filling Sgn", iLctopK0s)); fVariablesTreeSgn->Fill(); fHistoCodesSgn->Fill(bachCode, k0SCode); } @@ -1637,9 +1796,9 @@ Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *c } } - Double_t xn, xp, dca; + Double_t xn=0., xp=0.;//, dca; AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2)); - dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp); + // dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp); AliExternalTrackParam tr1(*esdv0Daugh1); AliExternalTrackParam tr2(*esdv0Daugh2); @@ -1717,11 +1876,13 @@ Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *c if (!tmpdaughv02 && labelsv0daugh[1] > 0){ AliDebug(2, "Could not access MC info for second daughter of V0, continuing"); } - Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays - Double_t yPionMC = tmpdaughv01->Yv(); - Double_t zPionMC = tmpdaughv01->Zv(); - //Printf("Got MC vtx for Pion"); - Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC); + if(tmpdaughv01){ + Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays + Double_t yPionMC = tmpdaughv01->Yv(); + Double_t zPionMC = tmpdaughv01->Zv(); + //Printf("Got MC vtx for Pion"); + Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC); + } } else { Printf("Not a true V0");