ReaderESDtree, MUON analysis, reading MUON data froESD in ReaderESD (Christian FINCK)
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 11 Aug 2004 13:15:41 +0000 (13:15 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 11 Aug 2004 13:15:41 +0000 (13:15 +0000)
ANALYSIS/ANALYSISLinkDef.h
ANALYSIS/AliMuonAnalysis.cxx [new file with mode: 0644]
ANALYSIS/AliMuonAnalysis.h [new file with mode: 0644]
ANALYSIS/AliReaderAOD.cxx
ANALYSIS/AliReaderAOD.h
ANALYSIS/AliReaderESD.cxx
ANALYSIS/AliReaderESD.h
ANALYSIS/AliReaderESDTree.cxx [new file with mode: 0644]
ANALYSIS/AliReaderESDTree.h [new file with mode: 0644]
ANALYSIS/WriteAOD.C
ANALYSIS/libANALYSIS.pkg

index 19943e6c05b819628b4f66ffe9821292addf0a22..fd0ffddaf3627131a72be24246ff04d3b2418b48 100644 (file)
 
 #pragma link C++ class AliReader+;
 #pragma link C++ class AliReaderESD+;
-#pragma link C++ class AliReaderAOD+;
+#pragma link C++ class AliReaderESDTree+;
 #pragma link C++ class AliReaderKineTree+;
 
 #pragma link C++ class AliFlowAnalysis+;
+#pragma link C++ class AliMuonAnalysis+;
 
 #pragma link C++ class AliEventCut+;
 #pragma link C++ class AliEventEmptyCut+;
diff --git a/ANALYSIS/AliMuonAnalysis.cxx b/ANALYSIS/AliMuonAnalysis.cxx
new file mode 100644 (file)
index 0000000..154df50
--- /dev/null
@@ -0,0 +1,167 @@
+#include "AliMuonAnalysis.h"
+//________________________________
+///////////////////////////////////////////////////////////
+//
+// class AliMuonAnalysis
+//
+// MUON Analysis
+//
+//
+// 
+// finck@subatech.in2p3.fr
+//
+///////////////////////////////////////////////////////////
+/*********************************************************/
+
+#include <TString.h>
+#include <TParticle.h>
+
+#include <AliStack.h>
+#include <AliAOD.h>
+#include <AliAODParticle.h>
+#include <AliAODParticleCut.h>
+
+#include <AliESDtrack.h>
+#include <AliESD.h>
+
+#include "TFile.h"
+#include "TH1.h"
+#include "TH2.h"
+
+ClassImp(AliMuonAnalysis)
+
+AliMuonAnalysis::AliMuonAnalysis():
+ fPartCut(0x0)
+{
+ //ctor
+}
+
+/*********************************************************/
+AliMuonAnalysis::~AliMuonAnalysis()
+{
+ //dtor
+  delete fPartCut;
+  delete fHistoFile;
+  delete fHPtMuon;
+  delete fHPtMuonPlus;
+  delete fHPtMuonMinus;
+  delete fHPMuon;
+  delete fHInvMassAll;
+  delete fHRapMuon;
+  delete fHRapResonance;
+  delete fHPtResonance;
+  delete fHInvMassAll_vs_Pt;
+}
+/*********************************************************/
+
+Int_t AliMuonAnalysis::Init()
+{
+  //Initilizes anaysis
+  Info("Init","Histo initialized for MUON Analysis");
+
+  fHistoFile         = new TFile("MUONmassPlot.root", "RECREATE");
+  fHPtMuon           = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
+  fHPMuon            = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
+  fHPtMuonPlus       = new TH1F("hPtMuonPlus", "Muon+ Pt (GeV/c)", 100, 0., 20.);
+  fHPtMuonMinus      = new TH1F("hPtMuonMinus", "Muon- Pt (GeV/c)", 100, 0., 20.);
+  fHInvMassAll       = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
+  fHRapMuon          = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
+  fHRapResonance     = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
+  fHPtResonance      = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
+  fHInvMassAll_vs_Pt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);
+
+  return 0; 
+}
+/*********************************************************/
+
+Int_t AliMuonAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
+{
+  if (aodrec) {
+    GetInvMass(aodrec);
+    //    Info("ProcessEvent","Inv Mass Rec");
+  }  
+
+  if (aodsim) {
+    //     Info("ProcessEvent","aodsim not implemented");
+  }  
+  
+  return 0;
+  
+}
+
+/*********************************************************/
+
+Int_t AliMuonAnalysis::Finish()
+{
+  //Finish analysis and writes results
+  Info("Finish","Histo writing for MUON Analysis");
+
+  fHistoFile->Write();
+  fHistoFile->Close();
+
+  return 0;
+}
+/*********************************************************/
+
+void AliMuonAnalysis::GetInvMass(AliAOD* aod)
+{
+
+  TLorentzVector lorV1, lorV2, lorVtot;
+  Float_t massMin = 9.17;
+  Float_t massMax = 9.77;
+  Int_t charge1, charge2;
+
+//returns flow parameters: v2 and event plane
+  if (aod == 0x0) {
+     Error("AliMuonAnalysis::GetInvMass","Pointer to AOD is NULL");
+     return;
+  }
+   
+  Int_t nPart = aod->GetNumberOfParticles();
+  
+  for (Int_t iPart1 = 0; iPart1 < nPart; iPart1++)  {
+    AliAODParticle* aodPart1 =  (AliAODParticle*)aod->GetParticle(iPart1);
+
+    if (aodPart1 == 0x0) {
+      Error("AliMuonAnalysis::GetInvMass","Cannot get particle %d", iPart1);
+      continue;
+    }
+
+    lorV1 = aodPart1->FourMomentum();
+
+    fHPtMuon->Fill(lorV1.Pt());
+    fHPMuon->Fill(lorV1.P());
+
+    charge1 = TMath::Sign(1,aodPart1->GetPdgCode());
+
+    if (charge1 > 0) {
+      fHPtMuonPlus->Fill(lorV1.Pt());
+    } else {
+      fHPtMuonMinus->Fill(lorV1.Pt());
+    }
+    fHRapMuon->Fill(lorV1.Rapidity());
+    for (Int_t iPart2 = iPart1 + 1; iPart2 < nPart; iPart2++)  {
+
+      AliAODParticle* aodPart2 =  (AliAODParticle*)aod->GetParticle(iPart2);
+
+      lorV2 = aodPart2->FourMomentum();
+      charge2 = TMath::Sign(1,aodPart2->GetPdgCode());
+
+      if ((charge1 * charge2) == -1) {
+
+       lorVtot = lorV1 + lorV2;
+       Float_t invMass = lorVtot.M();
+
+       fHInvMassAll->Fill(invMass);
+       fHInvMassAll_vs_Pt->Fill(invMass,lorVtot.Pt());
+
+       if (invMass > massMin && invMass < massMax) {
+         fHRapResonance->Fill(lorVtot.Rapidity());
+         fHPtResonance->Fill(lorVtot.Pt());
+       }
+      }
+
+    }
+  }
+}
diff --git a/ANALYSIS/AliMuonAnalysis.h b/ANALYSIS/AliMuonAnalysis.h
new file mode 100644 (file)
index 0000000..97bb6dc
--- /dev/null
@@ -0,0 +1,59 @@
+#ifndef ALIMUONANALYSIS_H
+#define ALIMUONANALYSIS_H
+//________________________________
+///////////////////////////////////////////////////////////
+//
+// class AliMuonAnalysis
+//
+// Flow Analysis
+//
+//
+// S.Radomski@gsi.de
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////
+
+#include "AliAnalysis.h"
+
+class AliESD;
+class AliAOD;
+class AliStack;
+class AliAODParticleCut;
+class TFile;
+class TH1F;
+class TH2F;
+
+class AliMuonAnalysis: public AliAnalysis
+{ 
+  public: 
+     AliMuonAnalysis();
+     virtual ~AliMuonAnalysis();
+
+    Int_t Init();
+    Int_t ProcessEvent(AliAOD* aodrec, AliAOD* aodsim);
+    Int_t Finish();
+   
+    void SetParticleCut(AliAODParticleCut* pcut){fPartCut = pcut;}
+
+    void GetInvMass(AliAOD* aod);
+
+  protected:
+    
+  private:
+
+    TFile *fHistoFile;
+    TH1F *fHPtMuon;
+    TH1F *fHPtMuonPlus;
+    TH1F *fHPtMuonMinus;
+    TH1F *fHPMuon;
+    TH1F *fHInvMassAll;
+    TH1F *fHRapMuon;
+    TH1F *fHRapResonance;
+    TH1F *fHPtResonance;
+    TH2F *fHInvMassAll_vs_Pt;
+
+    AliAODParticleCut* fPartCut;//Particle Cut
+    ClassDef(AliMuonAnalysis,1)
+};
+
+#endif
index 38e884f9da769b12a1db8177eb70f304d1f732f8..6575d6c6acdf38fbfa94710c5cde92c3d96228b6 100644 (file)
@@ -7,7 +7,7 @@ ClassImp(AliReaderAOD)
 #include <TTree.h>
 #include "AliAOD.h"
 
