]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Copy contructor and operator = (Christian)
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Jan 2004 15:25:44 +0000 (15:25 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Jan 2004 15:25:44 +0000 (15:25 +0000)
MUON/AliMUONData.cxx
MUON/AliMUONData.h
MUON/AliMUONDigit.cxx
MUON/AliMUONDigit.h
MUON/AliMUONTrack.cxx
MUON/AliMUONTrackParam.cxx
MUON/AliMUONTrackParam.h
MUON/MUONmassPlot_NewIO.C [new file with mode: 0644]

index 35dc26130be1db48c9c8b8855b1d7bf242392360..7bc2dd806dd942a698b76e3eb9dcd940deedf30d 100644 (file)
@@ -139,6 +139,15 @@ void AliMUONData::AddDigit(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digit
   new(ldigits[fNdigits[id]++]) AliMUONDigit(tracks,charges,digits);
 }
 //_____________________________________________________________________________
+void AliMUONData::AddDigit(Int_t id, const AliMUONDigit& digit)
+{
+  //
+  // Add a MUON digit to the list of Digits of the detection plane id
+  //
+  TClonesArray &ldigits = * Digits(id) ; 
+  new(ldigits[fNdigits[id]++]) AliMUONDigit(digit);
+}
+//_____________________________________________________________________________
 void AliMUONData::AddGlobalTrigger(Int_t *singlePlus, Int_t *singleMinus,
                                   Int_t *singleUndef,
                                   Int_t *pairUnlike, Int_t *pairLike)
@@ -149,6 +158,13 @@ void AliMUONData::AddGlobalTrigger(Int_t *singlePlus, Int_t *singleMinus,
     AliMUONGlobalTrigger(singlePlus, singleMinus,  singleUndef, pairUnlike, pairLike);
 }
 //_____________________________________________________________________________
+void AliMUONData::AddGlobalTrigger(const AliMUONGlobalTrigger& trigger )
+{
+  // add a MUON Global Trigger to the list (only one GlobalTrigger per event !)
+  TClonesArray &globalTrigger = *fGlobalTrigger;
+  new(globalTrigger[fNglobaltrigger++]) AliMUONGlobalTrigger(trigger);
+}
+//_____________________________________________________________________________
 void AliMUONData::AddHit(Int_t fIshunt, Int_t track, Int_t iChamber, 
                         Int_t idpart, Float_t X, Float_t Y, Float_t Z, 
                         Float_t tof, Float_t momentum, Float_t theta, 
@@ -177,12 +193,25 @@ void AliMUONData::AddHit(Int_t fIshunt, Int_t track, Int_t iChamber,
                                  Xref,Yref,Zref);
 }
 //____________________________________________________________________________
+void AliMUONData::AddHit(const AliMUONHit& hit)
+{
+  TClonesArray &lhits = *fHits;
+  new(lhits[fNhits++]) AliMUONHit(hit);
+}
+//____________________________________________________________________________
 void AliMUONData::AddLocalTrigger(Int_t *localtr)
 {
   // add a MUON Local Trigger to the list
   TClonesArray &localTrigger = *fLocalTrigger;
   new(localTrigger[fNlocaltrigger++]) AliMUONLocalTrigger(localtr);
 }
+//____________________________________________________________________________
+void AliMUONData::AddLocalTrigger(const  AliMUONLocalTrigger& trigger)
+{
+  // add a MUON Local Trigger to the list
+  TClonesArray &localTrigger = *fLocalTrigger;
+  new(localTrigger[fNlocaltrigger++]) AliMUONLocalTrigger(trigger);
+}
 //_____________________________________________________________________________
 void AliMUONData::AddRawCluster(Int_t id, const AliMUONRawCluster& c)
 {
@@ -200,6 +229,7 @@ void AliMUONData::AddRecTrack(const AliMUONTrack& track)
   //
   TClonesArray &lrectracks = *fRecTracks;
   new(lrectracks[fNrectracks++]) AliMUONTrack(track);
+  //  printf("TTTTTT %d ,\n",((AliMUONTrack*)fRecTracks->At(fNrectracks-1))->GetNTrackHits());
 }
 //____________________________________________________________________________
 TClonesArray*  AliMUONData::Digits(Int_t DetectionPlane) 
@@ -253,6 +283,7 @@ void AliMUONData::Fill(Option_t* option)
   const char *cRC  = strstr(option,"RC");  // RawCluster branches in TreeR
   const char *cGLT = strstr(option,"GLT"); // Global and Local Trigger branches in TreeR
   const char *cRT  = strstr(option,"RT");  // Reconstructed Track in TreeT
+
   //const char *cRP  = strstr(option,"RP");  // Reconstructed Particle in TreeP
   
   char branchname[30];
@@ -303,6 +334,10 @@ void AliMUONData::Fill(Option_t* option)
     sprintf(branchname,"%sTrack",GetName());  
     TreeT()->Fill();
   }
