MFT Analysis support classes added
authorauras <antonio.uras@cern.ch>
Tue, 25 Feb 2014 18:16:07 +0000 (19:16 +0100)
committerauras <antonio.uras@cern.ch>
Tue, 25 Feb 2014 18:16:07 +0000 (19:16 +0100)
MFT/AliAnalysisTaskMFTExample.cxx [new file with mode: 0755]
MFT/AliAnalysisTaskMFTExample.h [new file with mode: 0755]
MFT/AliMFTSupport.cxx [new file with mode: 0644]
MFT/AliMFTSupport.h [new file with mode: 0644]
MFT/RunAnalysisTaskMFTExample.C [new file with mode: 0755]

diff --git a/MFT/AliAnalysisTaskMFTExample.cxx b/MFT/AliAnalysisTaskMFTExample.cxx
new file mode 100755 (executable)
index 0000000..4768137
--- /dev/null
@@ -0,0 +1,251 @@
+#include "TH1D.h"
+#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliAnalysisTaskMFTExample.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 "AliMFTSupport.h"
+#include "TRandom.h"
+
+ClassImp(AliAnalysisTaskMFTExample)
+
+//====================================================================================================================================================
+
+AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample() : 
+  AliAnalysisTaskSE(),
+  fVertexMode(0),
+  fHistPtSingleMuons(0),
+  fHistPtSingleMuonsFromJpsi(0),
+  fHistPtDimuonsOS(0),
+  fHistMassDimuonsOS(0),
+  fHistPtDimuonsJpsi(0),
+  fHistMassDimuonsJpsi(0),
+  fHistResidualXVtxJpsi(0),
+  fHistResidualYVtxJpsi(0),
+  fHistResidualZVtxJpsi(0) {
+  
+  // Default constructor
+
+  for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
+  fVtxResolutionITS[0] = 5.e-4;
+  fVtxResolutionITS[1] = 5.e-4;
+  fVtxResolutionITS[2] = 4.e-4;
+
+}
+
+//====================================================================================================================================================
+
+AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample(const Char_t *name) : 
+  AliAnalysisTaskSE(name),
+  fVertexMode(0),
+  fHistPtSingleMuons(0),
+  fHistPtSingleMuonsFromJpsi(0),
+  fHistPtDimuonsOS(0),
+  fHistMassDimuonsOS(0),
+  fHistPtDimuonsJpsi(0),
+  fHistMassDimuonsJpsi(0),
+  fHistResidualXVtxJpsi(0),
+  fHistResidualYVtxJpsi(0),
+  fHistResidualZVtxJpsi(0) {
+
+  // Constructor
+
+  for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
+  fVtxResolutionITS[0] = 5.e-4;
+  fVtxResolutionITS[1] = 5.e-4;
+  fVtxResolutionITS[2] = 4.e-4;
+
+  // Define input and output slots here
+  DefineOutput(1, TH1D::Class());
+  DefineOutput(2, TH1D::Class());
+  DefineOutput(3, TH1D::Class());
+  DefineOutput(4, TH1D::Class());
+  DefineOutput(5, TH1D::Class());
+  DefineOutput(6, TH1D::Class());
+  DefineOutput(7, TH1D::Class());
+  DefineOutput(8, TH1D::Class());
+  DefineOutput(9, TH1D::Class());
+
+}
+
+//====================================================================================================================================================
+
+void AliAnalysisTaskMFTExample::UserCreateOutputObjects() {
+
+  // Called once
+
+  fHistPtSingleMuons = new TH1D("fHistPtSingleMuons","p_{T} of single muons (All)", 100, 0, 10);
+  fHistPtSingleMuons -> SetXTitle("p_{T}  [GeV/c]");
+  fHistPtSingleMuons -> Sumw2();
+
+  fHistPtSingleMuonsFromJpsi = new TH1D("fHistPtSingleMuonsFromJpsi", "p_{T} of single muons (from J/#psi)", 100, 0, 10);
+  fHistPtSingleMuonsFromJpsi -> SetXTitle("p_{T}  [GeV/c]");
+  fHistPtSingleMuonsFromJpsi -> Sumw2();
+
+  fHistMassDimuonsOS = new TH1D("fHistMassDimuonsOS", "Mass of OS dimuons", 500, 0, 10);
+  fHistMassDimuonsOS -> SetXTitle("Mass  [GeV/c^{2}]");
+  fHistMassDimuonsOS -> Sumw2();
+
+  fHistMassDimuonsJpsi = new TH1D("fHistMassDimuonsJpsi", "Mass of J/#psi dimuons", 500, 0, 10);
+  fHistMassDimuonsJpsi -> SetXTitle("Mass  [GeV/c^{2}]");
+  fHistMassDimuonsJpsi -> Sumw2();
+
+  fHistPtDimuonsOS = new TH1D("fHistPtDimuonsOS", "p_{T} of OS dimuons", 100, 0, 10);
+  fHistPtDimuonsOS -> SetXTitle("p_{T}  [GeV/c]");
+  fHistPtDimuonsOS -> Sumw2();
+
+  fHistPtDimuonsJpsi = new TH1D("fHistPtDimuonsJpsi", "p_{T} of J/#psi dimuons", 100, 0, 10);
+  fHistPtDimuonsJpsi -> SetXTitle("p_{T}  [GeV/c]");
+  fHistPtDimuonsJpsi -> Sumw2();
+
+  fHistResidualXVtxJpsi = new TH1D("fHistResidualXVtxJpsi", "J/#psi vertex residuals along x", 100, -100, 100);
+  fHistResidualXVtxJpsi -> SetXTitle("Residual  [#mum]");
+  fHistResidualXVtxJpsi -> Sumw2();
+
+  fHistResidualYVtxJpsi = new TH1D("fHistResidualYVtxJpsi", "J/#psi vertex residuals along y", 100, -100, 100);
+  fHistResidualYVtxJpsi -> SetXTitle("Residual  [#mum]");
+  fHistResidualYVtxJpsi -> Sumw2();
+
+  fHistResidualZVtxJpsi = new TH1D("fHistResidualZVtxJpsi", "J/#psi vertex residuals along z", 100, -1000, 1000);
+  fHistResidualZVtxJpsi -> SetXTitle("Residual  [#mum]");
+  fHistResidualZVtxJpsi -> Sumw2();
+
+  PostData(1, fHistPtSingleMuons);
+  PostData(2, fHistPtSingleMuonsFromJpsi);
+  PostData(3, fHistMassDimuonsOS);
+  PostData(4, fHistMassDimuonsJpsi);
+  PostData(5, fHistPtDimuonsOS);
+  PostData(6, fHistPtDimuonsJpsi);
+  PostData(7, fHistResidualXVtxJpsi);
+  PostData(8, fHistResidualYVtxJpsi);
+  PostData(9, fHistResidualZVtxJpsi);
+  
+}
+
+//====================================================================================================================================================
+
+void AliAnalysisTaskMFTExample::UserExec(Option_t *) {
+
+  // Main loop
+  // Called for each event     
+
+  AliAODEvent *aodEv = dynamic_cast<AliAODEvent*>(InputEvent());  
+  if (!aodEv) return;
+
+  aodEv->GetHeader()->InitMagneticField();
+  AliMUONTrackExtrap::SetField();
+
+  TClonesArray *stackMC = (TClonesArray*) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
+
+  // Getting primary vertex, either from the generation or from the reconstruction -------------------
+
+  if (fVertexMode == kGenerated) {
+    AliAODMCHeader *mcHeader = (AliAODMCHeader*) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+    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);
+  }
+
+  // -------------------------------------------------------------------------------------------------
+
+  AliAODTrack *recMuon1=0, *recMuon2=0;
+  AliAODMCParticle *mcMuon1=0, *mcMuon2=0;
+  AliAODMCParticle *mcMother1=0;
+  Bool_t isMuon1FromJpsi=0;
+
+  // Loop over MUON+MFT muons
+
+  //--------------------------------------------------------------------------------------------------
+
+  for (Int_t iTrack=0; iTrack<aodEv->GetNumberOfTracks(); iTrack++) { 
+
+    if (!(aodEv->GetTrack(iTrack)->IsMuonGlobalTrack())) continue;
+
+    recMuon1 = aodEv->GetTrack(iTrack);
+    isMuon1FromJpsi = kFALSE;
+    
+    // pt all muons
+    fHistPtSingleMuons -> Fill(recMuon1->Pt());
+
+    // pt muons from J/psi
+    if (recMuon1->GetLabel() >= 0) {
+      mcMuon1 = (AliAODMCParticle*) stackMC->At(recMuon1->GetLabel());
+      if (mcMuon1) {
+       if (mcMuon1->GetMother()>=0) {
+         mcMother1 = (AliAODMCParticle*) stackMC->At(mcMuon1->GetMother());
+         if (mcMother1->PdgCode()==443) {
+           isMuon1FromJpsi = kTRUE;
+           fHistPtSingleMuonsFromJpsi -> Fill(recMuon1->Pt());
+         }
+       }
+      }
+    }
+
+    for (Int_t jTrack=0; jTrack<iTrack; jTrack++) { 
+
+      if (!(aodEv->GetTrack(jTrack)->IsMuonGlobalTrack())) continue;
+      
+      recMuon2 = aodEv->GetTrack(jTrack);
+      
+      AliAODDimuon *dimuon = new AliAODDimuon(recMuon1, recMuon2);
+      if (dimuon->Charge()) continue;
+
+      // pt and mass all OS dimuons
+      fHistPtDimuonsOS   -> Fill(dimuon->Pt());
+      fHistMassDimuonsOS -> Fill(dimuon->Mass());
+
+      // pt and mass J/psi dimuons
+      if (!isMuon1FromJpsi) continue;
+      if (recMuon2->GetLabel() >= 0) {
+       mcMuon2 = (AliAODMCParticle*) stackMC->At(recMuon2->GetLabel());
+       if (mcMuon2) {
+         if (mcMuon2->GetMother() == mcMuon1->GetMother()) {
+           Double_t vertex[3]={0};
+           TLorentzVector kinem(0,0,0,0);
+           if (!AliMFTSupport::RefitAODDimuonWithCommonVertex(dimuon, vertex, kinem)) continue;
+           fHistPtDimuonsJpsi    -> Fill(kinem.Pt());
+           fHistMassDimuonsJpsi  -> Fill(kinem.M());
+           fHistResidualXVtxJpsi -> Fill(1.e4*(vertex[0] - fPrimaryVertex[0]));
+           fHistResidualYVtxJpsi -> Fill(1.e4*(vertex[1] - fPrimaryVertex[1]));
+           fHistResidualZVtxJpsi -> Fill(1.e4*(vertex[2] - fPrimaryVertex[2]));
+         }
+       }
+      }
+   
+    }   // end of loop on 2nd muon
+
+  }   // end of loop on 1st muon
+  //--------------------------------------------------------------------------------------------------
+   
+  PostData(1, fHistPtSingleMuons);
+  PostData(2, fHistPtSingleMuonsFromJpsi);
+  PostData(3, fHistMassDimuonsOS);
+  PostData(4, fHistMassDimuonsJpsi);
+  PostData(5, fHistPtDimuonsOS);
+  PostData(6, fHistPtDimuonsJpsi);
+  PostData(7, fHistResidualXVtxJpsi);
+  PostData(8, fHistResidualYVtxJpsi);
+  PostData(9, fHistResidualZVtxJpsi);
+
+}
+
+//====================================================================================================================================================
+
+void AliAnalysisTaskMFTExample::Terminate(Option_t *) {
+
+  // Draw result to the screen
+  // Called once at the end of the query
+  
+}
+
+//====================================================================================================================================================
diff --git a/MFT/AliAnalysisTaskMFTExample.h b/MFT/AliAnalysisTaskMFTExample.h
new file mode 100755 (executable)
index 0000000..f49c181
--- /dev/null
@@ -0,0 +1,58 @@
+#ifndef AliAnalysisTaskMFTExample_H
+#define AliAnalysisTaskMFTExample_H
+
+#include "AliAnalysisTaskSE.h"
+#include "TH1D.h"
+
+//====================================================================================================================================================
+
+class  AliAnalysisTaskMFTExample : public AliAnalysisTaskSE {
+
+public:
+
+  enum {kGenerated, kReconstructed};
+  AliAnalysisTaskMFTExample();
+  AliAnalysisTaskMFTExample(const char *name);
+
+  virtual ~AliAnalysisTaskMFTExample() {
+    delete fHistPtSingleMuons;
+    delete fHistPtSingleMuonsFromJpsi;
+    delete fHistPtDimuonsOS;
+    delete fHistMassDimuonsOS;
+    delete fHistPtDimuonsJpsi;
+    delete fHistMassDimuonsJpsi;
+    delete fHistResidualXVtxJpsi;
+    delete fHistResidualYVtxJpsi;
+    delete fHistResidualZVtxJpsi;
+  }
+  
+  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;
+  }
+
+  virtual void UserCreateOutputObjects();
+  virtual void UserExec(Option_t *option);
+  virtual void Terminate(Option_t *);
+
+private:
+
+  Double_t fPrimaryVertex[3], fVtxResolutionITS[3];
+  Int_t fVertexMode;
+
+  TH1D *fHistPtSingleMuons, *fHistPtSingleMuonsFromJpsi;
+  TH1D *fHistPtDimuonsOS, *fHistMassDimuonsOS;
+  TH1D *fHistPtDimuonsJpsi, *fHistMassDimuonsJpsi;
+
+  TH1D *fHistResidualXVtxJpsi, *fHistResidualYVtxJpsi, *fHistResidualZVtxJpsi;
+
+  ClassDef(AliAnalysisTaskMFTExample, 1) // example of analysis
+
+};
+
+//====================================================================================================================================================
+
+#endif
diff --git a/MFT/AliMFTSupport.cxx b/MFT/AliMFTSupport.cxx
new file mode 100644 (file)
index 0000000..b24c264
--- /dev/null
@@ -0,0 +1,185 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * 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.                  *
+ **************************************************************************/
+
+//====================================================================================================================================================
+//
+//      Support class for various common operation on MFT objects
+//
+//      Contact author: antonio.uras@cern.ch
+//
+//====================================================================================================================================================
+
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
+#include "AliAODTrack.h"
+#include "AliAODDimuon.h"
+#include "TLorentzVector.h"
+#include "AliMFTConstants.h"
+#include "TDatabasePDG.h"
+#include "TMath.h"
+#include "AliLog.h"
+
+#include "AliMFTSupport.h"
+
+ClassImp(AliMFTSupport)
+
+//====================================================================================================================================================
+
+Bool_t AliMFTSupport::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]) {
+
+  if (!(muon->Pz()!=0)) return kFALSE;
+
+  AliMUONTrackParam *param = new AliMUONTrackParam();
+
+  param -> SetNonBendingCoor(muon->XAtDCA());
+  param -> SetBendingCoor(muon->YAtDCA());
+  param -> SetZ(0.);
+  param -> SetNonBendingSlope(muon->Px()/muon->Pz());
+  param -> SetBendingSlope(muon->Py()/muon->Pz());
+  param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
+
+  AliMUONTrackExtrap::ExtrapToZ(param, z);
+  xy[0] = param->GetNonBendingCoor();
+  xy[1] = param->GetBendingCoor();
+
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTSupport::RefitAODDimuonWithCommonVertex(AliAODDimuon *dimuon, Double_t *vertex, TLorentzVector &kinem) {
+
+  Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0;
+
+  AliAODTrack *muon0 = dimuon->GetMu(0);
+  AliAODTrack *muon1 = dimuon->GetMu(1);
+
+  if (!(muon0->Pz()!=0 && muon1->Pz()!=0)) return kFALSE;
+
+  AliMUONTrackParam *param0 = new AliMUONTrackParam();
+  AliMUONTrackParam *param1 = new AliMUONTrackParam();
+
+  param0 -> SetNonBendingCoor(muon0->XAtDCA());
+  param1 -> SetNonBendingCoor(muon1->XAtDCA());
+
+  param0 -> SetBendingCoor(muon0->YAtDCA());
+  param1 -> SetBendingCoor(muon1->YAtDCA());
+
+  param0 -> SetZ(0.);
+  param1 -> SetZ(0.);
+  param0 -> SetNonBendingSlope(muon0->Px()/muon0->Pz());
+  param1 -> SetNonBendingSlope(muon1->Px()/muon1->Pz());
+
+  param0 -> SetBendingSlope(muon0->Py()/muon0->Pz());
+  param1 -> SetBendingSlope(muon1->Py()/muon1->Pz());
+
+  param0 -> SetInverseBendingMomentum( muon0->Charge() * (1./muon0->Pz()) / (TMath::Sqrt(1+TMath::Power(muon0->Py()/muon0->Pz(),2))) );
+  param1 -> SetInverseBendingMomentum( muon1->Charge() * (1./muon1->Pz()) / (TMath::Sqrt(1+TMath::Power(muon1->Py()/muon1->Pz(),2))) );
+
+  Double_t step = 1.;  // in cm
+  Double_t startPoint = 0.;
+
+  Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
+  
+  for (Int_t i=0; i<3; i++) {
+    AliMUONTrackExtrap::ExtrapToZ(param0, z[i]);
+    AliMUONTrackExtrap::ExtrapToZ(param1, z[i]);
+    Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
+    Double_t dY = param0->GetBendingCoor()    - param1->GetBendingCoor();
+    r[i] = TMath::Sqrt(dX*dX + dY*dY);
+  }
+  
+  Double_t researchDirection=0.;
+  
+  if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1.;   // towards z positive
+  else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1.;   // towards z negative
+  else if (r[0]<r[1] && r[1]>r[2]) { 
+    printf("E-AliMFTSupport::RefitAODDimuonWithCommonVertex: Point of closest approach cannot be found for dimuon (no minima)\n");
+    return kFALSE;
+  }
+  
+  while (TMath::Abs(researchDirection)>0.5) {
+    if (researchDirection>0.) {
+      z[0] = z[1];
+      z[1] = z[2];
+      z[2] = z[1]+researchDirection*step;
+    }
+    else {
+      z[2] = z[1];
+      z[1] = z[0];
+      z[0] = z[1]+researchDirection*step;
+    }
+    if (TMath::Abs(z[0])>900.) {
+      printf("E-AliMFTSupport::RefitAODDimuonWithCommonVertex: Point of closest approach cannot be found for dimuon (no minima)\n");
+      return kFALSE;
+    }
+    for (Int_t i=0; i<3; i++) {
+      AliMUONTrackExtrap::ExtrapToZ(param0, z[i]);
+      AliMUONTrackExtrap::ExtrapToZ(param1, z[i]);
+      Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
+      Double_t dY = param0->GetBendingCoor()    - param1->GetBendingCoor();
+      r[i] = TMath::Sqrt(dX*dX + dY*dY);
+    }
+    researchDirection=0.;
+    if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1.;   // towards z positive
+    else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1.;   // towards z negative
+  }
+    
+  step *= 0.5;
+  while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
+    z[0] = z[1]-step;
+    z[2] = z[1]+step;
+    for (Int_t i=0; i<3; i++) {
+      AliMUONTrackExtrap::ExtrapToZ(param0, z[i]);
+      AliMUONTrackExtrap::ExtrapToZ(param1, z[i]);
+      Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
+      Double_t dY = param0->GetBendingCoor()    - param1->GetBendingCoor();
+      r[i] = TMath::Sqrt(dX*dX + dY*dY);
+    }
+    if      (r[0]<r[1]) z[1] = z[0];
+    else if (r[2]<r[1]) z[1] = z[2];
+    else step *= 0.5;
+  }
+  
+  fZPointOfClosestApproach = z[1];
+  AliMUONTrackExtrap::ExtrapToZ(param0, fZPointOfClosestApproach);
+  AliMUONTrackExtrap::ExtrapToZ(param1, fZPointOfClosestApproach);  
+  fXPointOfClosestApproach = 0.5*(param0->GetNonBendingCoor() + param1->GetNonBendingCoor());
+  fYPointOfClosestApproach = 0.5*(param0->GetBendingCoor()    + param1->GetBendingCoor());
+  
+  //-------------------------------------------------------------------------------------------------------------------------------------------------
+
+  vertex[0] = fXPointOfClosestApproach;
+  vertex[1] = fYPointOfClosestApproach;
+  vertex[2] = fZPointOfClosestApproach;
+
+  Double_t pTot[3] = {0};
+  pTot[0] = param0->Px() + param1->Px();
+  pTot[1] = param0->Py() + param1->Py();
+  pTot[2] = param0->Pz() + param1->Pz();
+
+  Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
+
+  Double_t energy = TMath::Sqrt(massMu*massMu + pTot[0]*pTot[0] + pTot[1]*pTot[1] + pTot[2]*pTot[2]); 
+
+  kinem.SetPxPyPzE(pTot[0], pTot[1], pTot[2], energy);
+
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
diff --git a/MFT/AliMFTSupport.h b/MFT/AliMFTSupport.h
new file mode 100644 (file)
index 0000000..980fa2c
--- /dev/null
@@ -0,0 +1,49 @@
+#ifndef AliMFTSupport_H
+#define AliMFTSupport_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//====================================================================================================================================================
+//
+//      Support class for various common operation on MFT objects
+//
+//      Contact author: antonio.uras@cern.ch
+//
+//====================================================================================================================================================
+
+#include "TObject.h"
+#include "AliAODTrack.h"
+#include "AliAODDimuon.h"
+#include "TLorentzVector.h"
+#include "AliMFTConstants.h"
+#include "TDatabasePDG.h"
+#include "TMath.h"
+#include "AliLog.h"
+
+//====================================================================================================================================================
+
+class AliMFTSupport : public TObject {
+
+public:
+
+  AliMFTSupport() : TObject() {;}
+  virtual ~AliMFTSupport() {;}
+  
+  static Bool_t ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2]);
+  static Bool_t RefitAODDimuonWithCommonVertex(AliAODDimuon *dimuon, Double_t *vertex, TLorentzVector &kinem);
+
+  static Bool_t PlaneExists(AliAODTrack *muon, Int_t iPlane) { return muon->GetMFTClusterPattern() & (1<<iPlane); }
+
+  static Bool_t IsWrongCluster(AliAODTrack *muon, Int_t iPlane) { 
+    if (!PlaneExists(muon, iPlane)) return kFALSE;
+    else return !(muon->GetMFTClusterPattern() & (1<<(iPlane+AliMFTConstants::fNMaxPlanes)));
+  }
+
+  ClassDef(AliMFTSupport,1)
+    
+};
+
+//====================================================================================================================================================
+
+#endif
diff --git a/MFT/RunAnalysisTaskMFTExample.C b/MFT/RunAnalysisTaskMFTExample.C
new file mode 100755 (executable)
index 0000000..86d6d0a
--- /dev/null
@@ -0,0 +1,185 @@
+Bool_t RunAnalysisTaskMFTExample(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");
+  
+  AliAODInputHandler* inputHandler = new AliAODInputHandler();
+  mgr->SetInputEventHandler(inputHandler);
+  
+  if (!strcmp(runType,"grid")) mgr->SetGridHandler(CreateAlienHandler());
+  
+  gROOT->LoadMacro("AliAnalysisTaskMFTExample.cxx++g");   
+  AliAnalysisTaskMFTExample *task = new AliAnalysisTaskMFTExample("AliAnalysisTaskMFTExample");
+
+  // 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(kGenerated);
+
+  mgr->AddTask(task);
+  
+  // create output container(s)
+  AliAnalysisDataContainer *output1 = mgr->CreateContainer("output1",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  AliAnalysisDataContainer *output2 = mgr->CreateContainer("output2",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  AliAnalysisDataContainer *output3 = mgr->CreateContainer("output3",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  AliAnalysisDataContainer *output4 = mgr->CreateContainer("output4",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  AliAnalysisDataContainer *output5 = mgr->CreateContainer("output5",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  AliAnalysisDataContainer *output6 = mgr->CreateContainer("output6",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  AliAnalysisDataContainer *output7 = mgr->CreateContainer("output7",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  AliAnalysisDataContainer *output8 = mgr->CreateContainer("output8",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  AliAnalysisDataContainer *output9 = mgr->CreateContainer("output9",   TH1D::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root");
+  
+  // connect input and output
+  mgr->ConnectInput(task,  0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task, 1, output1);
+  mgr->ConnectOutput(task, 2, output2);
+  mgr->ConnectOutput(task, 3, output3);
+  mgr->ConnectOutput(task, 4, output4);
+  mgr->ConnectOutput(task, 5, output5);
+  mgr->ConnectOutput(task, 6, output6);
+  mgr->ConnectOutput(task, 7, output7);
+  mgr->ConnectOutput(task, 8, output8);
+  mgr->ConnectOutput(task, 9, output9);
+  
+  // 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() {
+
+  // Set up the analysis plugin in case of a grid analysis
+
+  TString rootVersion = "v5-34-02-1";
+  TString alirootVersion = "v5-04-33-AN";
+
+  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");
+  
+  // Location of Data and Working dir
+  analysisPlugin->SetGridDataDir("");
+  analysisPlugin->SetDataPattern("");
+  analysisPlugin->SetRunPrefix("");
+  analysisPlugin->SetGridWorkingDir("");
+  
+  // 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("AliAnalysisTaskMFTExample.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 AliAnalysisTaskMFTExample.h AliAnalysisTaskMFTExample.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("./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")       ;
+
+}
+
+//====================================================================================================================================================