]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDInput.cxx
First V0 MC Analysis from H.Ricaud
[u/mrichter/AliRoot.git] / FMD / AliFMDInput.cxx
index a8facd7819a7af66cc93290c549a346044dd0e0a..ae5d30f5f288b6fdafdd831a4c4ae65a203075a5 100644 (file)
@@ -29,7 +29,7 @@
 // Latest changes by Christian Holm Christensen
 //
 #include "AliFMDInput.h"       // ALIFMDHIT_H
-#include "AliLog.h"            // ALILOG_H
+#include "AliFMDDebug.h"       // ALIFMDDEBUG_H ALILOG_H
 #include "AliLoader.h"          // ALILOADER_H
 #include "AliRunLoader.h"       // ALIRUNLOADER_H
 #include "AliRun.h"             // ALIRUN_H
 #include "AliFMDRawReader.h"    // ALIFMDRAWREADER_H
 #include <AliESD.h>
 #include <AliESDFMD.h>
+#include <AliESDEvent.h>
 #include <AliCDBManager.h>
 #include <AliCDBEntry.h>
-#include <AliAlignObjAngles.h>
+#include <AliAlignObjParams.h>
 #include <TTree.h>              // ROOT_TTree
 #include <TChain.h>             // ROOT_TChain
 #include <TParticle.h>          // ROOT_TParticle
@@ -55,6 +56,9 @@
 #include <TGeoManager.h>        // ROOT_TGeoManager 
 #include <TSystemDirectory.h>   // ROOT_TSystemDirectory
 #include <Riostream.h>         // ROOT_Riostream
+#include <TFile.h>              // ROOT_TFile
+#include <TStreamerInfo.h>
+#include <TArrayF.h>
 
 //____________________________________________________________________
 ClassImp(AliFMDInput)
@@ -72,8 +76,8 @@ AliFMDInput::AliFMDInput()
     fFMDLoader(0), 
     fReader(0),
     fFMD(0),
-    fMainESD(0),
     fESD(0),
+    fESDEvent(0),
     fTreeE(0),
     fTreeH(0),
     fTreeD(0),
@@ -91,6 +95,7 @@ AliFMDInput::AliFMDInput()
     fTreeMask(0), 
     fIsInit(kFALSE)
 {
+
   // Constructor of an FMD input object.  Specify what data to read in
   // using the AddLoad member function.   Sub-classes should at a
   // minimum overload the member function Event.   A full job can be
@@ -108,8 +113,8 @@ AliFMDInput::AliFMDInput(const char* gAliceFile)
     fFMDLoader(0), 
     fReader(0),
     fFMD(0),
-    fMainESD(0),
     fESD(0),
+    fESDEvent(0),
     fTreeE(0),
     fTreeH(0),
     fTreeD(0),
@@ -127,6 +132,7 @@ AliFMDInput::AliFMDInput(const char* gAliceFile)
     fTreeMask(0), 
     fIsInit(kFALSE)
 {
+  
   // Constructor of an FMD input object.  Specify what data to read in
   // using the AddLoad member function.   Sub-classes should at a
   // minimum overload the member function Event.   A full job can be
@@ -200,7 +206,10 @@ AliFMDInput::Init()
       TString fname(file->GetName());
       if (fname.Contains("AliESDs")) fChainE->AddFile(fname.Data());
     }
-    fChainE->SetBranchAddress("ESD", &fMainESD);
+    fESDEvent = new AliESDEvent();
+    fESDEvent->ReadFromTree(fChainE);
+    //    fChainE->SetBranchAddress("ESD", &fMainESD);
+    
   }
     
   if (TESTBIT(fTreeMask, kRaw)) {
@@ -233,17 +242,17 @@ AliFMDInput::Init()
       else {
        Int_t nAlign = array->GetEntries();
        for (Int_t i = 0; i < nAlign; i++) {
-         AliAlignObjAngles* a = static_cast<AliAlignObjAngles*>(array->At(i));
+         AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
          if (!a->ApplyToGeometry()) {
            AliWarning(Form("Failed to apply alignment to %s", 
-                           a->GetVolPath()));
+                           a->GetSymName()));
          }
        }
       }
     }
   }
 
-  
+  fEventCount = 0;
   fIsInit = kTRUE;
   return fIsInit;
 }
@@ -262,61 +271,87 @@ AliFMDInput::Begin(Int_t event)
     AliError("Not initialized");
     return fIsInit;
   }
+
   // Get the event 
   if (fLoader->GetEvent(event)) return kFALSE;
-  AliInfo(Form("Now in event %d/%d", event, NEvents()));
+  AliInfo(Form("Now in event %8d/%8d", event, NEvents()));
 
   // Possibly load global kinematics information 
-  if (TESTBIT(fTreeMask, kKinematics)) {
-    AliInfo("Getting kinematics");
+  if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
+    // AliInfo("Getting kinematics");
     if (fLoader->LoadKinematics()) return kFALSE;
     fStack = fLoader->Stack();
   }