+//   if ( TreeT() && cRTT ) {
+//     sprintf(branchname,"%sTrackTrig",GetName());  
+//     TreeT()->Fill();
+//   }
 }
 //_____________________________________________________________________________
 void AliMUONData::MakeBranch(Option_t* option)
@@ -318,7 +353,7 @@ void AliMUONData::MakeBranch(Option_t* option)
   const char *cRC  = strstr(option,"RC");  // RawCluster branches in TreeR
   const char *cGLT = strstr(option,"GLT"); // Global and Local Trigger branches in TreeR
   const char *cRT  = strstr(option,"RT");  // Reconstructed Track in TreeT
-  const char *cRP  = strstr(option,"RP");  // Reconstructed Particle in TreeP
+  //const char *cRP  = strstr(option,"RP");  // Reconstructed Particle in TreeP
   
   TBranch * branch = 0x0;
   
@@ -446,10 +481,6 @@ void AliMUONData::MakeBranch(Option_t* option)
     branch = TreeT()->Branch(branchname,&fRecTracks,kBufferSize);
     Info("MakeBranch","Making Branch %s for tracks \n",branchname);
   }  
-
-  if (TreeP() && cRP ) {
-    Info("MakeBranch","Making Branch for TreeP is not yet ready. \n");
-  }
 }
 //____________________________________________________________________________
 TClonesArray*  AliMUONData::RawClusters(Int_t DetectionPlane)
@@ -606,9 +637,10 @@ void AliMUONData::SetTreeAddress(Option_t* option)
   }
 
   if ( TreeT() ) {
-    if (fRecTracks == 0x0 && cRT )  {
+    if (fRecTracks == 0x0 && cRT)  {
       fRecTracks  = new TClonesArray("AliMUONTrack",100);
     }
+
   }
   if ( TreeT() && fRecTracks && cRT ) {
     sprintf(branchname,"%sTrack",GetName());  
@@ -616,5 +648,6 @@ void AliMUONData::SetTreeAddress(Option_t* option)
     if (branch) branch->SetAddress(&fRecTracks);
     else Warning("SetTreeAddress","(%s) Failed for Tracks. Can not find branch in tree.",GetName());
   }
+
 }
 //_____________________________________________________________________________
index dd747076cdb7179bccf0f34be0f10ae9bc03fb43..b41863cd6ebb668449c75dcd6ad9ec74892d7f0e 100644 (file)
@@ -22,6 +22,10 @@ class TTree;
 class AliMUONConstants;
 class AliMUONRawCluster;
 class AliMUONTrack;
+class AliMUONDigit;
+class AliMUONHit;
+class AliMUONLocalTrigger;
+class AliMUONGlobalTrigger;
 
 //__________________________________________________________________
 /////////////////////////////////////////////////////////////////////