-Int_t AliReaderAOD::WriteAOD(AliReader* reader, const char* outfilename, Bool_t /*multcheck*/)
+Int_t AliReaderAOD::WriteAOD(AliReader* reader, const char* outfilename, const char* pclassname,  Bool_t /*multcheck*/)
 {
 //reads tracks from runs and writes them to file
   ::Info("AliReaderAOD::Write","________________________________________________________");
@@ -29,33 +29,38 @@ Int_t AliReaderAOD::WriteAOD(AliReader* reader, const char* outfilename, Bool_t
   TTree *tree = new TTree("TAOD","Tree with tracks");
   
   TBranch *recbranch = 0x0, *simbranch = 0x0;
-
-  AliAOD* eventsim = new AliAOD();
-  AliAOD* eventrec = new AliAOD;
   
-  eventsim->SetParticleClassName("AliAODParticle");
-  eventrec->SetParticleClassName("AliAODParticle");
   
-  if (reader->ReadsSim()) simbranch = tree->Branch("simulated","AliAOD",&eventsim,32000,99);
+  AliAOD* eventrec = new AliAOD();//must be created before Branch is called. Otherwise clones array is not splitted
+  AliAOD* eventsim = new AliAOD();//AOD together with fParticles clones array knowing exact type of particles
+  
+  eventrec->SetParticleClassName(pclassname);
+  eventsim->SetParticleClassName(pclassname);
+  
   if (reader->ReadsRec()) recbranch = tree->Branch("reconstructed","AliAOD",&eventrec,32000,99);
+  if (reader->ReadsSim()) simbranch = tree->Branch("simulated","AliAOD",&eventsim,32000,99);
 
+  delete eventsim;
+  delete eventrec;
+  
   reader->Rewind();
   while (reader->Next() == kFALSE)
    {
      
-     if (reader->ReadsSim())
-      {
-        eventsim = reader->GetEventSim();
-//        simbranch->SetAddress(&eventsim);
-      }
      if (reader->ReadsRec())
       {
         eventrec = reader->GetEventRec();
-//        recbranch->SetAddress(&eventrec);
+        recbranch->SetAddress(&eventrec);
+      }
+
+     if (reader->ReadsSim())
+      {
+        eventsim = reader->GetEventSim();
+        simbranch->SetAddress(&eventsim);
       }
+     eventrec->GetParticle(0)->Print();
+     eventsim->GetParticle(0)->Print();
      tree->Fill();
-     tree->Print();
    }
   
   ::Info("AliReaderAOD::Write","Written %d events",tree->GetEntries());
index 7ccd47a6719760cb9a15c3b8f1077e47f5b9cdee..52f8b340bc884ca2449d51f96683a7a65bec8cac 100644 (file)
@@ -14,7 +14,8 @@ class AliReaderAOD: public AliReader
     Bool_t        ReadsSim() const {return fReadSim;}
 
 
-    static Int_t WriteAOD(AliReader* reader, const char* outfilename = "AliAOD.root", Bool_t multcheck = kFALSE);//reads tracks from runs and writes them to file
+    static Int_t WriteAOD(AliReader* reader, const char* outfilename = "AliAOD.root", //reads tracks from runs and writes them to file
+                          const char* pclassname = "AliAODParticle", Bool_t multcheck = kFALSE);
     
   protected:
   private:
index 17e09a737d065511382c807bfc2a584a3a7cd28f..f49f9483acccd63f7165df7e7f20e4d4f09c527a 100644 (file)
@@ -50,6 +50,9 @@ AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename)
  fClusterMap(kFALSE),
  fITSTrackPoints(kFALSE),
  fMustTPC(kFALSE),
+ fReadCentralBarrel(kFALSE),
+ fReadMuon(kFALSE),
+ fReadPHOS(kFALSE),
  fNTPCClustMin(0),
  fNTPCClustMax(150),
  fTPCChi2PerClustMin(0.0),
@@ -99,6 +102,9 @@ AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char
  fClusterMap(kFALSE),
  fITSTrackPoints(kFALSE),
  fMustTPC(kFALSE),
+ fReadCentralBarrel(kFALSE),
+ fReadMuon(kFALSE),
+ fReadPHOS(kFALSE),
  fNTPCClustMin(0),
  fNTPCClustMax(150),
  fTPCChi2PerClustMin(0.0),
@@ -135,10 +141,11 @@ AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char
 AliReaderESD::~AliReaderESD()
 {
  //desctructor
-  delete fRunLoader;
-  delete fKeyIterator;
-  delete fFile;
+    delete fRunLoader;
+    delete fKeyIterator;
+    delete fFile;
 }
+
 /**********************************************************/
 Int_t AliReaderESD::ReadNext()
 {
@@ -155,89 +162,67 @@ Int_t AliReaderESD::ReadNext()
   fEventRec->Reset();
         
   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
-   {
-     if (fFile == 0x0)
-      {
-       fFile = OpenFile(fCurrentDir);//rl is opened here
-       if (fFile == 0x0)
-         {
-           Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
-           fCurrentDir++;
-           continue;
-         }
-       fCurrentEvent = 0;
-       fKeyIterator = new TIter(fFile->GetListOfKeys());  
-//       fFile->Dump();
-//       fFile->GetListOfKeys()->Print();
-      } 
-     TKey* key = (TKey*)fKeyIterator->Next();
-     if (key == 0x0)
-      {
-        if (AliVAODParticle::GetDebug() > 2 )
-          {
-            Info("ReadNext","No more keys.");
-          }
-        fCurrentDir++;
-        delete fKeyIterator;
-        fKeyIterator = 0x0;
-        delete fFile;//we have to assume there is no more ESD objects in the fFile
-        fFile = 0x0;
-        delete fRunLoader;
-        fRunLoader = 0x0;
-        continue;
-      }
-     //try to read
-     
-     
-//     TObject* esdobj = key->ReadObj();
-//     if (esdobj == 0x0)
-//      {
-//        if (AliVAODParticle::GetDebug() > 2 )
-//          {
-//            Info("ReadNext","Key read NULL. Key Name is %s",key->GetName());
-//            key->Dump();
-//          }
-//        continue;
-//      }
-//     esdobj->Dump();
-//     AliESD* esd = dynamic_cast<AliESD*>(esdobj);
-     
+    {
+      if (fFile == 0x0)
+       {
+         fFile = OpenFile(fCurrentDir);//rl is opened here
+         if (fFile == 0x0)
+           {
+             Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
+             fCurrentDir++;
+             continue;
+           }
+         fCurrentEvent = 0;
+       }
      TString esdname = "ESD";
      esdname+=fCurrentEvent;
      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
      if (esd == 0x0)
       {
-//        if (AliVAODParticle::GetDebug() > 2 )
-//          {
-//            Info("ReadNext","This key is not an AliESD object %s",key->GetName());
-//          }
         if (AliVAODParticle::GetDebug() > 2 )
           {
             Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
           }
         fCurrentDir++;
-        delete fKeyIterator;
-        fKeyIterator = 0x0;
         delete fFile;//we have to assume there is no more ESD objects in the fFile
         fFile = 0x0;
         delete fRunLoader;
         fRunLoader = 0x0;
         continue;
       }
-     
-     ReadESD(esd);
+      ReadESD(esd);
       
-     fCurrentEvent++;
-     fNEventsRead++;
-     delete esd;
-     return 0;//success -> read one event
-   }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
+      fCurrentEvent++;
+      fNEventsRead++;
+      delete esd;
+      return 0;//success -> read one event
+    }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
    
   return 1; //no more directories to read
 }
-/**********************************************************/
 
+/**********************************************************/
 Int_t AliReaderESD::ReadESD(AliESD* esd)
+{
+//Reads esd data
+  if (esd == 0x0)
+   {
+     Error("ReadESD","ESD is NULL");
+     return 1;
+   }
+  
+  // seperate each method
+  if (fReadCentralBarrel) ReadESDCentral(esd);
+
+  if (fReadMuon) ReadESDMuon(esd);
+
+  if (fReadPHOS) ReadESDPHOS(esd);
+
+  return 1;
+}
+
+/**********************************************************/
+Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
 {
   //****** Tentative particle type "concentrations"
   static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
@@ -248,12 +233,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
   Double_t pos[3];//position
   Double_t vertexpos[3];//vertex position
   //Reads one ESD
-  if (esd == 0x0)
-   {
-     Error("ReadESD","ESD is NULL");
-     return 1;
-   }
-  
+
   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
   if (pdgdb == 0x0)
    {
@@ -577,15 +557,89 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
   return 0;
 }
 
+/**********************************************************/
+Int_t AliReaderESD::ReadESDMuon(AliESD* esd)
+{
+
+  Double_t vertexpos[3];//vertex position, assuming no secondary decay
+
+  const AliESDVertex* vertex = esd->GetVertex();
+
+  if (vertex == 0x0) {
+    Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
+    vertexpos[0] = 0.0;
+    vertexpos[1] = 0.0;
+    vertexpos[2] = 0.0;
+  } else {
+    vertex->GetXYZ(vertexpos);
+  }
+
+ Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
+
+ if (AliVAODParticle::GetDebug() > 0) {
+   Info("ReadESD","Reading Event %d",fCurrentEvent);
+   Info("ReadESD","Found %d tracks.",nTracks);
+ }
+ // settings
+  Float_t Chi2Cut = 100.;
+  Float_t PtCutMin = 1.;
+  Float_t PtCutMax = 10000.;
+  Float_t muonMass = 0.105658389;
+  Int_t pdgcode = -13;
+  Double_t thetaX, thetaY, pYZ;
+  Double_t pxRec1, pyRec1, pzRec1, E1;
+  Int_t charge;
+
+  Int_t ntrackhits;
+  Double_t fitfmin;
+
+  TLorentzVector fV1;
+  fEventRec->Reset();
+  for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
+
+      AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
+
+      thetaX = muonTrack->GetThetaX();
+      thetaY = muonTrack->GetThetaY();
+
+      pYZ     =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
+      pzRec1  = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
+      pxRec1  = pzRec1 * TMath::Tan(thetaX);
+      pyRec1  = pzRec1 * TMath::Tan(thetaY);
+      charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
+      E1 = TMath::Sqrt(muonMass * muonMass + pxRec1 * pxRec1 + pyRec1 * pyRec1 + pzRec1 * pzRec1);
+      fV1.SetPxPyPzE(pxRec1, pyRec1, pzRec1, E1);
+
+      ntrackhits = muonTrack->GetNHit();
+      fitfmin    = muonTrack->GetChi2();
+
+      // transverse momentum
+      Float_t pt1 = fV1.Pt();
+
+      // chi2 per d.o.f.
+      Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
+
+      if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
+       AliAODParticle* track = new AliAODParticle(pdgcode*charge,1,iTrack, 
+                                                   pxRec1, pyRec1,pzRec1, E1,
+                                                   vertexpos[0], vertexpos[1], vertexpos[2], 0.);
+        fEventRec->AddParticle(track);
+      }
+  }
+  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
+  return 0;
+}
+
 /**********************************************************/
 
 void AliReaderESD::Rewind()
 {
   //rewinds reading 
-  delete fKeyIterator;
+  //  delete fKeyIterator;
   delete fFile;
   fFile = 0x0;
-  fKeyIterator = 0x0;
+  // fKeyIterator = 0x0;
   delete fRunLoader;
   fRunLoader = 0x0;
   fCurrentDir = 0;
@@ -596,7 +650,7 @@ void AliReaderESD::Rewind()
 
 TFile* AliReaderESD::OpenFile(Int_t n)
 {
-//opens fFile with kine tree
+//opens fFile with  tree
 
  const TString& dirname = GetDirName(n);
  if (dirname == "")
index 2636a9b46b497447d340c2dbeb393abacb2efd74..fb5f72bb17b86d794e8716c5b00a7eae8bbdd634 100644 (file)
@@ -54,15 +54,23 @@ class AliReaderESD: public AliReader
     void          SetITSTrackPoints(Bool_t flag = kTRUE){fITSTrackPoints = flag;}
     void          MustTPC(Bool_t flag){fMustTPC = flag;}
 
-    
+    void          SetReadCentralBarrel(Bool_t flag){fReadCentralBarrel = flag;}
+    void          SetReadMuon(Bool_t flag){fReadMuon = flag;}
+    void          SetReadPHOS(Bool_t flag){fReadPHOS = flag;}
+
     enum ESpecies {kESDElectron = 0, kESDMuon, kESDPion, kESDKaon, kESDProton, kNSpecies};
     static Int_t  GetSpeciesPdgCode(ESpecies spec);
     
     Int_t         ReadESD(AliESD* esd);
-    
+    Int_t         ReadESDCentral(AliESD* esd);
+    Int_t         ReadESDMuon(AliESD* esd);
+    Int_t         ReadESDPHOS(AliESD* /*esd*/){return 0;}
+
   protected:
-    Int_t         ReadNext();
-    TFile*        OpenFile(Int_t evno);//opens files to be read for given event
+    virtual Int_t         ReadNext();
+    virtual TFile*        OpenFile(Int_t evno);//opens files to be read for given event
+
     Bool_t        CheckTrack(AliESDtrack* t) const;
     
     TString       fESDFileName;//name of the file with tracks
@@ -85,7 +93,11 @@ class AliReaderESD: public AliReader
                      //used by anti-merging cut in non-id analysis
 
     Bool_t        fMustTPC;// must be reconstructed in TPC -> reject tracks reconstructed ITS stand alone
-    
+
+    Bool_t        fReadCentralBarrel; // Flag for reading ESD central track 
+    Bool_t        fReadMuon;// Flag for reading ESD Muon track 
+    Bool_t        fReadPHOS;// Flag for reading ESD Phos 
+
     //Cut Parameters specific to TPC tracks
         
     Int_t         fNTPCClustMin;//Number of clusters min value
diff --git a/ANALYSIS/AliReaderESDTree.cxx b/ANALYSIS/AliReaderESDTree.cxx
new file mode 100644 (file)
index 0000000..6c7650c
--- /dev/null
@@ -0,0 +1,153 @@
+#include "AliReaderESDTree.h"
+//_______________________________________________________________________
+/////////////////////////////////////////////////////////////////////////
+//
+// class AliReaderESDTree
+//
+// Reader for MUON ESD Tree (only for rec)
+//
+// finck@subatech.in2p3.fr
+//
+/////////////////////////////////////////////////////////////////////////
+
+#include <TString.h>
+#include <TTree.h>
+#include <TFile.h>
+
+
+#include <AliRun.h>
+#include <AliRunLoader.h>
+
+#include <AliESD.h>
+#include "AliAOD.h"
+
+ClassImp(AliReaderESDTree)
+
+AliReaderESDTree::AliReaderESDTree(const Char_t* esdfilename, const Char_t* galfilename):
+  AliReaderESD(esdfilename,galfilename),
+  fTree(0x0)
+{
+//ctor
+}
+
+/********************************************************************/
+AliReaderESDTree::~AliReaderESDTree()
+{
+//dtor 
+ delete fTree;
+}
+
+/**********************************************************/
+Int_t AliReaderESDTree::ReadNext()
+{
+//reads next event from fFile
+//fRunLoader is for reading Kine
+  
+  if (AliVAODParticle::GetDebug())
+    Info("ReadNext","Entered");
+    
+  if (fEventSim == 0x0)  fEventSim = new AliAOD();
+  if (fEventRec == 0x0)  fEventRec = new AliAOD();
+  
+  fEventSim->Reset();
+  fEventRec->Reset();
+        
+  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
+    {
+      if (fFile == 0x0)
+       {
+         fFile = OpenFile(fCurrentDir);//rl is opened here
+         if (fFile == 0x0)
+           {
+             Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
+             fCurrentDir++;
+             continue;
+           }
+         fCurrentEvent = 0;
+       }
+
+      static AliESD* esd = 0x0;
+      fTree->SetBranchAddress("ESD", &esd);
+      Int_t status = fTree->GetEvent(fCurrentEvent);
+
+      if (!status)
+       {
+         if (AliVAODParticle::GetDebug() > 2 )
+           {
+             Info("ReadNext","Can not find event# %d in Tree", fCurrentEvent);
+           }
+         fCurrentDir++;
+         delete fFile;//we have to assume there is no more ESD objects in the fFile
+         fFile = 0x0;
+         delete fRunLoader;
+         fRunLoader = 0x0;
+         continue;
+       }
+
+      ReadESD(esd);
+      
+      fCurrentEvent++;
+      fNEventsRead++;
+      return 0;//success -> read one event
+    }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
+   
+  return 1; //no more directories to read
+}
+
+/**********************************************************/
+TFile* AliReaderESDTree::OpenFile(Int_t n)
+{
+//opens fFile with kine tree
+
+ const TString& dirname = GetDirName(n);
+ if (dirname == "")
+  {
+   Error("OpenFiles","Can not get directory name");
+   return 0x0;
+  }
+ TString filename = dirname +"/"+ fESDFileName;
+ TFile *ret = TFile::Open(filename.Data()); 
+
+ if (ret == 0x0)
+  {
+    Error("OpenFiles","Can't open fFile %s",filename.Data());
+    return 0x0;
+  }
+ if (!ret->IsOpen())
+  {
+    Error("OpenFiles","Can't open fFile  %s",filename.Data());
+    return 0x0;
+  }
+ TString esdname = "esdTree";
+ fTree = dynamic_cast<TTree*> (ret->Get(esdname));
+
+ if (!fTree)
+  {
+    Error("OpenFiles","Can't open ESD Tree %s",esdname.Data());
+    return 0x0;
+
+  }
+ if (fReadSim )
+  {
+   fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
+   if (fRunLoader == 0x0)
+    {
+      Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
+      delete ret;
+      return 0x0;
+    }
+    
+   fRunLoader->LoadHeader();
+   if (fRunLoader->LoadKinematics())
+    {
+      Error("Next","Error occured while loading kinematics.");
+      delete fRunLoader;
+      delete ret;
+      return 0x0;
+    }
+  }
+   
+ return ret;
+}
diff --git a/ANALYSIS/AliReaderESDTree.h b/ANALYSIS/AliReaderESDTree.h
new file mode 100644 (file)
index 0000000..9c55b2b
--- /dev/null
@@ -0,0 +1,39 @@
+#ifndef AliReaderESDTree_H
+#define AliReaderESDTree_H
+//_______________________________________________________________________
+/////////////////////////////////////////////////////////////////////////
+//
+// class AliReaderESDTree
+//
+// Reader for ESD Tree 
+//
+// Ch. Finck
+//
+/////////////////////////////////////////////////////////////////////////
+#include "AliReaderESD.h"
+#include <TString.h>
+
+class TFile;
+class TTree;
+
+class AliReaderESDTree: public AliReaderESD
+ {
+   public:
+
+    AliReaderESDTree(const Char_t* esdfilename = "AliESDs.root", 
+                     const Char_t* galfilename = "galice.root");
+
+    virtual ~AliReaderESDTree();
+
+
+   protected:
+    Int_t         ReadNext();//reads tracks and particles and puts them in runs
+    TFile*        OpenFile(Int_t evno);//opens files to be read for given event
+   
+    TTree*        fTree;// tree pointer
+    
+   private:
+    ClassDef(AliReaderESDTree,1)
+ };
+
+#endif
index 6ac7008c562e3df26820526d32fab0d16f7ff134..bb85d39cb2477bf0b3fde0a4a41caabc9242e342 100644 (file)
@@ -95,15 +95,15 @@ void WriteAOD(Option_t* datatype, Option_t* processopt="TracksAndParticles",
    reader->SetDirs(dirs);
     
    AliAODParticleCut* readerpartcut= new AliAODParticleCut();
-   readerpartcut->SetPtRange(0.0,10000.0);
-   readerpartcut->SetPID(kKPlus);
-   AliAODPIDCut* pidcut = new AliAODPIDCut(kKPlus,0.5);
+   readerpartcut->SetPtRange(0.4,1.2);
+   readerpartcut->SetPID(kPiPlus);
+   AliAODPIDCut* pidcut = new AliAODPIDCut(kPiPlus,0.5);
    readerpartcut->AddBasePartCut(pidcut);
    
    reader->AddParticleCut(readerpartcut);//read this particle type with this cut
 
    cout<<"WriteAOD.C:   P R O C S E S S I N G .....\n\n";
-   AliReaderAOD::WriteAOD(reader,outfile,multcheck);
+   AliReaderAOD::WriteAOD(reader,outfile,"AliAODParticle",multcheck);
    cout<<"\n\nWriteAOD.C:   F I N I S H E D\n";
    
    if (dirs) delete dirs;
index 8085223f5e87225821c2ec072eddc9300d187cec..ab57524522e961ac4378c30ab9bb14cd7cd29ab8 100644 (file)
@@ -5,10 +5,10 @@ SRCS= AliAOD.cxx AliEventBuffer.cxx \
       AliAODParticleCut.cxx AliAODParticleBaseCut.cxx \
       AliAODPairCut.cxx AliAODPairBaseCut.cxx \
       AliEventCut.cxx AliEventBaseCut.cxx \
-      AliReader.cxx AliReaderESD.cxx AliReaderKineTree.cxx \
-      AliReaderAOD.cxx \
+      AliReader.cxx AliReaderESD.cxx AliReaderKineTree.cxx\
       AliTrackPoints.cxx AliClusterMap.cxx \
       AliD0toKpi.cxx  AliD0toKpiAnalysis.cxx AliFlowAnalysis.cxx \
+      AliReaderESDTree.cxx AliMuonAnalysis.cxx  \
 
 HDRS= $(SRCS:.cxx=.h) 
 
@@ -23,6 +23,7 @@ EXPORT:=AliAOD.h AliEventBuffer.h\
       AliEventCut.h         AliEventBaseCut.h \
       AliReader.h           AliReaderESD.h    \
       AliTrackPoints.h      AliClusterMap.h   \
-      AliFlowAnalysis.h 
+      AliFlowAnalysis.h     AliReaderESDTree.h \
+      AliMuonAnalysis.h
 
 EINCLUDE:= TPC CONTAINERS ITS