+
   // Possibly load FMD Hit information 
-  if (TESTBIT(fTreeMask, kHits)) {
-    AliInfo("Getting FMD hits");
-    if (fFMDLoader->LoadHits()) return kFALSE;
+  if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
+    // AliInfo("Getting FMD hits");
+    if (!fFMDLoader || fFMDLoader->LoadHits()) return kFALSE;
     fTreeH = fFMDLoader->TreeH();
     if (!fArrayH) fArrayH = fFMD->Hits(); 
   }
+
   // Possibly load FMD Digit information 
   if (TESTBIT(fTreeMask, kDigits)) {
-    AliInfo("Getting FMD digits");
-    if (fFMDLoader->LoadDigits()) return kFALSE;
+    // AliInfo("Getting FMD digits");
+    if (!fFMDLoader || fFMDLoader->LoadDigits()) return kFALSE;
     fTreeD = fFMDLoader->TreeD();
-    if (!fArrayD) fArrayD = fFMD->Digits();
+    if (fTreeD) {
+      if (!fArrayD) fArrayD = fFMD->Digits();
+    }
+    else {
+      fArrayD = 0;
+      AliWarning(Form("Failed to load FMD Digits"));
+    } 
   }
+
   // Possibly load FMD Sdigit information 
   if (TESTBIT(fTreeMask, kSDigits)) {
-    AliInfo("Getting FMD summable digits");
-    if (fFMDLoader->LoadSDigits()) return kFALSE;
+    // AliInfo("Getting FMD summable digits");
+    if (!fFMDLoader || fFMDLoader->LoadSDigits()) return kFALSE;
     fTreeS = fFMDLoader->TreeS();
     if (!fArrayS) fArrayS = fFMD->SDigits();
   }
+
   // Possibly load FMD RecPoints information 
   if (TESTBIT(fTreeMask, kRecPoints)) {
-    AliInfo("Getting FMD reconstructed points");
-    if (fFMDLoader->LoadRecPoints()) return kFALSE;
+    // AliInfo("Getting FMD reconstructed points");
+    if (!fFMDLoader || fFMDLoader->LoadRecPoints()) return kFALSE;
     fTreeR = fFMDLoader->TreeR();
     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
     fTreeR->SetBranchAddress("FMD",  &fArrayR);
-  }  // Possibly load FMD ESD information 
+  }  
+
+  // Possibly load FMD ESD information 
   if (TESTBIT(fTreeMask, kESD)) {
-    AliInfo("Getting FMD event summary data");
+    // AliInfo("Getting FMD event summary data");
     Int_t read = fChainE->GetEntry(event);
     if (read <= 0) return kFALSE;
-    fESD = fMainESD->GetFMDData();
+    fESD = fESDEvent->GetFMDData();
     if (!fESD) return kFALSE;
+#if 0
+    TFile* f = fChainE->GetFile();
+    if (f) {
+      TObject* o = f->GetStreamerInfoList()->FindObject("AliFMDMap");
+      if (o) {
+       TStreamerInfo* info = static_cast<TStreamerInfo*>(o);
+       std::cout << "AliFMDMap class version read is " 
+                 <<  info->GetClassVersion() << std::endl;
+      }
+    }
+    // fESD->CheckNeedUShort(fChainE->GetFile());
+#endif
   }
+
   // Possibly load FMD Digit information 
   if (TESTBIT(fTreeMask, kRaw)) {
-    AliInfo("Getting FMD raw data digits");
+    // AliInfo("Getting FMD raw data digits");
     if (!fReader->NextEvent()) return kFALSE;
     AliFMDRawReader r(fReader, 0);
     fArrayA->Clear();
     r.ReadAdcs(fArrayA);
   }
-  
+  fEventCount++;
   return kTRUE;
 }
 
@@ -335,6 +370,8 @@ AliFMDInput::Event()
   // 
   if (TESTBIT(fTreeMask, kHits)) 
     if (!ProcessHits()) return kFALSE; 
+  if (TESTBIT(fTreeMask, kTracks)) 
+    if (!ProcessTracks()) return kFALSE; 
   if (TESTBIT(fTreeMask, kDigits)) 
     if (!ProcessDigits()) return kFALSE;
   if (TESTBIT(fTreeMask, kSDigits)) 
@@ -344,7 +381,7 @@ AliFMDInput::Event()
   if (TESTBIT(fTreeMask, kRecPoints)) 
     if (!ProcessRecPoints()) return kFALSE;
   if (TESTBIT(fTreeMask, kESD))
-    if (!ProcessESD(fESD)) return kFALSE;
+    if (!ProcessESDs()) return kFALSE;
   
   return kTRUE;
 }
@@ -359,19 +396,23 @@ AliFMDInput::ProcessHits()
     AliError("No hit tree defined");
     return kFALSE;
   }