@@ -38,6 +42,7 @@ class AliMUONData : public TNamed {
     virtual ~AliMUONData();  
     virtual void   AddDigit(Int_t id, Int_t* tracks, Int_t* charges,
                             Int_t* digits); 
+    virtual void   AddDigit(Int_t, const AliMUONDigit& ); // use copy constructor
     virtual void   AddHit(Int_t fIshunt, Int_t track, Int_t iChamber, 
                          Int_t idpart, Float_t X, Float_t Y, Float_t Z, 
                          Float_t tof, Float_t momentum, Float_t theta, 
@@ -47,10 +52,16 @@ class AliMUONData : public TNamed {
                          Float_t tof, Float_t momentum, Float_t theta, 
                          Float_t phi, Float_t length, Float_t destep, 
                          Float_t Xref,Float_t Yref,Float_t Zref);
+    virtual void   AddHit(const AliMUONHit& ); // use copy constructor
+    
     virtual void   AddGlobalTrigger(Int_t *singlePlus, Int_t *singleMinus,
                                    Int_t *singleUndef, Int_t *pairUnlike, 
                                    Int_t *pairLike);
+    virtual void   AddGlobalTrigger(const AliMUONGlobalTrigger& trigger); // use copy constructor
+
     virtual void   AddLocalTrigger(Int_t* ltrigger);
+    virtual void   AddLocalTrigger(const AliMUONLocalTrigger& trigger); // use copy constructor
+
     virtual void   AddRawCluster(Int_t id, const AliMUONRawCluster& clust);
     virtual void   AddRecTrack(const AliMUONTrack& track);
 
index b97b51d29ed5602f47d25a3e82495f5827ebea4b..cde04b760a88478d4d74f73f1583e27559303b08 100644 (file)
 #include "AliMUONDigit.h"
 
 ClassImp(AliMUONDigit)
+//_____________________________________________________________________________
+ AliMUONDigit::AliMUONDigit(const AliMUONDigit& digits):TObject(digits)
+{
+// copy constructor
+  
+    fPadX        = digits.fPadX;
+    fPadY        = digits.fPadY;
+    fCathode     = digits.fCathode;
+    fSignal      = digits.fSignal;
+    fPhysics     = digits.fPhysics;
+    fHit         = digits.fHit;
+
+    for(Int_t i=0; i<kMAXTRACKS; i++) {
+       fTcharges[i]  = digits.fTcharges[i];
+       fTracks[i]    = digits.fTracks[i];
+    }
+}
+
 //_____________________________________________________________________________
 AliMUONDigit::AliMUONDigit(Int_t *digits)
 {
@@ -55,3 +73,24 @@ AliMUONDigit::~AliMUONDigit()
 {
     // Destructor 
 }
+
+//_____________________________________________________________________________
+AliMUONDigit& AliMUONDigit::operator=(const AliMUONDigit& digits)
+{
+    if (this == &digits)
+      return *this;
+    fPadX        = digits.fPadX;
+    fPadY        = digits.fPadY;
+    fCathode     = digits.fCathode;
+    fSignal      = digits.fSignal;
+    fPhysics     = digits.fPhysics;
+    fHit         = digits.fHit;
+
+    for(Int_t i=0; i<kMAXTRACKS; i++) {
+       fTcharges[i]  = digits.fTcharges[i];
+       fTracks[i]    = digits.fTracks[i];
+    }
+
+    return *this;
+}
index a0331c5677ef192be64ed6dcf6cb59acfb180545..13427e488f5afece28e58b836289ca03b2baab6e 100644 (file)
@@ -13,9 +13,12 @@ class AliMUONDigit : public TObject {
 
  public:
     AliMUONDigit() {}
+    AliMUONDigit(const AliMUONDigit& );
     AliMUONDigit(Int_t *digits);
     AliMUONDigit(Int_t *tracks, Int_t *charges, Int_t *digits);
     virtual ~AliMUONDigit();
+
+    AliMUONDigit& operator=(const AliMUONDigit& );
     
     virtual Int_t    PadX() const         {return fPadX;}
     virtual Int_t    PadY() const         {return fPadY;}
index 866467fdfe96eac2cb4e136aa7cd6bf00908858e..2a3e830045eaddd961c8c23a73b457f9e2d52698 100644 (file)
@@ -105,20 +105,31 @@ AliMUONTrack::~AliMUONTrack()
   //__________________________________________________________________________
 AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack):TObject(MUONTrack)
 {
-  fEventReconstructor =  MUONTrack.fEventReconstructor;
-  fTrackHitsPtr =  MUONTrack.fTrackHitsPtr;
+  fEventReconstructor = new AliMUONEventReconstructor(*MUONTrack.fEventReconstructor);
   fTrackParamAtVertex = MUONTrack.fTrackParamAtVertex;
+  fTrackHitsPtr =  new TObjArray(*MUONTrack.fTrackHitsPtr);
   fNTrackHits =  MUONTrack.fNTrackHits;
   fFitMCS     =  MUONTrack.fFitMCS;
   fFitNParam  =  MUONTrack.fFitNParam;
   fFitFMin    =  MUONTrack.fFitFMin;
+  fFitStart   =  MUONTrack.fFitStart;
 }
 
   //__________________________________________________________________________
-AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& /*MUONTrack*/)
+AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
 {
-// Dummy assignment operator
+  if (this == &MUONTrack)
     return *this;
+
+  fEventReconstructor = new AliMUONEventReconstructor(*MUONTrack.fEventReconstructor);
+  fTrackParamAtVertex = MUONTrack.fTrackParamAtVertex;
+  fTrackHitsPtr =  new TObjArray(*MUONTrack.fTrackHitsPtr);
+  fNTrackHits =  MUONTrack.fNTrackHits;
+  fFitMCS     =  MUONTrack.fFitMCS;
+  fFitNParam  =  MUONTrack.fFitNParam;
+  fFitFMin    =  MUONTrack.fFitFMin;
+  fFitStart   =  MUONTrack.fFitStart;
+  return *this;
 }
 
   //__________________________________________________________________________
index b3326209dc9af7f79fa1ce246a65ae0531dee3e5..a64ed48ae9e192c03c922ca98ddb48f0dd2b74cf 100644 (file)
@@ -70,6 +70,32 @@ extern "C" {
     Field[0] = b[0]; Field[1] = b[1]; Field[2] = b[2];
   }
 }
+  //_________________________________________________________________________
+
+AliMUONTrackParam& AliMUONTrackParam::operator=(const AliMUONTrackParam& MUONTrackParam)
+{
+  if (this == &MUONTrackParam)
+    return *this;
+
+  fInverseBendingMomentum =  MUONTrackParam.fInverseBendingMomentum; 
+  fBendingSlope           =  MUONTrackParam.fBendingSlope; 
+  fNonBendingSlope        =  MUONTrackParam.fNonBendingSlope; 
+  fZ                      =  MUONTrackParam.fZ; 
+  fBendingCoor            =  MUONTrackParam.fBendingCoor; 
+  fNonBendingCoor         =  MUONTrackParam.fNonBendingCoor;
+
+  return *this;
+}
+  //_________________________________________________________________________
+AliMUONTrackParam::AliMUONTrackParam(const AliMUONTrackParam& MUONTrackParam):TObject(MUONTrackParam)
+{
+  fInverseBendingMomentum =  MUONTrackParam.fInverseBendingMomentum; 
+  fBendingSlope           =  MUONTrackParam.fBendingSlope; 
+  fNonBendingSlope        =  MUONTrackParam.fNonBendingSlope; 
+  fZ                      =  MUONTrackParam.fZ; 
+  fBendingCoor            =  MUONTrackParam.fBendingCoor; 
+  fNonBendingCoor         =  MUONTrackParam.fNonBendingCoor;
+}
 
   //__________________________________________________________________________
 void AliMUONTrackParam::ExtrapToZ(Double_t Z)
index 5ba7d9a40fd9b1436ef15cce9d68a256066bbd41..1e027ce5406790d29b9cee3359aa775d2e6765c7 100644 (file)
@@ -19,7 +19,9 @@ class AliMUONTrackParam : public TObject {
   virtual ~AliMUONTrackParam(){
     // Destructor
     ;} // Destructor
-
+  
+  AliMUONTrackParam(const AliMUONTrackParam& );// copy constructor (should be added per default !)
+  AliMUONTrackParam& operator=(const  AliMUONTrackParam& );// (should be added per default !)
   // Get and Set methods for data
   Double_t GetInverseBendingMomentum(void) const {return fInverseBendingMomentum;}
   void SetInverseBendingMomentum(Double_t InverseBendingMomentum) {fInverseBendingMomentum = InverseBendingMomentum;}
@@ -38,7 +40,8 @@ class AliMUONTrackParam : public TObject {
   void ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam);
   void ExtrapToVertex();  // extrapolation to vertex through the absorber
   void BransonCorrection(); // makes Branson correction
-  Double_t TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta); // returns total momentum after energy loss correction in the absorber
+  // returns total momentum after energy loss correction in the absorber
+  Double_t TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta);
   void FieldCorrection(Double_t Z); // makes simple magnetic field correction through the absorber 
 
  protected:
