MFT analysis class
authorauras <antonio.uras@cern.ch>
Tue, 16 Dec 2014 16:44:36 +0000 (17:44 +0100)
committerauras <antonio.uras@cern.ch>
Tue, 16 Dec 2014 16:44:36 +0000 (17:44 +0100)
PWG/mftmuondep/CMakeLists.txt
PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.cxx [new file with mode: 0644]
PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.h [new file with mode: 0755]
PWG/mftmuondep/charmonium/CMakeLists.txt
PWG/mftmuondep/charmonium/PWGmftmuondepCharmoniumLinkDef.h [new file with mode: 0644]
PWG/mftmuondep/charmonium/RunAnalysisTaskDimuonBackground.C [new file with mode: 0755]

index ceda7e2..ac5ba0b 100644 (file)
@@ -15,6 +15,6 @@
 
 # Include the libraries
 add_subdirectory(charmonium)
-add_subdirectory(openHF)
+#add_subdirectory(openHF)
 
 message(STATUS "PWG mftmuondep enabled")
diff --git a/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.cxx b/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.cxx
new file mode 100644 (file)
index 0000000..03a5959
--- /dev/null
@@ -0,0 +1,525 @@
+#include "TH1D.h"
+#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliAnalysisTaskDimuonBackground.h"
+#include "TDatabasePDG.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
+#include "AliAODDimuon.h"
+#include "AliAODTrack.h"
+#include "AliAODHeader.h"
+#include "TClonesArray.h"
+#include "AliMFTConstants.h"
+#include "AliMFTAnalysisTools.h"
+#include "TRandom.h"
+#include "TList.h"
+
+ClassImp(AliAnalysisTaskDimuonBackground)
+
+//====================================================================================================================================================
+
+AliAnalysisTaskDimuonBackground::AliAnalysisTaskDimuonBackground() : 
+  AliAnalysisTaskSE(),
+  fVertexMode(0),
+  fMinTriggerMatch(0),
+  fSingleMuonMinEta(-9999), 
+  fSingleMuonMaxEta(9999), 
+  fSingleMuonMinPt(0),
+  fSingleMuonMaxChi2(9999),
+  fHistogramList(0),
+  fMainInputHandler(0),
+  fMixingInputHandler(0)
+{
+  
+  // Default constructor
+
+  fHistSingleMuonsChi2VsPtVsRapidityBeforeCut = 0;
+  fHistSingleMuonsChi2VsPtVsRapidityAfterCut = 0;
+
+  for (Int_t i=0; i<2; i++) fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] = 0;
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances = 0;
+  
+  for (Int_t i=0; i<3; i++) { fPrimaryVertex[i]=0; fPrimaryVertexMixEv[i]=0; fPrimaryVertexTrue[i]=0; fPrimaryVertexMixEvTrue[i]=0; }
+  fVtxResolutionITS[0] = 5.e-4;
+  fVtxResolutionITS[1] = 5.e-4;
+  fVtxResolutionITS[2] = 4.e-4;
+
+  fMassJpsi = TDatabasePDG::Instance()->GetParticle(443)->Mass();
+
+}
+
+//====================================================================================================================================================
+
+AliAnalysisTaskDimuonBackground::AliAnalysisTaskDimuonBackground(const Char_t *name) : 
+  AliAnalysisTaskSE(name),
+  fVertexMode(0),
+  fMinTriggerMatch(0),
+  fSingleMuonMinEta(-9999), 
+  fSingleMuonMaxEta(9999), 
+  fSingleMuonMinPt(0),
+  fSingleMuonMaxChi2(9999),
+  fHistogramList(0),
+  fMainInputHandler(0),
+  fMixingInputHandler(0) 
+{
+
+  // Constructor
+
+  fHistSingleMuonsChi2VsPtVsRapidityBeforeCut = 0;
+  fHistSingleMuonsChi2VsPtVsRapidityAfterCut = 0;
+
+  for (Int_t i=0; i<2; i++) fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] = 0;
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances = 0;
+
+  for (Int_t i=0; i<3; i++) { fPrimaryVertex[i]=0; fPrimaryVertexMixEv[i]=0; fPrimaryVertexTrue[i]=0; fPrimaryVertexMixEvTrue[i]=0; }
+  fVtxResolutionITS[0] = 5.e-4;
+  fVtxResolutionITS[1] = 5.e-4;
+  fVtxResolutionITS[2] = 4.e-4;
+
+  fMassJpsi = TDatabasePDG::Instance()->GetParticle(443)->Mass();
+
+  // Define input and output slots here
+  DefineOutput(1, TList::Class());
+
+}
+
+//====================================================================================================================================================
+
+void AliAnalysisTaskDimuonBackground::UserCreateOutputObjects() {
+
+  // Called once
+
+  fHistogramList = new TList();
+  fHistogramList->SetOwner(kTRUE);
+
+  // Max 100 dimensions
+  Int_t  nBins[100] = {0};
+  Double_t min[100] = {0};
+  Double_t max[100] = {0};
+
+  TString tag[2] = {"SingleEvents", "MixedEvents"};
+
+  //-------------------------------------------------------------------------------
+
+  nBins[0] = 100;  
+  nBins[1] = 1000;
+  nBins[2] = 1000;
+  nBins[3] = 100;
+  nBins[4] = 100;
+
+  min[0] = 0.;
+  min[1] = -10.;
+  min[2] = 0.;
+  min[3] = 0.;
+  min[4] = -4.5;
+
+  max[0] = 1.;
+  max[1] = 10.;
+  max[2] = 10.;
+  max[3] = 10.;
+  max[4] = -2.0;
+
+  for (Int_t i=0; i<2; i++) {
+    fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] = new THnSparseD(Form("fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_%s",tag[i].Data()), "", 5, nBins, min, max);
+    fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> Sumw2(); 
+    fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(0)->SetTitle("PCA Quality");
+    fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(1)->SetTitle("t_{z}  [ps]");
+    fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(2)->SetTitle("Mass  [GeV/c^{2}]");
+    fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(3)->SetTitle("p_{T}  [GeV/c]");
+    fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(4)->SetTitle("Rapidity");
+    fHistogramList->Add(fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i]);
+  }
+
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances = new THnSparseD("fHistPCAQualityVsPPDecayTimeVsMassVsPt_Resonances", "", 5, nBins, min, max);
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> Sumw2(); 
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(0)->SetTitle("PCA Quality");
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(1)->SetTitle("t_{z}  [ps]");
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(2)->SetTitle("Mass  [GeV/c^{2}]");
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(3)->SetTitle("p_{T}  [GeV/c]");
+  fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(4)->SetTitle("Rapidity");
+  fHistogramList->Add(fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances);
+
+  nBins[0] = 200;  
+  nBins[1] = 100;
+  nBins[2] = 100;
+
+  min[0] = 0.;
+  min[1] = 0.;
+  min[2] = -4.5;
+
+  max[0] = 10.;
+  max[1] = 10.;
+  max[2] = -2.0;
+  
+  fHistSingleMuonsChi2VsPtVsRapidityBeforeCut = new THnSparseD("fHistSingleMuonsChi2VsPtVsRapidityBeforeCut", "", 3, nBins, min, max);
+  fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> Sumw2(); 
+  fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> GetAxis(0)->SetTitle("p_{T}  [GeV/c]");
+  fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> GetAxis(1)->SetTitle("#chi^{2}");
+  fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> GetAxis(2)->SetTitle("Rapidity");
+  fHistogramList->Add(fHistSingleMuonsChi2VsPtVsRapidityBeforeCut);
+  
+  fHistSingleMuonsChi2VsPtVsRapidityAfterCut = new THnSparseD("fHistSingleMuonsChi2VsPtVsRapidityAfterCut", "", 3, nBins, min, max);
+  fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> Sumw2(); 
+  fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> GetAxis(0)->SetTitle("p_{T}  [GeV/c]");
+  fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> GetAxis(1)->SetTitle("#chi^{2}");
+  fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> GetAxis(2)->SetTitle("Rapidity");
+  fHistogramList->Add(fHistSingleMuonsChi2VsPtVsRapidityAfterCut);
+
+  //-------------------------------------------------------------------------------
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  SetMainInputHandler(mgr);
+  if (fMainInputHandler) SetMixingInputHandler(fMainInputHandler);
+
+  PostData(1, fHistogramList);
+
+}
+
+//====================================================================================================================================================
+
+void AliAnalysisTaskDimuonBackground::UserExec(Option_t *) {
+
+  // Main loop
+  // Called for each event
+
+  AliAODEvent *aodEv = dynamic_cast<AliAODEvent*>(GetMainEvent());
+
+  if (!aodEv) return;
+
+  if (!(AliMUONTrackExtrap::IsFieldON())) {
+    ((AliAODHeader*) aodEv->GetHeader())->InitMagneticField();
+    AliMUONTrackExtrap::SetField();
+  }
+
+  // Getting primary vertex, either from the generation or from the reconstruction -------------------
+
+  AliAODMCHeader *mcHeader = (AliAODMCHeader*) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+
+  if (fVertexMode == kGenerated) {
+    mcHeader->GetVertex(fPrimaryVertex);
+    for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = gRandom->Gaus(fPrimaryVertex[i], fVtxResolutionITS[i]);
+  }
+  else if (fVertexMode == kReconstructed) {
+    aodEv->GetPrimaryVertex()->GetXYZ(fPrimaryVertex);
+  }
+
+  TClonesArray *stackMC = (TClonesArray*) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
+
+  AliAODTrack *recMuon[2] = {0};
+  AliAODMCParticle *mcMuon[2]={0}, *mcMother[2]={0};
+  Double_t var[100] = {0};
+
+  //--------------------------------------------------------------------------------
+
+  for (Int_t iTrack=0; iTrack<aodEv->GetNumberOfTracks(); iTrack++) { 
+
+    recMuon[0] = (AliAODTrack*) aodEv->GetTrack(iTrack);
+    if (!(recMuon[0]->IsMuonGlobalTrack()))                                  continue;
+    if (AliMFTAnalysisTools::IsTrackInjected(recMuon[0], mcHeader, stackMC)) continue;    // Only HIJING tracks are considered to build the background
+
+    var[0] = recMuon[0]->Pt();
+    var[1] = recMuon[0]->Chi2perNDF();
+    var[2] = recMuon[0]->Y();
+
+    fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> Fill(var);
+
+    if (!IsSingleMuonCutPassed(recMuon[0]))                                  continue;
+
+    fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> Fill(var);
+
+    mcMuon[0]   = NULL;
+    mcMother[0] = NULL;
+
+    if (recMuon[0]->GetLabel()>=0)              mcMuon[0]   = (AliAODMCParticle*) stackMC->At(recMuon[0]->GetLabel());
+    if (mcMuon[0] && mcMuon[0]->GetMother()>=0) mcMother[0] = (AliAODMCParticle*) stackMC->At(mcMuon[0]->GetMother());
+
+    for (Int_t jTrack=0; jTrack<iTrack; jTrack++) { 
+
+      recMuon[1] = (AliAODTrack*) aodEv->GetTrack(jTrack);
+      if (!(recMuon[1]->IsMuonGlobalTrack()))                                  continue;
+      if (AliMFTAnalysisTools::IsTrackInjected(recMuon[1], mcHeader, stackMC)) continue;    // Only HIJING tracks are considered to build the background
+      if (!IsSingleMuonCutPassed(recMuon[1]))                                  continue;
+
+      mcMuon[1]   = NULL;
+      mcMother[1] = NULL;
+      
+      if (recMuon[1]->GetLabel()>=0)              mcMuon[1]   = (AliAODMCParticle*) stackMC->At(recMuon[1]->GetLabel());
+      if (mcMuon[1] && mcMuon[1]->GetMother()>=0) mcMother[1] = (AliAODMCParticle*) stackMC->At(mcMuon[1]->GetMother());
+      
+      TObjArray *dimuon = new TObjArray();
+      dimuon -> Add(recMuon[0]);
+      dimuon -> Add(recMuon[1]);       
+      
+      if (recMuon[0]->Charge() + recMuon[1]->Charge()) {
+       delete dimuon;
+       continue;
+      }
+
+      Bool_t isDimuonResonance = kFALSE;
+      if (mcMother[0] && mcMother[1] && mcMuon[0]->GetMother()==mcMuon[1]->GetMother()) {
+       if (AliMFTAnalysisTools::IsPDGResonance(mcMother[0]->GetPdgCode())) isDimuonResonance = kTRUE;
+      }
+
+      Double_t pca[3]={0};
+      Double_t pcaQuality=0;
+      TLorentzVector kinem(0,0,0,0);
+      if (!AliMFTAnalysisTools::CalculatePCA(dimuon, pca, pcaQuality, kinem)) {
+       delete dimuon;
+       continue;
+      }
+
+      Double_t ppDecayTime = AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(fPrimaryVertex[2], pca[2], fMassJpsi, kinem.Pz());
+
+      var[0] = pcaQuality;
+      var[1] = ppDecayTime;
+      var[2] = kinem.M();
+      var[3] = kinem.Pt();
+      var[4] = kinem.Rapidity();
+
+      if (!isDimuonResonance) fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[kSingleEvents] -> Fill(var);
+      else                    fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances     -> Fill(var);
+
+      delete dimuon;
+
+    }   // end of loop on 2nd muon
+
+  }   // end of loop on 1st muon
+
+  //--------------------------------------------------------------------------------
+  
+  PostData(1, fHistogramList);
+
+}
+
+//====================================================================================================================================================
+
+void AliAnalysisTaskDimuonBackground::UserExecMix(Option_t *) {
+
+  // Performing Mixing
+
+  if (!fMixingInputHandler) return;
+
+  Int_t bufferSize  = fMixingInputHandler->BufferSize();
+
+  AliAODEvent *aodEv = dynamic_cast<AliAODEvent*>(GetMainEvent());
+
+  if (!aodEv) return;
+
+  if (!(AliMUONTrackExtrap::IsFieldON())) {
+    ((AliAODHeader*) aodEv->GetHeader())->InitMagneticField();
+    AliMUONTrackExtrap::SetField();
+  }
+
+  // Getting primary vertex, either from the generation or from the reconstruction -------------------
+
+  AliAODMCHeader *mcHeader = (AliAODMCHeader*) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+  mcHeader->GetVertex(fPrimaryVertexTrue);
+
+  if (fVertexMode == kGenerated) {
+    mcHeader->GetVertex(fPrimaryVertex);
+    for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = gRandom->Gaus(fPrimaryVertex[i], fVtxResolutionITS[i]);
+  }
+  else if (fVertexMode == kReconstructed) {
+    aodEv->GetPrimaryVertex()->GetXYZ(fPrimaryVertex);
+  }
+
+  TClonesArray *stackMC = (TClonesArray*) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
+
+  // -------------------------------------------------------------------------------------------------
+
+  AliAODTrack *recMuon[2] = {0};
+  Double_t var[100] = {0};
+
+  for (Int_t iBuffer=0; iBuffer<bufferSize; iBuffer++) {
+
+    AliAODEvent *aodEvMix = dynamic_cast<AliAODEvent *>(GetMixedEvent(iBuffer));
+    if (!aodEvMix) continue;
+
+    // Getting primary vertex, either from the generation or from the reconstruction -------------------
+    
+    AliAODMCHeader *mcHeaderMixEv = (AliAODMCHeader*) (aodEvMix->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+    mcHeaderMixEv->GetVertex(fPrimaryVertexMixEvTrue);
+    if (fVertexMode == kGenerated) {
+      mcHeaderMixEv->GetVertex(fPrimaryVertexMixEv);
+      for (Int_t i=0; i<3; i++) fPrimaryVertexMixEv[i] = gRandom->Gaus(fPrimaryVertexMixEv[i], fVtxResolutionITS[i]);
+    }
+    else if (fVertexMode == kReconstructed) {
+      aodEvMix->GetPrimaryVertex()->GetXYZ(fPrimaryVertexMixEv);
+    }
+
+    TClonesArray *stackMCMixEv = (TClonesArray*) (aodEvMix->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
+  
+    // Loop over MUON+MFT muons of 1st and 2nd events
+    
+    //--------------------------------------------------------------------------------------------------
+    
+    for (Int_t iTrack=0; iTrack<aodEv->GetNumberOfTracks(); iTrack++) { 
+      
+      recMuon[0] = (AliAODTrack*) aodEv->GetTrack(iTrack);
+      if (!(recMuon[0]->IsMuonGlobalTrack()))                                  continue;
+      if (AliMFTAnalysisTools::IsTrackInjected(recMuon[0], mcHeader, stackMC)) continue;    // Only HIJING tracks are considered to build the background
+      if (!IsSingleMuonCutPassed(recMuon[0]))                                  continue;
+
+      for (Int_t jTrack=0; jTrack<aodEvMix->GetNumberOfTracks(); jTrack++) { 
+
+       recMuon[1] = (AliAODTrack*) aodEvMix->GetTrack(jTrack);
+       if (!(recMuon[1]->IsMuonGlobalTrack()))                                            continue;
+       if (AliMFTAnalysisTools::IsTrackInjected(recMuon[1], mcHeaderMixEv, stackMCMixEv)) continue;    // Only HIJING tracks are considered to build the background
+       if (!IsSingleMuonCutPassed(recMuon[1]))                                            continue;
+
+       // ----------- Translating muons to a common primary vertex (the origin). Preserve original tracks in case one wants to use them
+
+       AliAODTrack *recMuonTranslated[2] = {0};
+       for (Int_t i=0; i<2; i++) recMuonTranslated[i] = new AliAODTrack(*(recMuon[i]));
+
+       if (!(AliMFTAnalysisTools::TranslateMuonToOrigin(recMuonTranslated[0], fPrimaryVertexTrue)) ||
+           !(AliMFTAnalysisTools::TranslateMuonToOrigin(recMuonTranslated[1], fPrimaryVertexMixEvTrue))) {
+         printf("Error: tracks cannot be translated!!!\n");
+         delete recMuonTranslated[0];
+         delete recMuonTranslated[1];
+         continue;
+       }
+
+       TObjArray *dimuon = new TObjArray();
+       dimuon -> Add(recMuonTranslated[0]);
+       dimuon -> Add(recMuonTranslated[1]);    
+
+       if (recMuonTranslated[0]->Charge() + recMuonTranslated[1]->Charge()) {
+         delete dimuon;
+         delete recMuonTranslated[0];
+         delete recMuonTranslated[1];
+         continue;
+       }
+
+       Double_t pca[3]={0};
+       Double_t pcaQuality=0;
+       TLorentzVector kinem(0,0,0,0);
+       if (!AliMFTAnalysisTools::CalculatePCA(dimuon, pca, pcaQuality, kinem)) {
+         delete dimuon;
+         delete recMuonTranslated[0];
+         delete recMuonTranslated[1];
+         continue;
+       }
+       
+       Double_t ppDecayTime = AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(fPrimaryVertex[2]-fPrimaryVertexTrue[2], pca[2], fMassJpsi, kinem.Pz());
+       
+       var[0] = pcaQuality;
+       var[1] = ppDecayTime;
+       var[2] = kinem.M();
+       var[3] = kinem.Pt();
+       var[4] = kinem.Rapidity();
+
+       fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[kMixedEvents] -> Fill(var);
+       
+       delete dimuon;
+       delete recMuonTranslated[0];
+       delete recMuonTranslated[1];
+       
+      }   // end of loop on 2nd muon
+      
+    }   // end of loop on 1st muon
+    
+    //--------------------------------------------------------------------------------------------------
+   
+  }   // end of loop on events in buffer
+
+  PostData(1, fHistogramList);
+
+}
+
+//====================================================================================================================================================
+
+void AliAnalysisTaskDimuonBackground::Terminate(Option_t *) {
+
+  // Draw result to the screen
+  // Called once at the end of the query
+  
+}
+
+//====================================================================================================================================================
+
+AliVEvent *AliAnalysisTaskDimuonBackground::GetMainEvent() {
+
+  // Access to MainEvent
+  
+  AliMultiInputEventHandler *inEvHMainMulti = fMainInputHandler;
+  if (inEvHMainMulti) {
+    AliInputEventHandler *inEvMain = dynamic_cast<AliInputEventHandler *>(inEvHMainMulti->GetFirstInputEventHandler());
+    if (inEvMain) return inEvMain->GetEvent();
+  }
+  
+  return 0;
+
+}
+
+//====================================================================================================================================================
+
+AliVEvent *AliAnalysisTaskDimuonBackground::GetMixedEvent(Int_t buffId) {
+
+  // Access to Mixed event with buffer id
+  
+  AliMultiInputEventHandler *inEvHMain = fMainInputHandler;
+  
+  if (inEvHMain) {
+    
+    AliMixInputEventHandler *mixIH = fMixingInputHandler;
+    if (!mixIH) return 0;
+
+    if (mixIH->CurrentBinIndex() < 0) {
+      AliDebug(AliLog::kDebug + 1, "Current event mixEH->CurrentEntry() == -1");
+      return 0;
+    }
+    
+    AliMultiInputEventHandler *inEvHMixedCurrent = dynamic_cast<AliMultiInputEventHandler *>(mixIH->InputEventHandler(buffId));
+    if (!inEvHMixedCurrent) return 0;
+
+    AliInputEventHandler *ihMixedCurrent = inEvHMixedCurrent->GetFirstInputEventHandler();
+    if (ihMixedCurrent) return ihMixedCurrent->GetEvent();
+
+  }
+  
+  return 0;
+
+}
+
+//====================================================================================================================================================
+
+AliMultiInputEventHandler *AliAnalysisTaskDimuonBackground::SetMainInputHandler(AliAnalysisManager *mgr) {
+  
+  // Sets main input handler
+  
+  if (!fMainInputHandler) fMainInputHandler = dynamic_cast<AliMultiInputEventHandler *>(mgr->GetInputEventHandler());
+  
+  return fMainInputHandler;
+
+}
+
+//====================================================================================================================================================
+
+AliMixInputEventHandler *AliAnalysisTaskDimuonBackground::SetMixingInputHandler(AliMultiInputEventHandler *mainIH) {
+
+  // Sets mixing input handler
+  
+  if (!fMixingInputHandler) fMixingInputHandler = dynamic_cast<AliMixInputEventHandler *>(mainIH->GetFirstMultiInputHandler());
+  
+  return fMixingInputHandler;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliAnalysisTaskDimuonBackground::IsSingleMuonCutPassed(AliAODTrack *mu) {
+
+  if (mu->GetMatchTrigger() < fMinTriggerMatch)                   return kFALSE;
+  if (mu->Eta()<fSingleMuonMinEta || mu->Eta()>fSingleMuonMaxEta) return kFALSE;
+  if (mu->Pt() < fSingleMuonMinPt)                                return kFALSE;
+  if (mu->Chi2perNDF() > fSingleMuonMaxChi2)                      return kFALSE;
+
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
diff --git a/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.h b/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.h
new file mode 100755 (executable)
index 0000000..d4293bf
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef AliAnalysisTaskDimuonBackground_H
+#define AliAnalysisTaskDimuonBackground_H
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliMultiInputEventHandler.h"
+#include "AliMixInputEventHandler.h"
+#include "TH1D.h"
+#include "TList.h"
+#include "AliAODTrack.h"
+#include "THnSparse.h"
+
+//====================================================================================================================================================
+
+class  AliAnalysisTaskDimuonBackground : public AliAnalysisTaskSE {
+
+public:
+
+  enum {kGenerated, kReconstructed};
+  enum {kSingleEvents, kMixedEvents};
+  AliAnalysisTaskDimuonBackground();
+  AliAnalysisTaskDimuonBackground(const char *name);
+
+  virtual ~AliAnalysisTaskDimuonBackground() { if (fHistogramList) delete fHistogramList; }
+  
+  void SetVertexMode(Int_t vertexMode) { fVertexMode = vertexMode; }
+  void SetVtxResolutionITS(Double_t sigmaX, Double_t sigmaY, Double_t sigmaZ) {
+    fVtxResolutionITS[0] = sigmaX;
+    fVtxResolutionITS[1] = sigmaY;
+    fVtxResolutionITS[2] = sigmaZ;
+  }
+
+  void SetMinTriggerMatch(Int_t minTriggerMatch) { fMinTriggerMatch = minTriggerMatch; }
+
+  void SetSingleMuonMinPt(Double_t minPt)   { fSingleMuonMinPt  = minPt;  }
+  void SetSingleMuonMinEta(Double_t minEta) { fSingleMuonMinEta = minEta; }
+  void SetSingleMuonMaxEta(Double_t maxEta) { fSingleMuonMaxEta = maxEta; }
+  void SetSingleMuonMaxChi2(Double_t maxChi2) { fSingleMuonMaxChi2 = maxChi2; }
+
+  virtual void UserCreateOutputObjects();
+  virtual void UserExec(Option_t *option);
+  virtual void UserExecMix(Option_t *);
+  virtual void Terminate(Option_t *);
+
+  AliVEvent* GetMainEvent();
+  AliVEvent* GetMixedEvent(Int_t buffId=0);
+  
+  AliMultiInputEventHandler* SetMainInputHandler(AliAnalysisManager *mgr);
+  AliMixInputEventHandler* SetMixingInputHandler(AliMultiInputEventHandler *mainIH);
+
+  Bool_t IsSingleMuonCutPassed(AliAODTrack *mu);
+
+private:
+
+  Double_t fMassJpsi;
+
+  Double_t fPrimaryVertexTrue[3], fPrimaryVertexMixEvTrue[3], fPrimaryVertex[3], fPrimaryVertexMixEv[3], fVtxResolutionITS[3];
+  Int_t fVertexMode;
+
+  Int_t fMinTriggerMatch;
+  Double_t fSingleMuonMinEta, fSingleMuonMaxEta, fSingleMuonMinPt, fSingleMuonMaxChi2;
+
+  TList *fHistogramList;
+
+  THnSparse *fHistSingleMuonsChi2VsPtVsRapidityBeforeCut;                   //!
+  THnSparse *fHistSingleMuonsChi2VsPtVsRapidityAfterCut;                    //!
+  THnSparse *fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[2];           //!
+  THnSparse *fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances;   //!
+
+  AliMultiInputEventHandler *fMainInputHandler;    //! tmp pointer to main input handler
+  AliMixInputEventHandler   *fMixingInputHandler;  //! tmp pointer to mixing input handler
+  
+  ClassDef(AliAnalysisTaskDimuonBackground, 1) // example of analysis
+
+};
+
+//====================================================================================================================================================
+
+#endif
index ca37268..c6526d5 100644 (file)
 # * appear in the supporting documentation. The authors make no claims     *
 # * about the suitability of this software for any purpose. It is          *
 # * provided "as is" without express or implied warranty.                  *
-# **************************************************************************/
+# **************************************************************************
 
-#Module
-set(MODULE PWGSomething)
+# Module
+set (MODULE PWGmftmuondepCharmonium)
 
 # Module include folder
 include_directories(${AliRoot_SOURCE_DIR}/PWG/mftmuondep/charmonium)
 
-# Additional includes - alphabetical order except ROOT
-include_directories(${ROOT_INCLUDE_DIRS}
-  )
-
 # Sources - alphabetical order
 set(SRCS
+  AliAnalysisTaskDimuonBackground.cxx
   )
 
 # Headers from sources
 string(REPLACE ".cxx" ".h" HDRS "${SRCS}")
 
-# Generate the dictionary
-# It will create G_ARG1.cxx and G_ARG1.h / ARG1 = function first argument
-get_directory_property(incdirs INCLUDE_DIRECTORIES)
-generate_dictionary("${MODULE}" "${MODULE}LinkDef.h" "${HDRS}" "${incdirs}")
-
-# Add a shared library
-add_library(${MODULE} SHARED  ${SRCS} G__${MODULE}.cxx)
-
-# Generate the ROOT map
-# Dependecies
-set(LIBDEPS)
-generate_rootmap("${MODULE}" "${LIBDEPS}" "${CMAKE_CURRENT_SOURCE_DIR}/${MODULE}LinkDef.h")
-
-# Linking the library
-target_link_libraries(${MODULE} ${LIBDEPS})
-
-# Public include folders that will be propagated to the dependecies
-target_include_directories(${MODULE} PUBLIC ${incdirs})
-
-# System dependent: Modify the way the library is build
-if(${CMAKE_SYSTEM} MATCHES Darwin)
-    set_target_properties(${MODULE} PROPERTIES LINK_FLAGS "-undefined dynamic_lookup")
-endif(${CMAKE_SYSTEM} MATCHES Darwin)
-
-# Installation
-install(TARGETS ${MODULE} 
-  ARCHIVE DESTINATION lib
-  LIBRARY DESTINATION lib)
-install(FILES ${HDRS} DESTINATION include)
\ No newline at end of file
+message(STATUS "PWG mftmuondep charmonium enabled")
diff --git a/PWG/mftmuondep/charmonium/PWGmftmuondepCharmoniumLinkDef.h b/PWG/mftmuondep/charmonium/PWGmftmuondepCharmoniumLinkDef.h
new file mode 100644 (file)
index 0000000..ba5509f
--- /dev/null
@@ -0,0 +1,11 @@
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliAnalysisTaskDimuonBackground.cxx
+
+#endif
+
+
diff --git a/PWG/mftmuondep/charmonium/RunAnalysisTaskDimuonBackground.C b/PWG/mftmuondep/charmonium/RunAnalysisTaskDimuonBackground.C
new file mode 100755 (executable)
index 0000000..dab9875
--- /dev/null
@@ -0,0 +1,297 @@
+//====================================================================================================================================================
+
+Bool_t RunAnalysisTaskDimuonBackground(const Char_t *runType="local",       // "grid" and "local" types have been tested
+                                      const Char_t *runMode="full") {
+  
+  //  enum {kGenerated, kReconstructed};
+
+  LoadLibs();
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (mgr) delete mgr;
+  mgr = new AliAnalysisManager("AM","Analysis Manager");
+  
+  AliMultiInputEventHandler *handler = new AliMultiInputEventHandler();
+  handler->AddInputEventHandler(new AliAODInputHandler());
+
+  //-------------------- Setup of the mixed event handler ---------------
+  
+  const Int_t mixNum = 20;      // Number of times the main event is used for a mixing procedure
+  const Int_t bufferSize = 1;   // Number of events involved in each mixing procedure with the main one (typically one)
+
+  AliMixInputEventHandler *mixHandler = new AliMixInputEventHandler(bufferSize, mixNum);
+  mixHandler->SetInputHandlerForMixing(dynamic_cast<AliMultiInputEventHandler*> handler);
+  AliMixEventPool *evPool = new AliMixEventPool();
+
+  AliMixEventCutObj *zVertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, -10, 10, 1);
+
+  evPool->AddCut(zVertex);
+  mixHandler->SetEventPool(evPool);
+
+  mixHandler->SelectCollisionCandidates(AliVEvent::kAny);
+  handler->AddInputEventHandler(mixHandler);
+
+  //---------------------------------------------------------------------
+
+  mgr->SetInputEventHandler(handler);
+  
+  if (!strcmp(runType,"grid")) mgr->SetGridHandler(CreateAlienHandler(runMode));
+  
+  gROOT->LoadMacro("AliAnalysisTaskDimuonBackground.cxx++g");   
+  AliAnalysisTaskDimuonBackground *task = new AliAnalysisTaskDimuonBackground("AliAnalysisTaskDimuonBackground");
+
+  // in cm, taken from Fig. 7.4 of the ITS Upgrade TDR, in the hypothesis of ~80 tracks participating to the vtx
+  // task -> SetVtxResolutionITS(5.e-4, 5.e-4, 4.e-4);
+
+  task -> SetVertexMode(AliAnalysisTaskDimuonBackground::kReconstructed);
+
+  task -> SetMinTriggerMatch(1);
+  task -> SetSingleMuonMinPt(1.0);
+  task -> SetSingleMuonMinEta(-3.6);
+  task -> SetSingleMuonMaxEta(-2.5);
+  task -> SetSingleMuonMaxChi2(9999.);
+  
+  // create output container(s)
+  AliAnalysisDataContainer *histogramList = mgr->CreateContainer("list", TList::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  
+  // connect input and output
+  mgr->ConnectInput(task,  0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task, 1, histogramList);
+
+  // RUN ANALYSIS
+
+  TStopwatch timer;
+  timer.Start();
+  mgr->InitAnalysis();
+  mgr->PrintStatus();
+  
+  if (!strcmp(runType,"grid")) {
+    printf("Starting MFT analysis on the grid");
+    mgr->StartAnalysis("grid");
+  }
+
+  else if (!strcmp(runType,"local")) {
+    printf("Starting MFT analysis locally");
+    mgr->StartAnalysis("local", GetInputLocalData());
+  }
+
+  else AliError(Form("Specified run type %s is not supported", runType));
+  
+  timer.Stop();
+  timer.Print();
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+AliAnalysisGrid* CreateAlienHandler(const Char_t *runMode) {
+
+  // Set up the analysis plugin in case of a grid analysis
+
+  TString rootVersion = "v5-34-08-6";
+  TString alirootVersion = "vAN-20141128";
+
+  AliAnalysisAlien *analysisPlugin = new AliAnalysisAlien();
+  if (!analysisPlugin) { Printf("Error : analysisPlugin is null !!!"); return kFALSE; }
+  analysisPlugin->SetAPIVersion("V1.1x");
+  analysisPlugin->SetROOTVersion(rootVersion.Data());
+  analysisPlugin->SetAliROOTVersion(alirootVersion.Data());
+  analysisPlugin->SetExecutableCommand("aliroot -b -q");
+
+  // Overwrite all generated files, datasets and output results from a previous session
+  analysisPlugin->SetOverwriteMode();
+  // Set the run mode (can be "full", "test", "offline", "submit" or "terminate")
+  analysisPlugin->SetRunMode(runMode);  // VERY IMPORTANT 
+  
+  analysisPlugin->SetAdditionalRootLibs("CORRFW");
+  analysisPlugin->SetAdditionalRootLibs("PWGmuon");
+  
+  // Method1: Location of Data and Working dir
+  // analysisPlugin->SetGridDataDir("/alice/sim/2014/LHC14j5_new_plus/");
+  // analysisPlugin->SetDataPattern("AliAOD.root");
+  // analysisPlugin->AddRunNumber(137161);
+
+  // Method2: Use existing xml collection
+
+  analysisPlugin->AddDataFile("../../collections/137844_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/137848_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/138190_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/138192_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/138197_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/138201_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/138225_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/138275_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/138364_1.xml"); 
+  analysisPlugin->AddDataFile("../../collections/138396_1.xml");
+
+  // analysisPlugin->AddDataFile("137230_2.xml");
+  // analysisPlugin->AddDataFile("137162_2.xml");
+  // analysisPlugin->AddDataFile("137161_2.xml");
+  // analysisPlugin->AddDataFile("137231_2.xml");
+  // analysisPlugin->AddDataFile("137430_2.xml");
+  // analysisPlugin->AddDataFile("137366_2.xml");
+  // analysisPlugin->AddDataFile("137243_2.xml");
+  // analysisPlugin->AddDataFile("137236_2.xml");
+  // analysisPlugin->AddDataFile("137232_2.xml");
+  // analysisPlugin->AddDataFile("137440_2.xml");
+  // analysisPlugin->AddDataFile("137434_2.xml");
+  // analysisPlugin->AddDataFile("137432_2.xml");
+  // analysisPlugin->AddDataFile("137431_2.xml");
+  // analysisPlugin->AddDataFile("137443_2.xml");
+  // analysisPlugin->AddDataFile("137441_2.xml");
+  // analysisPlugin->AddDataFile("137544_2.xml");
+  // analysisPlugin->AddDataFile("137541_2.xml");
+  // analysisPlugin->AddDataFile("137539_2.xml");
+  // analysisPlugin->AddDataFile("137549_2.xml");
+  // analysisPlugin->AddDataFile("137608_2.xml");
+  // analysisPlugin->AddDataFile("137595_2.xml");
+  // analysisPlugin->AddDataFile("137686_2.xml");
+  // analysisPlugin->AddDataFile("137639_2.xml");
+  // analysisPlugin->AddDataFile("137638_2.xml");
+  // analysisPlugin->AddDataFile("137692_2.xml");
+  // analysisPlugin->AddDataFile("137691_2.xml");
+  // analysisPlugin->AddDataFile("137718_2.xml");
+  // analysisPlugin->AddDataFile("137704_2.xml");
+  // analysisPlugin->AddDataFile("137751_2.xml");
+  // analysisPlugin->AddDataFile("137724_2.xml");
+  // analysisPlugin->AddDataFile("137722_2.xml");
+  // analysisPlugin->AddDataFile("137752_2.xml");
+  // analysisPlugin->AddDataFile("137844_2.xml");
+  // analysisPlugin->AddDataFile("138190_2.xml");
+  // analysisPlugin->AddDataFile("137848_2.xml");
+  // analysisPlugin->AddDataFile("138197_2.xml");
+  // analysisPlugin->AddDataFile("138192_2.xml");
+  // analysisPlugin->AddDataFile("138225_2.xml");
+  // analysisPlugin->AddDataFile("138201_2.xml");
+  // analysisPlugin->AddDataFile("138275_2.xml");
+  // analysisPlugin->AddDataFile("138364_2.xml");
+  // analysisPlugin->AddDataFile("138396_2.xml");
+  
+  // analysisPlugin->AddDataFile("137161_3.xml");
+  // analysisPlugin->AddDataFile("137162_3.xml");
+  // analysisPlugin->AddDataFile("137230_3.xml");
+  // analysisPlugin->AddDataFile("137231_3.xml");
+  // analysisPlugin->AddDataFile("137232_3.xml");
+  // analysisPlugin->AddDataFile("137236_3.xml");
+  // analysisPlugin->AddDataFile("137243_3.xml");
+  // analysisPlugin->AddDataFile("137366_3.xml");
+  // analysisPlugin->AddDataFile("137430_3.xml");
+  // analysisPlugin->AddDataFile("137431_3.xml");
+  // analysisPlugin->AddDataFile("137432_3.xml");
+  // analysisPlugin->AddDataFile("137434_3.xml");
+  // analysisPlugin->AddDataFile("137440_3.xml");
+  // analysisPlugin->AddDataFile("137441_3.xml");
+  // analysisPlugin->AddDataFile("137443_3.xml");
+  // analysisPlugin->AddDataFile("137539_3.xml");
+  // analysisPlugin->AddDataFile("137541_3.xml");
+  // analysisPlugin->AddDataFile("137544_3.xml");
+  // analysisPlugin->AddDataFile("137549_3.xml");
+  // analysisPlugin->AddDataFile("137595_3.xml");
+  // analysisPlugin->AddDataFile("137608_3.xml");
+  // analysisPlugin->AddDataFile("137638_3.xml");
+  // analysisPlugin->AddDataFile("137639_3.xml");
+  // analysisPlugin->AddDataFile("137686_3.xml");
+  // analysisPlugin->AddDataFile("137691_3.xml");
+  // analysisPlugin->AddDataFile("137692_3.xml");
+  // analysisPlugin->AddDataFile("137704_3.xml");
+  // analysisPlugin->AddDataFile("137718_3.xml");
+  // analysisPlugin->AddDataFile("137722_3.xml");
+  // analysisPlugin->AddDataFile("137724_3.xml");
+  // analysisPlugin->AddDataFile("137751_3.xml");
+  // analysisPlugin->AddDataFile("137752_3.xml");
+  // analysisPlugin->AddDataFile("137844_3.xml");
+  // analysisPlugin->AddDataFile("137848_3.xml");
+  // analysisPlugin->AddDataFile("138190_3.xml");
+  // analysisPlugin->AddDataFile("138192_3.xml");
+  // analysisPlugin->AddDataFile("138197_3.xml");
+  // analysisPlugin->AddDataFile("138201_3.xml");
+  // analysisPlugin->AddDataFile("138225_3.xml");
+  // analysisPlugin->AddDataFile("138275_3.xml");
+  // analysisPlugin->AddDataFile("138364_3.xml");
+  // analysisPlugin->AddDataFile("138396_3.xml");
+  
+  analysisPlugin->SetGridWorkingDir("MFT/analysis/LHC14cj5_PbPb_0_5/charmonia/background");
+  
+  // Declare alien output directory. Relative to working directory.
+  analysisPlugin->SetGridOutputDir("output"); // In this case will be $HOME/work/output
+  // Declare the analysis source files names separated by blancs. To be compiled runtime using ACLiC on the worker nodes.
+  analysisPlugin->SetAnalysisSource("AliAnalysisTaskDimuonBackground.cxx");
+  // Declare all libraries (other than the default ones for the framework. These will be
+  // loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here.
+  analysisPlugin->SetAdditionalLibs("libCORRFW.so libPWGmuon.so libEventMixing.so AliAnalysisTaskDimuonBackground.h AliAnalysisTaskDimuonBackground.cxx");
+  
+  analysisPlugin->AddIncludePath("-I.");
+
+  analysisPlugin->SetOutputToRunNo();
+  
+  // Optionally set a name for the generated analysis macro (default MyAnalysis.C)
+  analysisPlugin->SetAnalysisMacro("MFTAnalysis.C");
+
+  // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore)
+  analysisPlugin->SetSplitMaxInputFileNumber(100);
+  // Number of runs per master job
+  analysisPlugin->SetNrunsPerMaster(1);
+  
+  // Optionally modify the executable name (default analysis.sh)
+  analysisPlugin->SetExecutable("MFTAnalysis.sh");
+
+  // Optionally set number of failed jobs that will trigger killing waiting sub-jobs.
+  analysisPlugin->SetMaxInitFailed(5);
+  // Optionally resubmit threshold.
+  analysisPlugin->SetMasterResubmitThreshold(90);
+  // Optionally set time to live (default 30000 sec)
+  analysisPlugin->SetTTL(30000);
+  // Optionally set input format (default xml-single)
+  analysisPlugin->SetInputFormat("xml-single");
+  // Optionally modify the name of the generated JDL (default analysis.jdl)
+  analysisPlugin->SetJDLName("MFTAnalysis.jdl");
+  // Optionally modify job price (default 1)
+  analysisPlugin->SetPrice(1);      
+  // Optionally modify split mode (default 'se')    
+  analysisPlugin->SetSplitMode("se");
+  
+  analysisPlugin->SetNtestFiles(5);
+  //   analysisPlugin->SetMergeViaJDL(1);
+  analysisPlugin->SetOverwriteMode(kTRUE);
+  
+  return analysisPlugin;
+
+}
+
+//====================================================================================================================================================
+
+TChain* GetInputLocalData() {
+
+  // Set up the chain of input events in case of a local analysis
+
+  TChain *chain = new TChain("aodTree");
+  chain->Add("./test_03Nov2014_zVtx_Rndm/AliAOD.Muons.root");
+
+  return chain;
+
+}
+
+//====================================================================================================================================================
+
+void LoadLibs() {
+
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include ");
+  
+  gSystem->Load("libTree.so")          ;
+  gSystem->Load("libGeom.so")          ;
+  gSystem->Load("libVMC.so")           ;
+  gSystem->Load("libMinuit.so")        ;
+  gSystem->Load("libPhysics.so")       ;
+  gSystem->Load("libSTEERBase.so")     ;
+  gSystem->Load("libESD.so")           ;
+  gSystem->Load("libAOD.so")           ;
+  gSystem->Load("libANALYSIS.so")      ;
+  gSystem->Load("libOADB.so")          ;
+  gSystem->Load("libANALYSISalice.so") ;
+  gSystem->Load("libCORRFW.so")        ;
+  gSystem->Load("libPWGmuon.so")       ;
+
+}
+
+//====================================================================================================================================================