+  if (!fArrayH) {
+    AliError("No hit array defined");
+    return kFALSE;
+  }
+
   Int_t nTracks = fTreeH->GetEntries();
   for (Int_t i = 0; i < nTracks; i++) {
     Int_t hitRead  = fTreeH->GetEntry(i);
     if (hitRead <= 0) continue;
-    if (!fArrayH) {
-      AliError("No hit array defined");
-      return kFALSE;
-    }
+
     Int_t nHit = fArrayH->GetEntries();
     if (nHit <= 0) continue;
+  
     for (Int_t j = 0; j < nHit; j++) {
       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
       if (!hit) continue;
+
       TParticle* track = 0;
       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
        Int_t trackno = hit->Track();
@@ -383,6 +424,50 @@ AliFMDInput::ProcessHits()
   return kTRUE;
 }
 
+//____________________________________________________________________
+Bool_t 
+AliFMDInput::ProcessTracks()
+{
+  // Read the hit tree, and pass each hit to the member function
+  // ProcessTrack.
+  if (!fStack) {
+    AliError("No track tree defined");
+    return kFALSE;
+  }
+  if (!fTreeH) {
+    AliError("No hit tree defined");
+    return kFALSE;
+  }
+  if (!fArrayH) {
+    AliError("No hit array defined");
+    return kFALSE;
+  }
+
+  // Int_t nTracks = fStack->GetNtrack();
+  Int_t nTracks = fTreeH->GetEntries();
+  for (Int_t i = 0; i < nTracks; i++) {
+    Int_t      trackno = nTracks - i - 1;
+    TParticle* track   = fStack->Particle(trackno);
+    if (!track) continue;
+
+    // Get the hits for this track. 
+    Int_t hitRead  = fTreeH->GetEntry(i);
+    Int_t nHit     = fArrayH->GetEntries();
+    if (nHit == 0 || hitRead <= 0) { 
+      // Let user code see the track, even if there's no hits. 
+      if (!ProcessTrack(trackno, track, 0)) return kFALSE;
+      continue;
+    }
+
+    // Loop over the hits corresponding to this track.
+    for (Int_t j = 0; j < nHit; j++) {
+      AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
+      if (!ProcessTrack(trackno, track, hit)) return kFALSE;
+    }   
+  }
+  return kTRUE;
+}
+
 //____________________________________________________________________
 Bool_t 
 AliFMDInput::ProcessDigits()
@@ -461,6 +546,31 @@ AliFMDInput::ProcessRecPoints()
   return kTRUE;
 }
 
+//____________________________________________________________________
+Bool_t 
+AliFMDInput::ProcessESDs()
+{
+  // Process event summary data
+  if (!fESD) return kFALSE;
+  for (UShort_t det = 1; det <= 3; det++) {
+    Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
+    for (Char_t* rng = rings; *rng != '\0'; rng++) {
+      UShort_t nsec = (*rng == 'I' ?  20 :  40);
+      UShort_t nstr = (*rng == 'I' ? 512 : 256);
+      for (UShort_t sec = 0; sec < nsec; sec++) {
+       for (UShort_t str = 0; str < nstr; str++) {
+         Float_t eta  = fESD->Eta(det,*rng,sec,str);
+         Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
+         if (!fESD->IsAngleCorrected()) 
+           mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
+         if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
+       }
+      }
+    }
+  }
+  return kTRUE;
+}
+
 //____________________________________________________________________
 Bool_t
 AliFMDInput::End()
@@ -476,13 +586,13 @@ AliFMDInput::End()
     return fIsInit;
   }
   // Possibly unload global kinematics information 
-  if (TESTBIT(fTreeMask, kKinematics)) {
+  if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
     fLoader->UnloadKinematics();
     // fTreeK = 0;
     fStack = 0;
   }
   // Possibly unload FMD Hit information 
-  if (TESTBIT(fTreeMask, kHits)) {
+  if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
     fFMDLoader->UnloadHits();
     fTreeH = 0;
   }
@@ -501,7 +611,7 @@ AliFMDInput::End()
     fFMDLoader->UnloadRecPoints();
     fTreeR = 0;
   }
-  AliInfo("Now out event");
+  // AliInfo("Now out event");
   return kTRUE;
 }
 
@@ -525,6 +635,31 @@ AliFMDInput::Run()
   return retval;
 }
 
+//__________________________________________________________________
+TArrayF 
+AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max) 
+{
+  // Service function to define a logarithmic axis. 
+  // Parameters: 
+  //   n    Number of bins 
+  //   min  Minimum of axis 
+  //   max  Maximum of axis 
+  TArrayF bins(n+1);
+  bins[0]      = min;
+  if (n <= 20) {
+    for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
+    return bins;
+  }
+  Float_t dp   = n / TMath::Log10(max / min);
+  Float_t pmin = TMath::Log10(min);
+  for (Int_t i = 1; i < n+1; i++) {
+    Float_t p = pmin + i / dp;
+    bins[i]   = TMath::Power(10, p);
+  }
+  return bins;
+}
+
+
 
 //____________________________________________________________________
 //