diff --git a/MUON/MUONmassPlot_NewIO.C b/MUON/MUONmassPlot_NewIO.C
new file mode 100644 (file)
index 0000000..fc7a992
--- /dev/null
@@ -0,0 +1,266 @@
+// ROOT includes
+#include "TBranch.h"
+#include "TClonesArray.h"
+#include "TLorentzVector.h"
+#include "TFile.h"
+#include "TH1.h"
+#include "TParticle.h"
+#include "TTree.h"
+
+// STEER includes
+#include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliHeader.h"
+#include "AliLoader.h"
+#include "AliStack.h"
+
+// MUON includes
+#include "AliMUON.h"
+#include "AliMUONData.h"
+#include "AliMUONHit.h"
+#include "AliMUONConstants.h"
+#include "AliMUONDigit.h"
+#include "AliMUONRawCluster.h"
+#include "AliMUONGlobalTrigger.h"
+#include "AliMUONLocalTrigger.h"
+#include "AliMUONTrack.h"
+#include "AliMUONTrackParam.h"
+#include "AliESDMuonTrack.h"
+
+//
+// Macro MUONmassPlot.C for new I/O
+// Ch. Finck, Subatech, Jan. 2004
+//
+
+// macro to make invariant mass plots
+// for combinations of 2 muons with opposite charges,
+// from root file "MUON.tracks.root" containing the result of track reconstruction.
+// Histograms are stored on the "MUONmassPlot.root" file.
+// introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...)
+// using Invariant Mass for rapidity.
+
+// Arguments:
+//   FirstEvent (default 0)
+//   LastEvent (default 0)
+//   ResType (default 553)
+//      553 for Upsilon, anything else for J/Psi
+//   Chi2Cut (default 100)
+//      to keep only tracks with chi2 per d.o.f. < Chi2Cut
+//   PtCutMin (default 1)
+//      to keep only tracks with transverse momentum > PtCutMin
+//   PtCutMax (default 10000)
+//      to keep only tracks with transverse momentum < PtCutMax
+//   massMin (default 9.17 for Upsilon) 
+//      &  massMax (default 9.77 for Upsilon) 
+//         to calculate the reconstruction efficiency for resonances with invariant mass
+//         massMin < mass < massMax.
+
+// Add parameters and histograms for analysis 
+
+void MUONmassPlot(char* filename="galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, 
+                  Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,
+                  Float_t massMin = 9.17,Float_t massMax = 9.77)
+{
+  cout << "MUONmassPlot " << endl;
+  cout << "FirstEvent " << FirstEvent << endl;
+  cout << "LastEvent " << LastEvent << endl;
+  cout << "ResType " << ResType << endl;
+  cout << "Chi2Cut " << Chi2Cut << endl;
+  cout << "PtCutMin " << PtCutMin << endl;
+  cout << "PtCutMax " << PtCutMax << endl;
+  cout << "massMin " << massMin << endl;
+  cout << "massMax " << massMax << endl;
+
+
+  //Reset ROOT and connect tree file
+  gROOT->Reset();
+
+
+  // File for histograms and histogram booking
+  TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
+  TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
+  TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
+  TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
+  TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
+  TH1F *hInvMassRes;
+
+  if (ResType == 553) {
+    hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
+  } else {
+    hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
+  }
+
+  TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
+  TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
+  TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
+  TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
+
+
+  // settings
+  Int_t EventInMass = 0;
+  Float_t muonMass = 0.105658389;
+//   Float_t UpsilonMass = 9.46037;
+//   Float_t JPsiMass = 3.097;
+
+  Double_t bendingSlope, nonBendingSlope, pYZ;
+  Double_t fPxRec1, fPyRec1, fPzRec1, fZRec1, fE1;
+  Double_t fPxRec2, fPyRec2, fPzRec2, fZRec2, fE2;
+  Int_t fCharge, fCharge2;
+
+  Int_t ntrackhits, nevents;
+  Double_t fitfmin;
+
+  TClonesArray * recTracksArray;
+  TLorentzVector fV1, fV2, fVtot;
+  
+  // Creating Run Loader and openning file containing Hits
+  AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
+  if (RunLoader == 0x0) {
+    printf(">>> Error : Error Opening %s file \n",filename);
+    return;
+  }
+  
+  AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
+  MUONLoader->LoadTracks("READ");
+
+  // Creating MUON data container
+  AliMUONData muondata(MUONLoader,"MUON","MUON");
+  
+  nevents = RunLoader->GetNumberOfEvents();
+  
+  AliMUONTrack * rectrack;
+  AliMUONTrackParam *trackParam;
+        
+  // Loop over events
+  for (Int_t ievent = FirstEvent; ievent <= TMath::Min(LastEvent, nevents - 1); ievent++) {
+
+    // get current event
+    RunLoader->GetEvent(ievent);
+    
+    muondata.SetTreeAddress("RT");
+    muondata.GetRecTracks();
+    recTracksArray = muondata.RecTracks();
+
+    Int_t nrectracks = (Int_t) recTracksArray->GetEntriesFast(); //
+
+    printf("\n Nb of events analysed: %d\r",ievent);
+    //   cout << " number of tracks: " << nrectracks  <<endl;
+  
+    // loop over all reconstructed tracks (also first track of combination)
+     for (Int_t irectracks = 0; irectracks <  nrectracks;  irectracks++) {
+
+      rectrack = (AliMUONTrack*) recTracksArray->At(irectracks);
+
+      trackParam = rectrack->GetTrackParamAtVertex();
+      bendingSlope   = trackParam->GetBendingSlope();
+      nonBendingSlope = trackParam->GetNonBendingSlope();
+
+      pYZ     = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
+      fPzRec1  = - pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope); // spectro. (z<0)
+      fPxRec1  = fPzRec1 * nonBendingSlope;
+      fPyRec1  = fPzRec1 * bendingSlope;
+      fZRec1   = trackParam->GetZ();
+      fCharge = Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
+
+      fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
+      fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
+
+      ntrackhits = rectrack->GetNTrackHits();
+      fitfmin = rectrack->GetFitFMin();
+
+      // transverse momentum
+      Float_t pt1 = fV1.Pt();
+
+      // total momentum
+      Float_t p1 = fV1.P();
+
+      // Rapidity
+      Float_t rapMuon1 = fV1.Rapidity();
+
+      // chi2 per d.o.f.
+      Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
+//      printf(" px %f py %f pz %f NHits %d  Norm.chi2 %f charge %d\n", 
+//          fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge);
+
+      // condition for good track (Chi2Cut and PtCut)
+
+      if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
+
+       // fill histos hPtMuon and hChi2PerDof
+       hPtMuon->Fill(pt1);
+       hPMuon->Fill(p1);
+       hChi2PerDof->Fill(ch1);
+       hRapMuon->Fill(rapMuon1);
+
+       // loop over second track of combination
+       for (Int_t irectracks2 = irectracks + 1; irectracks2 < nrectracks; irectracks2++) {
+
+         rectrack = (AliMUONTrack*) recTracksArray->At(irectracks2);
+        
+         trackParam = rectrack->GetTrackParamAtVertex();
+         bendingSlope    = trackParam->GetBendingSlope();
+         nonBendingSlope = trackParam->GetNonBendingSlope();
+         
+         pYZ      = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
+         fPzRec2  = - pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope); // spectro. (z<0)
+         fPxRec2  = fPzRec2 * nonBendingSlope;
+         fPyRec2  = fPzRec2 * bendingSlope;
+         fZRec2   = trackParam->GetZ();
+         fCharge2 = Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
+
+         fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
+         fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
+
+         ntrackhits = rectrack->GetNTrackHits();
+         fitfmin = rectrack->GetFitFMin();
+
+         // transverse momentum
+         Float_t pt2 = fV2.Pt();
+
+         // chi2 per d.o.f.
+         Float_t ch2 = fitfmin  / (2.0 * ntrackhits - 5);
+
+         // condition for good track (Chi2Cut and PtCut)
+         if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
+
+           // condition for opposite charges
+           if ((fCharge * fCharge2) == -1) {
+
+             // invariant mass
+             fVtot = fV1 + fV2;
+             Float_t invMass = fVtot.M();
+                   
+             // fill histos hInvMassAll and hInvMassRes
+             hInvMassAll->Fill(invMass);
+             hInvMassRes->Fill(invMass);
+
+             if (invMass > massMin && invMass < massMax) {
+               EventInMass++;
+               hRapResonance->Fill(fVtot.Rapidity());
+               hPtResonance->Fill(fVtot.Pt());
+             }
+
+           } // if (fCharge * fCharge2) == -1)
+         } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
+       } //  for (Int_t irectracks2 = irectracks + 1; irectracks2 < irectracks; irectracks2++)
+      } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
+    } // for (Int_t irectracks = 0; irectracks < nrectracks; irectracks++)
+
+    hNumberOfTrack->Fill(nrectracks);
+  } // for (Int_t ievent = FirstEvent;
+
+  histoFile->Write();
+  histoFile->Close();
+
+  cout << "MUONmassPlot " << endl;
+  cout << "FirstEvent " << FirstEvent << endl;
+  cout << "LastEvent " << LastEvent << endl;
+  cout << "ResType " << ResType << endl;
+  cout << "Chi2Cut " << Chi2Cut << endl;
+  cout << "PtCutMin " << PtCutMin << endl;
+  cout << "PtCutMax " << PtCutMax << endl;
+  cout << "massMin " << massMin << endl;
+  cout << "massMax " << massMax << endl;
+  cout << "EventInMass " << EventInMass << endl;
+}
+