A lot of changes after detector review:
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Jan 2009 16:07:40 +0000 (16:07 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Jan 2009 16:07:40 +0000 (16:07 +0000)
* AliFMDBaseDigit (and AliFMDDigit and AliFMDSDigit
  by extension) derive from silly AliDigit.
* Up to 3 track labels propegated from hits to
  digits and sdigits.
* Implemented AliFMD::Raw2SDigits (see also
  scripts/TestRaw2SDigits.C).
* Added survey data of FMD1 and 2.
* Added FMD2 survey markers to geomtry.
* Survey to alignment is almost done.  The code
  is in AliFMDSurveyToAlignObjs.  Note, FMD2
  survey can only give FMD2 z-offset, while FMD1
  survey can in priniple also give x,y offsets
  as well as rotations.  FMD2 z-offset is determined
  as a best-fit of a plane to the available data.
* FMD2 and 3 are still roughly 4mm off from design. I
  need to talk to ITS people about this.

36 files changed:
FMD/AliFMD.cxx
FMD/AliFMD.h
FMD/AliFMDBackgroundCorrection.cxx
FMD/AliFMDBaseDigit.cxx
FMD/AliFMDBaseDigit.h
FMD/AliFMDBaseDigitizer.cxx
FMD/AliFMDBaseDigitizer.h
FMD/AliFMDDigit.cxx
FMD/AliFMDDigit.h
FMD/AliFMDDigitizer.cxx
FMD/AliFMDGeometryBuilder.cxx
FMD/AliFMDHitDigitizer.cxx
FMD/AliFMDInput.cxx
FMD/AliFMDInput.h
FMD/AliFMDRawReader.cxx
FMD/AliFMDRawReader.h
FMD/AliFMDReconstructor.cxx
FMD/AliFMDReconstructor.h
FMD/AliFMDSDigit.cxx
FMD/AliFMDSDigit.h
FMD/AliFMDStripIndex.h
FMD/AliFMDSurveyToAlignObjs.cxx
FMD/AliFMDSurveyToAlignObjs.h
FMD/Config.C
FMD/Reconstruct.C
FMD/Simulate.C
FMD/Survey_943928_FMD.txt [new file with mode: 0644]
FMD/Survey_976326_FMD.txt [new file with mode: 0644]
FMD/scripts/DrawBothDigits.C [new file with mode: 0644]
FMD/scripts/DrawDigits.C
FMD/scripts/DrawESD.C
FMD/scripts/DrawHits.C
FMD/scripts/DrawHitsDigits.C
FMD/scripts/DrawHitsSDigits.C [new file with mode: 0644]
FMD/scripts/DrawSDigits.C
FMD/scripts/TestRaw2SDigits.C [new file with mode: 0644]

index 0f404539c22b25408ab6f8ccc42b93dd6167cbb9..8e778889ebac0a9cf2009c34aba46dcec3d9884f 100644 (file)
 //#endif
 // #include "AliFMDGeometryBuilder.h"
 #include "AliFMDRawWriter.h"   // ALIFMDRAWWRITER_H
+#include "AliFMDRawReader.h"   // ALIFMDRAWREADER_H
 #include "AliTrackReference.h" 
 #include "AliFMDStripIndex.h"
+#include "AliFMDParameters.h"
+#include "AliFMDReconstructor.h"
 
 //____________________________________________________________________
 ClassImp(AliFMD)
@@ -696,14 +699,15 @@ AliFMD::AddDigit(Int_t* digits, Int_t*)
 
 //____________________________________________________________________
 void 
-AliFMD::AddDigitByFields(UShort_t detector, 
-                        Char_t   ring, 
-                        UShort_t sector, 
-                        UShort_t strip, 
-                        UShort_t count1, 
-                        Short_t  count2,
-                        Short_t  count3, 
-                        Short_t  count4)
+AliFMD::AddDigitByFields(UShort_t       detector, 
+                        Char_t         ring, 
+                        UShort_t       sector, 
+                        UShort_t       strip, 
+                        UShort_t       count1, 
+                        Short_t        count2,
+                        Short_t        count3, 
+                        Short_t        count4,
+                        const TArrayI& refs)
 {
   // add a real digit - as coming from data
   // 
@@ -719,11 +723,13 @@ AliFMD::AddDigitByFields(UShort_t detector,
   TClonesArray& a = *(DigitsArray());
   
   new (a[fNdigits++]) 
-    AliFMDDigit(detector, ring, sector, strip, count1, count2, count3, count4);
-  AliFMDDebug(15, ("Adding digit # %5d/%5d for FMD%d%c[%2d,%3d]=(%d,%d,%d,%d)",
+    AliFMDDigit(detector, ring, sector, strip, 
+               count1, count2, count3, count4, refs);
+  AliFMDDebug(15, ("Adding digit # %5d/%5d for FMD%d%c[%2d,%3d]"
+                  "=(%d,%d,%d,%d) with %d tracks",
                   fNdigits-1, a.GetEntriesFast(),
                   detector, ring, sector, strip, 
-                  count1, count2, count3, count4));
+                  count1, count2, count3, count4, refs.fN));
   
 }
 
@@ -907,6 +913,45 @@ AliFMD::Digits2Raw()
   writer.Exec();
 }
 
+//====================================================================
+//
+// Raw data reading 
+//
+//__________________________________________________________________
+Bool_t
+AliFMD::Raw2SDigits(AliRawReader* reader) 
+{
+  // Turn digits into raw data. 
+  // 
+  // This uses the class AliFMDRawWriter to do the job.   Please refer
+  // to that class for more information. 
+  AliFMDParameters::Instance()->Init();
+  MakeTree("S");
+  MakeBranch("S");
+  
+  TClonesArray*       sdigits = SDigits();
+  AliFMDReconstructor rec;
+  
+  // The two boolean arguments
+  //   Make sdigits instead of digits 
+  //   Subtract the pedestal off the signal
+  rec.Digitize(reader, sdigits);
+  // 
+  // Bool_t ret = fmdReader.ReadAdcs(sdigits, kTRUE, kTRUE);
+  // sdigits->ls();
+  UShort_t ns = sdigits->GetEntriesFast();
+  for (UShort_t i = 0; i < ns; i++) 
+    sdigits->At(i)->Print("pl");
+  
+  AliFMDDebug(1, ("Got a total of %d SDigits", ns));
+
+  fLoader->TreeS()->Fill();
+  ResetSDigits();
+  fLoader->WriteSDigits("OVERWRITE");
+
+  return kTRUE;
+}
+
 
 //====================================================================
 //
index b1c16407adb49a1e3365f6634241c34d21ac2be3..0e3c03728195ef514939017342e4f50614cc913f 100644 (file)
@@ -433,14 +433,15 @@ public:
       @param count1    ADC count (a 10-bit word)
       @param count2    ADC count (a 10-bit word), or -1 if not used
       @param count3    ADC count (a 10-bit word), or -1 if not used */
-  virtual        void   AddDigitByFields(UShort_t detector=0, 
-                                        Char_t   ring='\0', 
-                                        UShort_t sector=0, 
-                                        UShort_t strip=0, 
-                                        UShort_t count1=0, 
-                                        Short_t  count2=-1, 
-                                        Short_t  count3=-1, 
-                                        Short_t  count4=-1);
+  virtual        void   AddDigitByFields(UShort_t       detector=0, 
+                                        Char_t         ring='\0', 
+                                        UShort_t       sector=0, 
+                                        UShort_t       strip=0, 
+                                        UShort_t       count1=0, 
+                                        Short_t        count2=-1, 
+                                        Short_t        count3=-1, 
+                                        Short_t        count4=-1, 
+                                        const TArrayI& refs=TArrayI());
   /** Add a digit to the Digit tree 
       @param digits
       - digits[0]  [UShort_t] Detector #
@@ -495,16 +496,30 @@ public:
 
   /** @{ */
   /** @name Raw data */
-  /** Turn digits into raw data. This uses the class AliFMDRawWriter
-      to do the job.   Please refer to that class for more
-      information. */
-  virtual        void   Digits2Raw();
+  /** 
+   * Turn digits into raw data. This uses the class AliFMDRawWriter to
+   * do the job.  Please refer to that class for more information.
+   */
+  virtual void   Digits2Raw();
+  /** 
+   * Convert raw data to sdigits
+   * 
+   * @param reader Raw reader
+   * 
+   * @return @c true on success
+   */
+  virtual Bool_t Raw2SDigits(AliRawReader* reader);
   /** @}*/
 
   /** @{ */
-  /** @name Utility */
-  /** Browse this object 
-      @param b Browser to show this object in */
+  /** 
+   * @name Utility 
+   */
+  /** 
+   * Browse this object 
+   *
+   * @param b Browser to show this object in 
+   */
   void   Browse(TBrowser* b);
   /** @}*/
 protected:
index 94daa71a60debaf8c2246dc1d955d63af84ccc06..17a5770cf60a00fa4a2aed8c66c9127e0ede540a 100644 (file)
@@ -303,37 +303,48 @@ void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
 }
 
 //_____________________________________________________________________
-Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit( AliFMDHit* h, TParticle* p) {
+Bool_t 
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h, 
+                                                     TParticle* /*p*/) 
+{
   
   if(!h)
     return kTRUE;
-  Bool_t retval = ProcessEvent(h->Detector(),h->Ring(),h->Sector(),h->Strip(),h->Track(),h->Q());
+  Bool_t retval = ProcessEvent(h->Detector(),
+                              h->Ring(),
+                              h->Sector()
+                              ,h->Strip(),
+                              h->Track(),
+                              h->Q());
   
   return retval;
 }
 //_____________________________________________________________________
-Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef( AliTrackReference* tr, TParticle* p) {
-  
+Bool_t 
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr, 
+                                                          TParticle* p) 
+{
   if(!tr)
     return kTRUE;
   UShort_t det,sec,strip;
   Char_t   ring;
   AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip);
-  Int_t    nTrack = tr->GetTrack();
-  TDatabasePDG* pdgDB = TDatabasePDG::Instance();
+  Int_t         nTrack  = tr->GetTrack();
+  TDatabasePDG* pdgDB   = TDatabasePDG::Instance();
   TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
-  Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
-  Bool_t retval  = ProcessEvent(det,ring,sec,strip,nTrack,charge);
+  Float_t       charge  = (pdgPart ? pdgPart->Charge() : 0);
+  Bool_t        retval  = ProcessEvent(det,ring,sec,strip,nTrack,charge);
   return retval;
   
 }
 //_____________________________________________________________________
-Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
-                                                              Char_t   ring, 
-                                                              UShort_t sec, 
-                                                              UShort_t strip,
-                                                              Int_t    nTrack,
-                                                              Float_t  charge)
+Bool_t 
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
+                                                       Char_t   ring, 
+                                                       UShort_t sec, 
+                                                       UShort_t strip,
+                                                       Int_t    nTrack,
+                                                       Float_t  charge)
 {
   
   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
@@ -343,12 +354,14 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
   if(TMath::Abs(vertex.At(2)) > fZvtxCut) 
     return kTRUE;
   
-  Double_t delta = 2*fZvtxCut/fNvtxBins;
+  Double_t delta           = 2 * fZvtxCut / fNvtxBins;
   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
-  Int_t vertexBin = (Int_t)vertexBinDouble;
+  Int_t    vertexBin       = Int_t(vertexBinDouble);
   
-  if(charge !=  0 && (nTrack != fPrevTrack || det != fPrevDetector || ring != fPrevRing)) {
-    
+  if(charge !=  0 && 
+     (nTrack != fPrevTrack || 
+      det != fPrevDetector || 
+      ring != fPrevRing)) {
     Double_t x,y,z;
     AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
     fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
@@ -359,11 +372,11 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
     TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
     TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
     
-    Float_t phi   = TMath::ATan2(y,x);
-    if(phi<0)
-      phi = phi+2*TMath::Pi();
-    Float_t theta = TMath::ATan2(TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2)),z+vertex.At(2));
-    Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
+    Float_t   phi   = TMath::ATan2(y,x);
+    if(phi<0) phi   = phi+2*TMath::Pi();
+    Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
+    Float_t   theta = TMath::ATan2(r,z+vertex.At(2));
+    Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
     hHits->Fill(eta,phi);
     fHits++;
     fPrevDetector = det;
@@ -376,8 +389,8 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
 }
 
 //_____________________________________________________________________
-Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() {
-  
+Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() 
+{
   fPrimaryArray.SetOwner();
   fPrimaryArray.SetName("FMD_primary");
   
@@ -401,49 +414,40 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() {
   fHitArray.SetOwner();
   fHitArray.SetName("FMD_hits");
    
-  for(Int_t det =1; det<=3;det++)
-    {
-      TObjArray* detArrayHits = new TObjArray();
-      detArrayHits->SetName(Form("FMD%d",det));
-      fHitArray.AddAtAndExpand(detArrayHits,det);
-      Int_t nRings = (det==1 ? 1 : 2);
-      for(Int_t ring = 0;ring<nRings;ring++)
-       {
-         Int_t nSec = (ring == 1 ? 40 : 20);
-         Char_t ringChar = (ring == 0 ? 'I' : 'O');
-         TObjArray* vtxArrayHits = new TObjArray();
-         vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
-         detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
-         for(Int_t v=0; v<fNvtxBins;v++)
-           {
-             
-             TH2F* hHits          = new TH2F(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v),
-                                             Form("hHits_FMD%d%c_vtx%d",det,ringChar,v),
-                                             fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
-             hHits->Sumw2();
-             vtxArrayHits->AddAtAndExpand(hHits,v);
-                     
-           } 
-       }
-      
+  for(Int_t det =1; det<=3;det++) {
+    TObjArray* detArrayHits = new TObjArray();
+    detArrayHits->SetName(Form("FMD%d",det));
+    fHitArray.AddAtAndExpand(detArrayHits,det);
+    Int_t nRings = (det==1 ? 1 : 2);
+    for(Int_t ring = 0;ring<nRings;ring++) {
+      Int_t nSec = (ring == 1 ? 40 : 20);
+      Char_t ringChar = (ring == 0 ? 'I' : 'O');
+      TObjArray* vtxArrayHits = new TObjArray();
+      vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
+      detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
+      for(Int_t v=0; v<fNvtxBins;v++) {
+       TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
+                              Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
+                              fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
+       hHits->Sumw2();
+       vtxArrayHits->AddAtAndExpand(hHits,v);
+       
+      } 
     }
+  }
 
-
-  
   AliFMDInput::Init();
   
-  
   return kTRUE;
 }
 //_____________________________________________________________________
 
-Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) {
+Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) 
+{
 
-  Bool_t retVal = AliFMDInput::Begin(event); 
-  
-  AliStack* partStack=fLoader->Stack();
-  
-  Int_t nTracks=partStack->GetNtrack();
+  Bool_t             retVal    = AliFMDInput::Begin(event); 
+  AliStack*          partStack = fLoader->Stack();
+  Int_t              nTracks   = partStack->GetNtrack();
   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
   TArrayF vertex;
   genHeader->PrimaryVertex(vertex);
@@ -451,37 +455,35 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) {
   if(TMath::Abs(vertex.At(2)) > fZvtxCut) 
     return kTRUE;
   
-  Double_t delta = 2*fZvtxCut/fNvtxBins;
+  Double_t delta           = 2*fZvtxCut/fNvtxBins;
   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
-  Int_t vertexBin = (Int_t)vertexBinDouble;
-  
-  for(Int_t j=0;j<nTracks;j++)
-       {
-         TParticle* p  = partStack->Particle(j);
-         TDatabasePDG* pdgDB = TDatabasePDG::Instance();
-         TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
-         Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
-         Float_t phi = TMath::ATan2(p->Py(),p->Px());
-         
-         if(phi<0)
-           phi = phi+2*TMath::Pi();
-         if(p->Theta() == 0) continue;
-             Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*p->Theta()));
-         
-             Bool_t primary = partStack->IsPhysicalPrimary(j);
-             //(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01);
-             if(primary && charge !=0) {
-               
-               fPrim++;
-               TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
-               TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
-               
-               TH2F* hPrimaryInner  = (TH2F*)innerArray->At(vertexBin);
-               TH2F* hPrimaryOuter  = (TH2F*)outerArray->At(vertexBin);
-               hPrimaryInner->Fill(eta,phi);
-               hPrimaryOuter->Fill(eta,phi);
-             }
-       }
+  Int_t    vertexBin       = (Int_t)vertexBinDouble;
+  
+  for(Int_t j=0;j<nTracks;j++) {
+    TParticle* p  = partStack->Particle(j);
+    TDatabasePDG* pdgDB   = TDatabasePDG::Instance();
+    TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
+    Float_t       charge  = (pdgPart ? pdgPart->Charge() : 0);
+    Float_t       phi     = TMath::ATan2(p->Py(),p->Px());
+    
+    if(phi<0) phi = phi+2*TMath::Pi();
+    if(p->Theta() == 0) continue;
+    Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*p->Theta()));
+    
+    Bool_t primary = partStack->IsPhysicalPrimary(j);
+    //(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01);
+    if(primary && charge !=0) {
+      
+      fPrim++;
+      TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
+      TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
+      
+      TH2F* hPrimaryInner  = (TH2F*)innerArray->At(vertexBin);
+      TH2F* hPrimaryOuter  = (TH2F*)outerArray->At(vertexBin);
+      hPrimaryInner->Fill(eta,phi);
+      hPrimaryOuter->Fill(eta,phi);
+    }
+  }
   
   return retVal;
 }
index 5ac0b767e61d073d6211b23398a8b39e2892b988..41d69200d486a0c6d952585a046cc5939c730bd3 100644 (file)
@@ -64,6 +64,7 @@
 //////////////////////////////////////////////////////////////////////
 
 #include "AliFMDBaseDigit.h"   // ALIFMDDIGIT_H
+#include "AliFMDStripIndex.h"
 #include "Riostream.h"         // ROOT_Riostream
 // #include <TString.h>
 // #include <AliLog.h>
@@ -89,7 +90,8 @@ AliFMDBaseDigit::AliFMDBaseDigit(UShort_t detector,
                                 Char_t   ring, 
                                 UShort_t sector, 
                                 UShort_t strip)
-  : fDetector(detector), 
+  : AliDigit(), 
+    fDetector(detector), 
     fRing(ring), 
     fSector(sector), 
     fStrip(strip),
@@ -106,12 +108,37 @@ AliFMDBaseDigit::AliFMDBaseDigit(UShort_t detector,
   //    strip    Strip # (For inner/outer rings: 0-511/0-255)
 }
 
+//____________________________________________________________________
+AliFMDBaseDigit::AliFMDBaseDigit(Int_t*   tracks, 
+                                UShort_t detector, 
+                                Char_t   ring, 
+                                UShort_t sector, 
+                                UShort_t strip)
+  : AliDigit(tracks), 
+    fDetector(detector), 
+    fRing(ring), 
+    fSector(sector), 
+    fStrip(strip),
+    fName("")
+{
+  //
+  // Creates a base data digit object
+  //
+  // Parameters 
+  //
+  //    tracks    Array of 3 track labels
+  //    detector  Detector # (1, 2, or 3)                      
+  //    ring     Ring ID ('I' or 'O')
+  //    sector   Sector # (For inner/outer rings: 0-19/0-39)
+  //    strip    Strip # (For inner/outer rings: 0-511/0-255)
+}
+
 //____________________________________________________________________
 void
 AliFMDBaseDigit::Print(Option_t* /* option*/) const 
 {
   // Print digit to standard out 
-  cout << ClassName() << ": FMD" << GetName() << flush;
+  cout << ClassName() << ": " << GetName() << flush;
 }
 
 //____________________________________________________________________
@@ -132,9 +159,13 @@ ULong_t
 AliFMDBaseDigit::Hash() const
 {
   // Calculate a hash value based on the detector coordinates. 
+#if 1  
+  return AliFMDStripIndex::Pack(fDetector, fRing, fSector, fStrip);
+#else
   size_t ringi = (fRing == 'I' ||  fRing == 'i' ? 0 : 1);
   return fStrip + fMaxStrips * 
     (fSector + fMaxSectors * (ringi + fMaxRings * (fDetector - 1)));
+#endif
 }
 
 
@@ -164,6 +195,28 @@ AliFMDBaseDigit::Compare(const TObject* o) const
   return 1;
 }
 
+//____________________________________________________________________
+void
+AliFMDBaseDigit::AddTrack(Int_t track)
+{
+  if      (fTracks[0] == -1) fTracks[0] = track;
+  else if (fTracks[1] == -1) fTracks[1] = track;
+  else if (fTracks[2] == -1) fTracks[2] = track;
+  else 
+    AliWarning(Form("While adding track label to %s for %s: "
+                   "All 3 track labels used, cannot add reference to track %d",
+                   ClassName(), GetName(), track));
+}
+
+//____________________________________________________________________
+UShort_t
+AliFMDBaseDigit::GetNTrack() const
+{
+  for (Int_t i = 3; i > 0; i--) 
+    if (fTracks[i-1] != -1) return i;
+  return 0;
+}
+
 
 //____________________________________________________________________
 //
index fe9b7e70d2ddf0496d68a073167e5b59c847de07..a6d88a7b440111652f08447b0606ba87d92a8b6f 100644 (file)
 //  AliFMDDigit     - Normal (smeared) digit             
 //  AliFMDSDigit    - Summable (non-smeared) digit             
 //
-#ifndef ROOT_TObject
-# include <TObject.h>
+#ifndef ALIDIGIT_H
+# include <AliDigit.h>
 #endif
 #ifndef ROOT_TString
 # include <TString.h>
 #endif
 
 //____________________________________________________________________
-/** @class AliFMDBaseDigit AliFMDDigit.h <FMD/AliFMDDigit.h>
-    @brief base class for digits 
-    @ingroup FMD_base
+/** 
+ * @class AliFMDBaseDigit AliFMDDigit.h <FMD/AliFMDDigit.h>
+ * 
+ * @brief base class for digits 
+ * 
+ * @ingroup FMD_base
  */
-class AliFMDBaseDigit : public TObjec
+class AliFMDBaseDigit : public AliDigi
 {
 public: 
-  /** CTOR */
+  /** 
+   * CTOR 
+   */
   AliFMDBaseDigit();
-  /** Constrctor 
-      @param detector Detector 
-      @param ring     Ring
-      @param sector   Sector
-      @param strip    Strip */
+  /** 
+   * Constrctor 
+   * 
+   * @param detector Detector 
+   * @param ring     Ring
+   * @param sector   Sector
+   * @param strip    Strip 
+   */
   AliFMDBaseDigit(UShort_t detector, 
                  Char_t   ring='\0', 
                  UShort_t sector=0, 
                  UShort_t strip=0);
-  /** DTOR */
+  /** 
+   * Constrctor 
+   *
+   * @param tracks   Array of 3 track indicies
+   * @param detector Detector 
+   * @param ring     Ring
+   * @param sector   Sector
+   * @param strip    Strip 
+   */
+  AliFMDBaseDigit(Int_t*   tracks, 
+                 UShort_t detector, 
+                 Char_t   ring='\0', 
+                 UShort_t sector=0, 
+                 UShort_t strip=0);  
+  /** 
+   * DTOR 
+   */
   virtual ~AliFMDBaseDigit() {}
-  /** @return Detector # */
+  /** 
+   * 
+   * @return Detector # 
+   */
   UShort_t     Detector()         const { return fDetector; }
-  /** @return Ring ID */
+  /** 
+   * 
+   * @return Ring ID 
+   */
   Char_t       Ring()             const { return fRing;     }
-  /** @return sector # */
+  /** 
+   * 
+   * @return sector # 
+   */
   UShort_t     Sector()                   const { return fSector;   }
-  /** @return strip # */
+  /** 
+   * 
+   * @return strip # 
+   */
   UShort_t     Strip()            const { return fStrip;    }
-  /** Print information 
-      @param opt Not used */
+  /** 
+   * Print information 
+   *
+   * @param opt Not used 
+   */
   virtual void Print(Option_t* opt="") const;
-  /** @return Name */
+  /** 
+   * 
+   * @return Name 
+   */
   const char*  GetName() const;
-  /** @param rhs Other digit to compare to 
-      @return -1 if this is less than  @a rhs, 0 if the refer to the
-      same, and 1 if @a rhs is larger than this */
+  /** 
+   * @param rhs Other digit to compare to 
+   * 
+   * @return -1 if this is less than  @a rhs, 0 if the refer to the
+   * same, and 1 if @a rhs is larger than this 
+   */
   Int_t Compare(const TObject* o) const;
-  /** @return Always true */ 
+  /** 
+   * 
+   * @return Always true 
+   */ 
   Bool_t IsSortable() const { return kTRUE; }
+  
+  /** 
+   * Add a track referenc
+   * 
+   * @param trackno The track number
+   */  
+  void AddTrack(Int_t trackno);
+  
+  /** 
+   * Get the number of track references (max 3)
+   * 
+   * 
+   * @return Number of valid track references. 
+   */
+  UShort_t GetNTrack() const;
+  
+  /** 
+   * Set the count value 
+   * 
+   * @param s Sample number 
+   * @param c Counts 
+   */
+  virtual void SetCount(UShort_t s, Short_t c) = 0;
 protected:
+  /** 
+   * Calculate the hash value
+   * 
+   * 
+   * @return Hash value 
+   */  
   ULong_t  Hash() const;
   UShort_t fDetector;  // (Sub) Detector # (1,2, or 3)
   Char_t   fRing;      // Ring ID ('I' or 'O')
   UShort_t fSector;    // Sector # (phi division)
   UShort_t fStrip;     // Strip # (radial division)
   mutable TString  fName;      //! Name (cached, but not stored) 
-  ClassDef(AliFMDBaseDigit, 2) // Base class for FMD digits 
+  ClassDef(AliFMDBaseDigit, 3) // Base class for FMD digits 
 };
 
 #endif
index 780b27c3f1d48d4f0a0e4caa80c5d0c58b72ae0b..36bf25549b796ba25057b8af693d5b7dec665c04 100644 (file)
@@ -231,7 +231,7 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer()
          AliFMDMap::kMaxSectors, 
          AliFMDMap::kMaxStrips),
     fShapingTime(6),
-    fStoreTrackRefs(kFALSE)
+    fStoreTrackRefs(kTRUE)
 {
   AliFMDDebug(1, ("Constructed"));
   // Default ctor - don't use it
@@ -247,7 +247,7 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager)
          AliFMDMap::kMaxSectors, 
          AliFMDMap::kMaxStrips), 
     fShapingTime(6),
-    fStoreTrackRefs(kFALSE)
+    fStoreTrackRefs(kTRUE)
 {
   // Normal CTOR
   AliFMDDebug(1, ("Constructed"));
@@ -265,7 +265,7 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name,
          AliFMDMap::kMaxSectors, 
          AliFMDMap::kMaxStrips),
     fShapingTime(6),
-    fStoreTrackRefs(kFALSE)
+    fStoreTrackRefs(kTRUE)
 {
   // Normal CTOR
   AliFMDDebug(1, (" Constructed"));
@@ -314,7 +314,8 @@ AliFMDBaseDigitizer::AddContribution(UShort_t detector,
                                     UShort_t strip, 
                                     Float_t  edep, 
                                     Bool_t   isPrimary,
-                                    Int_t    trackno)
+                                    Int_t    nTrack,
+                                    Int_t*   tracknos)
 {
   // Add edep contribution from (detector,ring,sector,strip) to cache
   AliFMDParameters* param = AliFMDParameters::Instance();
@@ -341,15 +342,17 @@ AliFMDBaseDigitizer::AddContribution(UShort_t detector,
       
   // Sum energy deposition
   entry.fEdep += edep;
-  entry.fN    += 1;
-  if (isPrimary) entry.fNPrim += 1;
+  entry.fN    += nTrack;
+  if (isPrimary) entry.fNPrim += nTrack;
   if (fStoreTrackRefs) { 
+    Int_t oldN = entry.fLabels.fN;
     entry.fLabels.Set(entry.fN);
-    entry.fLabels[entry.fN-1] = trackno;
+    for (Int_t i = 0; i < nTrack; i++) 
+      entry.fLabels[oldN + i] = tracknos[i];
   }
-  AliFMDDebug(15, ("Adding contribution %f to FMD%d%c[%2d,%3d] (%f)", 
+  AliFMDDebug(15, ("Adding contribution %f to FMD%d%c[%2d,%3d] (%f) track %d", 
                   edep, detector, ring, sector, strip,
-                  entry.fEdep));
+                  entry.fEdep, (nTrack > 0 ? tracknos[0] : -1)));
   
 }
 
@@ -428,7 +431,7 @@ AliFMDBaseDigitizer::DigitizeHits() const
                   Short_t(counts[2]), Short_t(counts[3]), 
                   ntot, nprim, labels);
          AliFMDDebug(15, ("   Adding digit in FMD%d%c[%2d,%3d]=%d", 
-                         detector,ring,sector,strip,counts[0]));
+                          detector,ring,sector,strip,counts[0]));
 #if 0
          // This checks if the digit created will give the `right'
          // number of particles when reconstructed, using a naiive
@@ -550,12 +553,12 @@ AliFMDBaseDigitizer::AddDigit(UShort_t        detector,
                              Short_t         count4,
                              UShort_t        /* ntot */, 
                              UShort_t        /* nprim */,
-                             const TArrayI&  /* refs */) const
+                             const TArrayI&  refs) const
 {
   // Add a digit or summable digit
   
   fFMD->AddDigitByFields(detector, ring, sector, strip, 
-                        count1, count2, count3, count4);
+                        count1, count2, count3, count4, refs);
 }
 
 //____________________________________________________________________
index b5d77f8d9adfdb1147f34e013f8899dc8fe0eefc..1508eec3646d7f1597bb03113b84faeea23e6a5e 100644 (file)
@@ -230,7 +230,8 @@ protected:
                               UShort_t strip, 
                               Float_t  edep, 
                               Bool_t   isPrimary,
-                              Int_t    trackno);
+                              Int_t    nTrackno,
+                              Int_t*   tracknos);
   /** Add a digit to output */
   virtual void     AddDigit(UShort_t       detector, 
                            Char_t         ring,
index 2de6a9a039a81f75cc15a9af3c4d207e6ef90770..5941b9c0e5189031f8e4650012a90f795f65ba6c 100644 (file)
@@ -83,14 +83,15 @@ AliFMDDigit::AliFMDDigit()
 }
 
 //____________________________________________________________________
-AliFMDDigit::AliFMDDigit(UShort_t detector, 
-                        Char_t   ring, 
-                        UShort_t sector, 
-                        UShort_t strip, 
-                        UShort_t count1,
-                        Short_t  count2, 
-                        Short_t  count3,
-                        Short_t  count4)
+AliFMDDigit::AliFMDDigit(UShort_t       detector, 
+                        Char_t         ring, 
+                        UShort_t       sector, 
+                        UShort_t       strip, 
+                        UShort_t       count1,
+                        Short_t        count2, 
+                        Short_t        count3,
+                        Short_t        count4, 
+                        const TArrayI& refs)
   : AliFMDBaseDigit(detector, ring, sector, strip), 
     fCount1(count1),
     fCount2(count2),
@@ -109,6 +110,7 @@ AliFMDDigit::AliFMDDigit(UShort_t detector,
   //    count1    ADC count (a 10-bit word)
   //    count2    ADC count (a 10-bit word) -1 if not used
   //    count3    ADC count (a 10-bit word) -1 if not used
+  for (Int_t i = 0; i < refs.fN; i++) AddTrack(refs.fArray[i]);
 }
 
 //____________________________________________________________________
@@ -123,13 +125,23 @@ AliFMDDigit::GetTitle() const
 
 //____________________________________________________________________
 void
-AliFMDDigit::Print(Option_t* /* option*/) const 
+AliFMDDigit::Print(Option_t* option) const 
 {
   // Print digit to standard out 
   AliFMDBaseDigit::Print();
-  cout << "\t" 
-       << fCount1 << " (" << fCount2 << "," << fCount3 << "," << fCount4 
-       << ") = " << Counts() << endl;
+  std::cout << "\t" 
+           << std::setw(4) << fCount1 << " (" 
+           << std::setw(4) << fCount2 << "," 
+           << std::setw(4) << fCount3 << "," 
+           << std::setw(4) << fCount4 << ") = " 
+           << std::setw(4) << Counts() << std::flush;
+  TString opt(option);
+  if (opt.Contains("l", TString::kIgnoreCase)) {
+    std::cout << " ";
+    for (Int_t i = 0; i < 3; i++) 
+      std::cout << (i == 0 ? "" : ",") << std::setw(5) << fTracks[i];
+  }
+  std::cout << std::endl;
 }
 
 //____________________________________________________________________
index da905fbda17e4a9581f42ec7fb532bea5e830dcd..815f90e5d91cffda8ee98f22bdc3f02e6cb69962 100644 (file)
 #ifndef ALIFMDBASEDIGIT_H
 # include <AliFMDBaseDigit.h>
 #endif
+#ifndef ROOT_TArrayI
+# include <TArrayI.h>
+#endif
+
 
 //____________________________________________________________________
 /** @class AliFMDDigit AliFMDDigit.h <FMD/AliFMDDigit.h>
@@ -26,45 +30,79 @@ class AliFMDDigit : public AliFMDBaseDigit
 public:
   /** CTOR */
   AliFMDDigit();
-  /** Constrctor 
-      @param detector Detector 
-      @param ring     Ring
-      @param sector   Sector
-      @param strip    Strip 
-      @param count    ADC (first sample)
-      @param count2   ADC (second sample, or -1 if not used)
-      @param count3   ADC (third sample, or -1 if not used) */
-  AliFMDDigit(UShort_t detector, 
-             Char_t   ring='\0', 
-             UShort_t sector=0, 
-             UShort_t strip=0, 
-             UShort_t count=0, 
-             Short_t  count2=-1, 
-             Short_t  count3=-1, 
-             Short_t  count4=-1);
-  /** DTOR */
+  /** 
+   * Constrctor 
+   *
+   * @param detector Detector 
+   * @param ring     Ring
+   * @param sector   Sector
+   * @param strip    Strip 
+   * @param count    ADC (first sample)
+   * @param count2   ADC (second sample, or -1 if not used)
+   * @param count3   ADC (third sample, or -1 if not used) 
+   * @param refs     Track references
+   */
+  AliFMDDigit(UShort_t       detector, 
+             Char_t         ring='\0', 
+             UShort_t       sector=0, 
+             UShort_t       strip=0, 
+             UShort_t       count=0, 
+             Short_t        count2=-1, 
+             Short_t        count3=-1, 
+             Short_t        count4=-1, 
+             const TArrayI& refs=TArrayI());
+  /** 
+   * DTOR 
+   */
   virtual ~AliFMDDigit() {}
-  /** @param i # of sample to get 
-      @return sample # @a i */
+  /** 
+   * @param i # of sample to get 
+   * 
+   * @return sample # @a i 
+   */
   Int_t Count(UShort_t i=0) const;
-  /** @return ADC count (first sample) */
-  UShort_t Count1()                const { return fCount1;   }
-  /** @return ADC count (second sample, or -1 if not used) */
-  Short_t  Count2()                const { return fCount2;   }
-  /** @return ADC count (third sample, or -1 if not used) */
-  Short_t  Count3()                const { return fCount3;   }
-  /** @return ADC count (third sample, or -1 if not used) */
-  Short_t  Count4()                const { return fCount4;   }
-  /** @return Canonical ADC counts */
-  UShort_t Counts()                const;
-  /** Print info 
-      @param opt Not used */
+  /** 
+   * 
+   * @return ADC count (first sample) 
+   */
+  UShort_t Count1() const { return fCount1;   }
+  /** 
+   * 
+   * @return ADC count (second sample, or -1 if not used) 
+   */
+  Short_t  Count2() const { return fCount2;   }
+  /** 
+   * 
+   * @return ADC count (third sample, or -1 if not used) 
+   */
+  Short_t  Count3() const { return fCount3;   }
+  /** 
+   * 
+   * @return ADC count (third sample, or -1 if not used) 
+   */
+  Short_t  Count4() const { return fCount4;   }
+  /** 
+   * 
+   * @return Canonical ADC counts 
+   */
+  UShort_t Counts() const;
+  /** 
+   * Print info 
+   * 
+   * @param opt Not used 
+   */
   void     Print(Option_t* opt="") const;
-  /** @return Title */
+  /** 
+   * 
+   * @return Title 
+   */
   const char* GetTitle() const;
-  /** Set the count value 
-      @param s Sample number 
-      @param c Counts */
+  /** 
+   * Set the count value 
+   * 
+   * @param s Sample number 
+   * @param c Counts 
+   */
   void SetCount(UShort_t s, Short_t c);
 protected:
   UShort_t fCount1;     // Digital signal 
index 91e4ff8270ab638db8b2fbf527546ad237ce1e91..c17a265073221a93ba75a8e3f49887dbe8a28a26 100644 (file)
@@ -270,7 +270,7 @@ AliFMDDigitizer::Exec(Option_t*)
                    folderName.Data()));
       return;
     }
-    runLoader->LoadgAlice();
+    if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
     runLoader->GetAliRun();
   }
   if (!gAlice) { 
@@ -339,6 +339,7 @@ AliFMDDigitizer::Exec(Option_t*)
     fFMD->SetTreeAddress();
 
     // Sum contributions from the sdigits
+    AliFMDDebug(3, ("Will now sum contributions from SDigits"));
     SumContributions(sdigitsBranch);
 
     // Unload the sdigits
@@ -398,6 +399,16 @@ AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch)
       // Get the sdigit number `sdigit'
       AliFMDSDigit* fmdSDigit = 
        static_cast<AliFMDSDigit*>(fmdSDigits->UncheckedAt(sdigit));
+
+      AliFMDDebug(5, ("Adding contribution of %d tracks", 
+                     fmdSDigit->GetNTrack()));
+      AliFMDDebug(15, ("Contrib from FMD%d%c[%2d,%3d] (%s) from track %d", 
+                      fmdSDigit->Detector(),
+                      fmdSDigit->Ring(),
+                      fmdSDigit->Sector(),
+                      fmdSDigit->Strip(),
+                      fmdSDigit->GetName(), 
+                      fmdSDigit->GetTrack(0)));
       
       // Extract parameters 
       AddContribution(fmdSDigit->Detector(),
@@ -406,7 +417,8 @@ AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch)
                      fmdSDigit->Strip(),
                      fmdSDigit->Edep(), 
                      kTRUE,
-                     -1);
+                     fmdSDigit->GetNTrack(),
+                     fmdSDigit->GetTracks());
     }  // sdigit loop
   } // event loop
 
index ecb994f64f5ff9a6d34ed6a59dd7cd5d8825960c..85a39f187af386a74a447f564255c6da535c5083 100644 (file)
@@ -654,9 +654,12 @@ AliFMDGeometryBuilder::FMD1Geometry(AliFMD1* fmd1,
   AliFMDRing* r             = fmd1->GetInner();
   Double_t    z             = fmd1->GetInnerZ();  
   
+  // `Top' or `Outside' master volume
   TGeoVolume* fmd1TopVolume = new TGeoVolumeAssembly(Form(fgkFMDName, 
                                                          fmd1->GetId(), 'T'));
   fmd1TopVolume->SetTitle("FMD1 top half");
+
+  // `Bottom' or `Inside' master volume
   TGeoVolume* fmd1BotVolume = new TGeoVolumeAssembly(Form(fgkFMDName, 
                                                          fmd1->GetId(), 'B'));
   fmd1BotVolume->SetTitle("FMD1 bottom half");
@@ -728,11 +731,57 @@ AliFMDGeometryBuilder::FMD1Geometry(AliFMD1* fmd1,
   TGeoRotation* rot = new TGeoRotation("FMD1 rotatation");
   rot->RotateZ(90);
   TGeoMatrix* matrix = new TGeoCombiTrans("FMD1 trans", 0, 0, z, rot);
+
   AliFMDDebug(5, ("Placing volumes %s and %s in ALIC at z=%f", 
                   fmd1TopVolume->GetName(), fmd1BotVolume->GetName(), z));
   top->AddNode(fmd1TopVolume, fmd1->GetId(), matrix);
   top->AddNode(fmd1BotVolume, fmd1->GetId(), matrix);
+
+
+  // Survey points on V0A (screw holes for the FMD) 
+  const Double_t icb[] = { +12.700, -21.997, 324.670 };
+  const Double_t ict[] = { +12.700, +21.997, 324.670 };
+  const Double_t ocb[] = { -12.700, -21.997, 324.670 };
+  const Double_t oct[] = { -12.700, +21.997, 324.670 };
+
+  TGeoTube* surveyShape = new TGeoTube("FMD1_survey_marker", 
+                                       0, .2, .001);
+
+  TGeoMatrix* outMat = matrix;
+#if 0
+  if (gGeoManager->cd("/ALIC_1/F1MT_1")) 
+    outMat = gGeoManager->GetCurrentMatrix();
+  else 
+    AliWarning("Couldn't cd to /ALIC_1/F1MT_1");
+#endif
+
+  Double_t loct[3], locb[3];
+  outMat->MasterToLocal(oct, loct);
+  outMat->MasterToLocal(ocb, locb);
+  TGeoVolume* vOct = new TGeoVolume("V0L_OCT", surveyShape, fPlastic);
+  TGeoVolume* vOcb = new TGeoVolume("V0L_OCB", surveyShape, fPlastic);
+  
+  fmd1TopVolume->AddNode(vOct, 1, new TGeoTranslation(loct[0],loct[1],loct[2]));
+  fmd1TopVolume->AddNode(vOcb, 1, new TGeoTranslation(locb[0],locb[1],locb[2]));
+    
+  
+  TGeoMatrix* inMat = matrix;
+#if 0
+  if (gGeoManager->cd("/ALIC_1/F1MT_1")) 
+    inMat = gGeoManager->GetCurrentMatrix();
+  else 
+    AliWarning("Couldn't cd to /ALIC_1/F1MT_1");
+#endif
+
+  Double_t lict[3], licb[3];
+  inMat->MasterToLocal(ict, lict);
+  inMat->MasterToLocal(icb, licb);
+  TGeoVolume* vIct = new TGeoVolume("V0L_ICT", surveyShape, fPlastic);
+  TGeoVolume* vIcb = new TGeoVolume("V0L_ICB", surveyShape, fPlastic);
   
+  fmd1BotVolume->AddNode(vIct, 1, new TGeoTranslation(lict[0],lict[1],lict[2]));
+  fmd1BotVolume->AddNode(vIcb, 1, new TGeoTranslation(licb[0],licb[1],licb[2]));
+
   return 0;
 }
 
index 20714d60f4fc715291604f5f1b950b4e1edb5ecf..50132d81b4e82fdcce779f954762f58adf242838 100644 (file)
@@ -430,7 +430,8 @@ AliFMDHitDigitizer::SumContributions(TBranch* hitsBranch)
                      fmdHit->Strip(),
                      fmdHit->Edep(), 
                      isPrimary, 
-                     trackno);
+                     1, 
+                     &trackno);
     }  // hit loop
   } // track loop
   AliFMDDebug(5, ("Size of cache: %d bytes, read %d bytes", 
@@ -469,6 +470,9 @@ AliFMDHitDigitizer::AddDigit(UShort_t        detector,
 {
   // Add a digit or summable digit
   if (fOutput == kDigits) { 
+    AliFMDDebug(15,("Adding digit for FMD%d%c[%2d,%3d] = (%x,%x,%x,%x)",
+                   detector, ring, sector, strip, 
+                   count1, count2, count3, count4));
     AliFMDBaseDigitizer::AddDigit(detector, ring, sector, strip, 0,
                                  count1, count2, count3, count4, 0, 0, refs);
     return;
@@ -484,9 +488,10 @@ AliFMDHitDigitizer::AddDigit(UShort_t        detector,
                    detector, ring, sector, strip));
     return;
   }
-  AliFMDDebug(15, ("Adding digit for FMD%d%c[%2d,%3d] = (%x,%x,%x,%x) [%d/%d]",
+  AliFMDDebug(15, ("Adding sdigit for FMD%d%c[%2d,%3d] = "
+                  "(%x,%x,%x,%x) [%d/%d] %d",
                   detector, ring, sector, strip, 
-                  count1, count2, count3, count4, nprim, ntotal));
+                  count1, count2, count3, count4, nprim, ntotal, refs.fN));
   fFMD->AddSDigitByFields(detector, ring, sector, strip, edep,
                          count1, count2, count3, count4, 
                          ntotal, nprim, refs);
index a819d920afb0d2fe09dbd6c2e734f48cfeefd429..a8056824c602dc61db4f8b120685fd8a9d6ea816 100644 (file)
@@ -93,6 +93,7 @@ AliFMDInput::AliFMDInput()
     fChainE(0),
     fArrayE(0),
     fArrayH(0),
+    fArrayTR(0), 
     fArrayD(0),
     fArrayS(0), 
     fArrayR(0), 
@@ -136,6 +137,7 @@ AliFMDInput::AliFMDInput(const char* gAliceFile)
     fChainE(0),
     fArrayE(0),
     fArrayH(0),
+    fArrayTR(0), 
     fArrayD(0),
     fArrayS(0), 
     fArrayR(0), 
@@ -160,6 +162,7 @@ AliFMDInput::NEvents() const
 {
   // Get number of events
   if (TESTBIT(fTreeMask, kRaw)) return fReader->GetNumberOfEvents();
+  if (fChainE) return fChainE->GetEntriesFast();
   if (fTreeE) return fTreeE->GetEntries();
   return -1;
 }
@@ -175,38 +178,71 @@ AliFMDInput::Init()
     AliWarning("Already initialized");
     return fIsInit;
   }
+  Info("Init","Initialising w/mask 0x%04x\n"
+       "\tHits:          %d\n"
+       "\tKinematics:    %d\n"
+       "\tDigits:        %d\n"
+       "\tSDigits:       %d\n"
+       "\tHeader:        %d\n"
+       "\tRecPoints:     %d\n"
+       "\tESD:           %d\n"
+       "\tRaw:           %d\n"
+       "\tGeometry:      %d\n"
+       "\tTracks:        %d\n"
+       "\tTracksRefs:    %d",
+       fTreeMask,
+       TESTBIT(fTreeMask, kHits),
+       TESTBIT(fTreeMask, kKinematics),
+       TESTBIT(fTreeMask, kDigits),
+       TESTBIT(fTreeMask, kSDigits),
+       TESTBIT(fTreeMask, kHeader),
+       TESTBIT(fTreeMask, kRecPoints),
+       TESTBIT(fTreeMask, kESD),
+       TESTBIT(fTreeMask, kRaw),
+       TESTBIT(fTreeMask, kGeometry),
+       TESTBIT(fTreeMask, kTracks),
+       TESTBIT(fTreeMask, kTrackRefs));
   // Get the run 
-  if (!TESTBIT(fTreeMask, kRaw) && 
-      !TESTBIT(fTreeMask, kESD)) { 
-    if (fGAliceFile.IsNull()) fGAliceFile = "galice.root";
-    // Get the loader
-    fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
-    if (!fLoader) {
-      AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
-      return kFALSE;
+  if (TESTBIT(fTreeMask, kDigits)     ||
+      TESTBIT(fTreeMask, kSDigits)    || 
+      TESTBIT(fTreeMask, kKinematics) || 
+      TESTBIT(fTreeMask, kTrackRefs)  || 
+      TESTBIT(fTreeMask, kTracks)     || 
+      TESTBIT(fTreeMask, kHeader)) {
+    if (!gSystem->FindFile(".:/", fGAliceFile)) {
+      AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data()));
     }
-  
-    if  (fLoader->LoadgAlice()) return kFALSE;
-    fRun = fLoader->GetAliRun();
-  
-    // Get the FMD 
-    fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
-    if (!fFMD) {
-      AliError("Failed to get detector FMD from loader");
-      return kFALSE;
-    }
-  
-    // Get the FMD loader
-    fFMDLoader = fLoader->GetLoader("FMDLoader");
-    if (!fFMDLoader) {
-      AliError("Failed to get detector FMD loader from loader");
-      return kFALSE;
-    }
-    if (fLoader->LoadHeader()) { 
-      AliError("Failed to get event header information from loader");
-      return kFALSE;
+    else {
+      fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
+      if (!fLoader) {
+       AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
+       return kFALSE;
+      }
+      AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data()));
+
+      if  (fLoader->LoadgAlice()) return kFALSE;
+      
+      fRun = fLoader->GetAliRun();
+      
+      // Get the FMD 
+      fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
+      if (!fFMD) {
+       AliError("Failed to get detector FMD from loader");
+       return kFALSE;
+      }
+      
+      // Get the FMD loader
+      fFMDLoader = fLoader->GetLoader("FMDLoader");
+      if (!fFMDLoader) {
+       AliError("Failed to get detector FMD loader from loader");
+       return kFALSE;
+      }
+      if (fLoader->LoadHeader()) { 
+       AliError("Failed to get event header information from loader");
+       return kFALSE;
+      }
+      fTreeE = fLoader->TreeE();
     }
-    fTreeE = fLoader->TreeE();
   }
 
   // Optionally, get the ESD files
@@ -361,7 +397,10 @@ AliFMDInput::Begin(Int_t event)
   // Possibly load FMD Sdigit information 
   if (TESTBIT(fTreeMask, kSDigits)) {
     // AliInfo("Getting FMD summable digits");
-    if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) return kFALSE;
+    if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) { 
+      AliWarning("Failed to load SDigits!");
+      return kFALSE;
+    }
     fTreeS = fFMDLoader->TreeS();
     if (!fArrayS) fArrayS = fFMD->SDigits();
   }
@@ -428,10 +467,10 @@ AliFMDInput::Event()
     if (!ProcessTrackRefs()) return kFALSE; 
   if (TESTBIT(fTreeMask, kTracks)) 
     if (!ProcessTracks()) return kFALSE; 
-  if (TESTBIT(fTreeMask, kDigits)) 
-    if (!ProcessDigits()) return kFALSE;
   if (TESTBIT(fTreeMask, kSDigits)) 
     if (!ProcessSDigits()) return kFALSE;
+  if (TESTBIT(fTreeMask, kDigits)) 
+    if (!ProcessDigits()) return kFALSE;
   if (TESTBIT(fTreeMask, kRaw)) 
     if (!ProcessRawDigits()) return kFALSE;
   if (TESTBIT(fTreeMask, kRecPoints)) 
@@ -606,9 +645,13 @@ AliFMDInput::ProcessSDigits()
   Int_t nEv = fTreeS->GetEntries();
   for (Int_t i = 0; i < nEv; i++) {
     Int_t sdigitRead  = fTreeS->GetEntry(i);
-    if (sdigitRead <= 0) continue;
-    Int_t nSdigit = fArrayS->GetEntries();
+    if (sdigitRead <= 0) { 
+      AliInfo(Form("Read nothing from tree"));
+      continue;
+    }
+    Int_t nSdigit = fArrayS->GetEntriesFast();
     AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
+    AliInfo(Form("Got %5d digits for this event", nSdigit));
     if (nSdigit <= 0) continue;
     for (Int_t j = 0; j < nSdigit; j++) {
       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
index 355c26dd1776678a58d70ce6c4bc15b36762c198..9fbe3dc2096cbd2a6de311360915bce16f2a5b9e 100644 (file)
@@ -327,7 +327,8 @@ protected:
 };
 
 inline Bool_t AliFMDInput::ProcessHit(AliFMDHit*,TParticle*) { return kTRUE; }
-inline Bool_t AliFMDInput::ProcessTrackRef(AliTrackReference* trackRef, TParticle* track) { return kTRUE; }
+inline Bool_t AliFMDInput::ProcessTrackRef(AliTrackReference*, 
+                                          TParticle*) { return kTRUE; }
 inline Bool_t AliFMDInput::ProcessTrack(Int_t,TParticle*,
                                        AliFMDHit*) { return kTRUE; }
 inline Bool_t AliFMDInput::ProcessDigit(AliFMDDigit*) { return kTRUE; }
index 332c31991f103ad334ef84c4ce3efaf3556abd55..74c44e133a6950f90d304a7fb026f95edafbe00a 100644 (file)
@@ -50,6 +50,7 @@
 #include "AliFMDDebug.h" // Better debug macros
 #include "AliFMDParameters.h"  // ALIFMDPARAMETERS_H
 #include "AliFMDDigit.h"       // ALIFMDDIGIT_H
+#include "AliFMDSDigit.h"      // ALIFMDSDIGIT_H
 #include "AliFMDRawStream.h"   // ALIFMDRAWSTREAM_H 
 #include "AliRawReader.h"      // ALIRAWREADER_H 
 #include "AliFMDRawReader.h"   // ALIFMDRAWREADER_H 
@@ -77,12 +78,17 @@ AliFMDRawReader::AliFMDRawReader(AliRawReader* reader, TTree* tree)
   : TTask("FMDRawReader", "Reader of Raw ADC values from the FMD"),
     fTree(tree),
     fReader(reader), 
-    fSampleRate(1),
+    // fSampleRate(1),
     fData(0),
     fNbytes(0), 
     fSeen()
 {
   // Default CTOR
+  for (Int_t i = 0; i < 3; i++) { 
+    fSampleRate[i]   = 0;
+    fZeroSuppress[i] = kFALSE;
+    fNoiseFactor[i]  = 1;
+  }
 }
 
 //____________________________________________________________________
@@ -104,13 +110,11 @@ AliFMDRawReader::Exec(Option_t*)
                   array->GetEntriesFast(), nWrite));
 }
 
-
 //____________________________________________________________________
 Bool_t
-AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng, 
-                           UShort_t& sec, UShort_t& str, 
-                           Short_t&  adc, Bool_t&   zs, 
-                           UShort_t& fac)
+AliFMDRawReader::NextSample(UShort_t& det, Char_t&   rng, UShort_t& sec, 
+                           UShort_t& str, UShort_t& sam, UShort_t& rat, 
+                           Short_t&  adc, Bool_t&   zs,  UShort_t& fac)
 {
   // Scan current event for next signal.   It returns kFALSE when
   // there's no more data in the event. 
@@ -118,12 +122,14 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng,
   static AliFMDParameters*   pars     = 0;
   static AliFMDAltroMapping* map      = 0;
   static Int_t               ddl      = -1;
-  static UInt_t              rate     = 0;
+  // static UInt_t              rate     = 0;
   static UShort_t            tdet     = 0;
   static Char_t              trng     = '\0';
   static UShort_t            tsec     = 0;
   static Short_t             tstr     = 0;   
   static Short_t             bstr     = -1;
+  static Short_t             tsam     = -1;   
+  static UInt_t              trate    = 0;
   static Int_t               hwaddr   = -1;
   static UShort_t            stripMin = 0;
   static UShort_t            stripMax = 0; // 127;
@@ -146,11 +152,12 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng,
 
     // Reset variables
     ddl    = -1;  
-    rate   = 0;   
+    trate  = 0;   
     tdet   = 0;   
     trng   = '\0';
     tsec   = 0;   
     tstr   = 0;  
+    tsam   = -1;
     hwaddr = -1;
   }
   do { 
@@ -163,14 +170,14 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng,
     Int_t thisDDL = stream.GetDDLNumber();
     AliFMDDebug(10, ("RCU @ DDL %d", thisDDL));
     if (thisDDL != ddl) { 
-      ddl  = thisDDL;
-      zs   = stream.GetZeroSupp();
-      fac  = stream.GetNPostsamples();
-      tdet = map->DDL2Detector(ddl);
-      rate = 0; // stream.GetNPresamples();
+      ddl   = thisDDL;
+      fZeroSuppress[ddl] = zs    = stream.GetZeroSupp();
+      fNoiseFactor[ddl]  = fac   = stream.GetNPostsamples();
+      fSampleRate[ddl]   = trate = 0; // stream.GetNPresamples();
+      tdet  = map->DDL2Detector(ddl);
       AliFMDDebug(10, ("RCU @ DDL %d zero suppression: %s",ddl, zs?"yes":"no"));
       AliFMDDebug(10, ("RCU @ DDL %d noise factor: %d", ddl,fac));
-      AliFMDDebug(10, ("RCU @ DDL %d sample rate: %d", ddl, rate));
+      AliFMDDebug(10, ("RCU @ DDL %d sample rate: %d", ddl, trate));
     }
     Int_t thisAddr = stream.GetHWAddress();
     AliFMDDebug(10, ("RCU @ DDL %d, Address 0x%03x", ddl, thisAddr));    
@@ -185,9 +192,10 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng,
       stripMin = pars->GetMinStrip(tdet, trng, tsec, bstr);
       stripMax = pars->GetMaxStrip(tdet, trng, tsec, bstr);
       preSamp  = pars->GetPreSamples(tdet, trng, tsec, bstr);
-      if (rate == 0) rate = pars->GetSampleRate(tdet, trng, tsec, bstr);
+      if (trate == 0) 
+       fSampleRate[ddl] = trate = pars->GetSampleRate(tdet, trng, tsec, bstr);
       AliFMDDebug(10, ("RCU @ DDL %d, Address 0x%03x sample rate: %d", 
-                     ddl, hwaddr, rate));
+                     ddl, hwaddr, trate));
       
       Int_t nChAddrMismatch = stream.GetNChAddrMismatch();
       Int_t nChLenMismatch  = stream.GetNChLengthMismatch();
@@ -220,22 +228,12 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng,
 
     Short_t  strOff = 0;
     UShort_t samp   = 0;
-    map->Timebin2Strip(tsec, t, preSamp, rate, strOff, samp);
+    map->Timebin2Strip(tsec, t, preSamp, trate, strOff, samp);
     tstr = bstr + strOff;
+    tsam = samp;
     AliFMDDebug(20, ("0x%04x/0x%03x/%04d maps to FMD%d%c[%2d,%3d]-%d", 
                     ddl, hwaddr, t, tdet, trng, tsec, tstr, samp));
     
-    Bool_t take = kFALSE;
-    switch (rate) { 
-    case 1:                      take = kTRUE; break;
-    case 2:  if (samp == 1)      take = kTRUE; break;
-    case 3:  if (samp == 1)      take = kTRUE; break; 
-    case 4:  if (samp == 2)      take = kTRUE; break;
-    default: if (samp == rate-2) take = kTRUE; break;
-    }
-    if (!take) continue;
-
-
     // Local strip number to channel
     Short_t l = (tstr > bstr ? tstr - bstr : bstr - tstr);
     AliFMDDebug(10, ("Checking if strip %d in range [%d,%d]", 
@@ -252,6 +250,8 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng,
     rng = trng;
     sec = tsec;
     str = tstr;
+    sam = tsam;
+    rat = trate;
     // adc = stream.GetSignal();
     
     break;
@@ -259,6 +259,49 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng,
   return kTRUE;
 }
 
+//____________________________________________________________________
+Bool_t
+AliFMDRawReader::NextSignal(UShort_t& det, Char_t&   rng, 
+                           UShort_t& sec, UShort_t& str, 
+                           Short_t&  adc, Bool_t&   zs, 
+                           UShort_t& fac)
+{
+  
+  do { 
+    UShort_t samp, rate;
+    if (!NextSample(det, rng, sec, str, samp, rate, adc, zs, fac)) 
+      return kFALSE;
+
+    Bool_t take = kFALSE;
+    switch (rate) { 
+    case 1:                      take = kTRUE; break;
+    case 2:  if (samp == 1)      take = kTRUE; break;
+    case 3:  if (samp == 1)      take = kTRUE; break; 
+    case 4:  if (samp == 2)      take = kTRUE; break;
+    default: if (samp == rate-2) take = kTRUE; break;
+    }
+    if (!take) continue;
+    break;
+  } while (true);
+  return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDRawReader::SelectSample(UShort_t samp, UShort_t rate) 
+{
+  Bool_t take = kFALSE;
+  switch (rate) { 
+  case 1:                      take = kTRUE; break;
+  case 2:  if (samp == 1)      take = kTRUE; break;
+  case 3:  if (samp == 1)      take = kTRUE; break; 
+  case 4:  if (samp == 2)      take = kTRUE; break;
+  default: if (samp == rate-2) take = kTRUE; break;
+  }
+  
+  return take;
+}
+  
 #if 1
 //____________________________________________________________________
 Bool_t
@@ -291,7 +334,7 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
 
   Int_t  oldddl = -1;
   UInt_t ddl    = 0;
-  UInt_t rate   = 0;
+  // UInt_t rate   = 0;
   UInt_t last   = 0;
   UInt_t hwaddr = 0;
   // Data array is approx twice the size needed. 
@@ -318,6 +361,12 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
       fNoiseFactor[ddl]  = input.GetNPostsamples();
       AliFMDDebug(20, ("RCU @ DDL %d noise factor: %d", ddl,fNoiseFactor[ddl]));
 
+      // WARNING: We store the noise factor in the 2nd baseline
+      // filters excluded post samples, since we'll never use that
+      // mode. 
+      fSampleRate[ddl]     = input.GetNPretriggerSamples();
+      AliFMDDebug(20, ("RCU @ DDL %d Sample rate: %d", ddl,fNoiseFactor[ddl]));
+
       Int_t nChAddrMismatch = input.GetNChAddrMismatch();
       Int_t nChLenMismatch  = input.GetNChLengthMismatch();
       if (nChAddrMismatch != 0) 
@@ -343,16 +392,14 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
                    "hardware address 0x%03x", ddl, hwaddr));
       continue;
     }
-    AliFMDDebug(1, ("Board: 0x%02x, Altro: 0x%x, Channel: 0x%x, Length: %4d", 
+    AliFMDDebug(5, ("Board: 0x%02x, Altro: 0x%x, Channel: 0x%x, Length: %4d", 
                    board, chip, channel, last));
 
     stripMin = pars->GetMinStrip(det, ring, sec, strbase);
     stripMax = pars->GetMaxStrip(det, ring, sec, strbase);
     preSamp  = pars->GetPreSamples(det, ring, sec, strbase);
-    // WARNING: We use the number of pre-samples to store the
-    // oversampling rate in. 
-    rate     = input.GetNPretriggerSamples();
-    if (rate == 0) rate = pars->GetSampleRate(det, ring, sec, strbase);
+    if (fSampleRate[ddl] == 0) 
+      fSampleRate[ddl] = pars->GetSampleRate(det, ring, sec, strbase);
     
     // Loop over the `timebins', and make the digits
     for (size_t i = 0; i < last; i++) {
@@ -360,7 +407,7 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
       AliFMDDebug(15, ("0x%04x/0x%03x/%04d %4d", ddl, hwaddr, i, data[i]));
 
       Short_t  stroff = 0;
-      map->Timebin2Strip(sec, i, preSamp, rate, stroff, samp);
+      map->Timebin2Strip(sec, i, preSamp, fSampleRate[ddl], stroff, samp);
       Short_t  str    = strbase + stroff;
       
       AliFMDDebug(10, ("0x%04x/0x%03x/%04d maps to FMD%d%c[%2d,%3d]-%d", 
@@ -371,7 +418,7 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
        continue;
       }
       
-      Short_t lstrip = (i - preSamp) / rate + stripMin;
+      Short_t lstrip = (i - preSamp) / fSampleRate[ddl] + stripMin;
       
       AliFMDDebug(15, ("Checking if strip %d (%d) in range [%d,%d]", 
                      lstrip, str, stripMin, stripMax));
@@ -381,6 +428,9 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
        data[i] = 0; // Reset cache 
        continue;
       }
+      // Possibly do pedestal subtraction of signal 
+      Int_t counts = data[i];
+       
       
       // Check the cache of indicies
       Int_t idx = fSeen(det, ring, sec, str);
@@ -391,11 +441,11 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
                       det, ring, sec, str, samp, i));
        new ((*array)[idx]) AliFMDDigit(det, ring, sec, str);
       }
-      AliFMDDigit* digit = static_cast<AliFMDDigit*>(array->At(idx));
+      AliFMDBaseDigit* digit = static_cast<AliFMDBaseDigit*>(array->At(idx));
       AliFMDDebug(10,
-                 ("Setting from FMD%d%c[%2d,%3d]-%d from timebin %4d = %4d", 
-                  det, ring, sec, str, samp, i, data[i]));
-      digit->SetCount(samp, data[i]);
+                 ("Setting FMD%d%c[%2d,%3d]-%d from timebin %4d=%4d (%4d)", 
+                  det, ring, sec, str, samp, i, counts, data[i]));
+      digit->SetCount(samp, counts);
       data[i] = 0; // Reset cache 
     }
   }
index b1e2a537fddf2e7d18cf5ec117d2edc1cec5655e..2499ff8326c659455ebe6158f1fcf253eb019cf0 100644 (file)
@@ -47,18 +47,34 @@ class AliFMDCalibStripRange;
 class AliFMDRawReader : public TTask 
 {
 public:
-  /** CTOR 
-      @param reader Raw reader
-      @param array  Output tree */
+  /** 
+   * CTOR 
+   *
+   * @param reader Raw reader
+   * @param array  Output tree 
+   */
   AliFMDRawReader(AliRawReader* reader, TTree* array);
-  /** DTOR */
+  /** 
+   * DTOR 
+   */
   virtual ~AliFMDRawReader() {}
-  /** Read in, and store in output tree 
-      @param option Not used */
+  /** 
+   * Read in, and store in output tree 
+   *
+   * @param option Not used 
+   */
   virtual void   Exec(Option_t* option="");
-  /** Read ADC's into a TClonesArray of AliFMDDigit objects. 
-      @param array Array to read into 
-      @return @c true on success */
+  /**
+   * Read ADC's into a TClonesArray of AliFMDDigit objects. 
+   *
+   * @param array       Array to read into 
+   * @param summable    Create SDigits rather than digits 
+   * @param pedSub      Whether to do pedestal subtraction.
+   * @param noiseFactor If we do pedestal subtraction, then this is
+   *        the number we use to suppress remenants of the noise.
+   * 
+   * @return @c true on success 
+   */
   virtual Bool_t ReadAdcs(TClonesArray* array);
   /** 
    * Read SOD event into passed objects.
@@ -92,16 +108,55 @@ public:
    */
   UShort_t NoiseFactor(UShort_t ddl) const { return fNoiseFactor[ddl]; }
 
+  /** 
+   * Get the next signal
+   * 
+   * @param det  On return, the detector
+   * @param rng  On return, the ring
+   * @param sec  On return, the sector
+   * @param str  On return, the strip
+   * @param sam  On return, the sample
+   * @param rat  On return, the sample rate
+   * @param adc  On return, the ADC value
+   * @param zs   On return, whether zero-supp. is enabled
+   * @param fac  On return, the usd noise factor
+   * 
+   * @return true if valid data is returned
+   */
+  Bool_t NextSample(UShort_t& det, Char_t&   rng, UShort_t& sec, UShort_t& str,
+                   UShort_t& sam, UShort_t& rat, Short_t&  adc, 
+                   Bool_t&   zs,  UShort_t& fac);
+  /** 
+   * Get the next signal
+   * 
+   * @param det  On return, the detector
+   * @param rng  On return, the ring
+   * @param sec  On return, the sector
+   * @param str  On return, the strip
+   * @param adc  On return, the ADC value
+   * @param zs   On return, whether zero-supp. is enabled
+   * @param fac  On return, the usd noise factor
+   * 
+   * @return true if valid data is returned
+   */
   Bool_t NextSignal(UShort_t& det, Char_t&   rng, 
                    UShort_t& sec, UShort_t& str, 
                    Short_t&  adc, Bool_t&   zs, 
                    UShort_t& fac);
+  /** 
+   * Whether to keep a sample based on the rate used. 
+   * 
+   * @param samp Sample number 
+   * @param rate Over sampling rate
+   * 
+   * @return Whether to keep the sample or not
+   */
+  static Bool_t SelectSample(UShort_t samp, UShort_t rate);
 protected:
   AliFMDRawReader(const AliFMDRawReader& o) 
     : TTask(o), 
       fTree(0), 
       fReader(0), 
-      fSampleRate(0),
       fData(0),
       fNbytes(0), 
       fSeen()
@@ -112,7 +167,7 @@ protected:
   Int_t GetHalfringIndex(UShort_t det, Char_t ring, UShort_t board);
   TTree*          fTree;       //! Pointer to tree to read into 
   AliRawReader*   fReader;     //! Pointer to raw reader 
-  UShort_t        fSampleRate; // The sample rate (if 0, inferred from data)
+  UShort_t        fSampleRate[3]; // The sample rate (if 0, inferred from data)
   UChar_t*        fData; 
   ULong_t        fNbytes; 
   Bool_t          fZeroSuppress[3];
index 0bd09864de5deed6c4e14c04d1423618a5b458bb..5e9cab6d9ea9d9c52129cb2a2db084df79abd7ec 100644 (file)
@@ -39,6 +39,7 @@
 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
 #include "AliFMDAltroMapping.h"            // ALIFMDALTROMAPPING_H
 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
+#include "AliFMDSDigit.h"                  // ALIFMDDIGIT_H
 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
 #include "AliFMDRecPoint.h"               // ALIFMDMULTNAIIVE_H
@@ -205,7 +206,7 @@ AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
                                   TTree* digitsTree) const
 {
   // Convert Raw digits to AliFMDDigit's in a tree 
-  AliFMDDebug(2, ("Reading raw data into digits tree"));
+  AliFMDDebug(1, ("Reading raw data into digits tree"));
   AliFMDRawReader rawRead(reader, digitsTree);
   // rawRead.SetSampleRate(fFMD->GetSampleRate());
   rawRead.Exec();
@@ -276,6 +277,34 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
   }
 }
 
+//____________________________________________________________________
+void 
+AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
+{
+  // Reconstruct directly from raw data (no intermediate output on
+  // digit tree or rec point tree).  
+  // 
+  // Parameters: 
+  //   reader  Raw event reader 
+  //   ctree    Not used. 
+  AliFMDRawReader rawReader(reader, 0);
+
+  UShort_t det, sec, str, sam, rat, fac;
+  Short_t  adc, oldDet = -1;
+  Bool_t   zs;
+  Char_t   rng;
+    
+  while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) { 
+    if (!rawReader.SelectSample(sam, rat)) continue;
+    if (det != oldDet) { 
+      fZS[det-1]       = zs;
+      fZSFactor[det-1] = fac;
+      oldDet           = det;
+    }
+    DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
+  }
+}
+
 //____________________________________________________________________
 void 
 AliFMDReconstructor::Reconstruct(TTree* digitsTree, 
@@ -311,7 +340,7 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   AliFMDDebug(5, ("Getting entry 0 from digit branch"));
   digitBranch->GetEntry(0);
   
-  AliFMDDebug(5, ("Processing digits"));
+  AliFMDDebug(1, ("Processing digits"));
   ProcessDigits(digits);
 
   Int_t written = clusterTree->Fill();
@@ -418,6 +447,99 @@ AliFMDReconstructor::ProcessSignal(UShort_t det,
 
 }
 
+//____________________________________________________________________
+void
+AliFMDReconstructor::DigitizeSignal(TClonesArray* sdigits, 
+                                   UShort_t det, 
+                                   Char_t   rng, 
+                                   UShort_t sec, 
+                                   UShort_t str, 
+                                   UShort_t /* sam */,
+                                   Short_t  adc) const
+{
+  // Process the signal from a single strip 
+  // 
+  // Parameters: 
+  //    det    Detector ID
+  //    rng    Ring ID
+  //    sec    Sector ID
+  //    rng    Strip ID
+  //    adc     ADC counts
+  // 
+  AliFMDParameters* param  = AliFMDParameters::Instance();
+  // Check that the strip is not marked as dead 
+  if (param->IsDead(det, rng, sec, str)) {
+    AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
+    return;
+  }
+  
+  // Substract pedestal. 
+  UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
+  if(counts == USHRT_MAX || counts == 0) return;
+  
+    // Gain match digits. 
+  Double_t edep     = Adc2Energy(det, rng, sec, str, counts);
+  // Get rid of nonsense energy
+  if(edep < 0)  return;
+
+  Int_t n = sdigits->GetEntriesFast();
+  // AliFMDSDigit* sdigit = 
+  new ((*sdigits)[n]) 
+    AliFMDSDigit(det, rng, sec, str, edep, counts, counts, counts, counts);
+  // sdigit->SetCount(sam, counts);
+}
+
+//____________________________________________________________________
+UShort_t
+AliFMDReconstructor::SubtractPedestal(UShort_t det, 
+                                     Char_t   rng, 
+                                     UShort_t sec, 
+                                     UShort_t str, 
+                                     UShort_t adc, 
+                                     Float_t  noiseFactor,
+                                     Bool_t   zsEnabled, 
+                                     UShort_t zsNoiseFactor) const
+{
+  AliFMDParameters* param  = AliFMDParameters::Instance();
+  Float_t           ped    = (zsEnabled ? 0 : 
+                               param->GetPedestal(det, rng, sec, str));
+  Float_t           noise  = param->GetPedestalWidth(det, rng, sec, str);
+  if(ped < 0 || noise < 0) { 
+    AliWarningClass(Form("Invalid pedestal (%f) or noise (%f) "
+                        "for FMD%d%c[%02d,%03d]", 
+                   ped, noise, det, rng, sec, str));
+    return USHRT_MAX;
+  }
+  AliDebugClass(15, Form("Subtracting pedestal for FMD%d%c[%2d,%3d]=%4d "
+                        "(%s w/factor %d, noise factor %f, "
+                        "pedestal %8.2f+/-%8.2f)",
+                        det, rng, sec, str, adc, 
+                        (zsEnabled ? "zs'ed" : "straight"), 
+                        zsNoiseFactor, noiseFactor, ped, noise));
+
+  Int_t counts = adc + Int_t(zsEnabled ? zsNoiseFactor * noise : - ped);
+  counts =  TMath::Max(Int_t(counts), 0);
+  // Calculate the noise factor for suppressing remenants of the noise
+  // peak.  If we have done on-line zero suppression, we only check
+  // for noise signals that are larger than the suppressed noise.  If
+  // the noise factor used on line is larger than the factor used
+  // here, we do not do this check at all.  
+  // 
+  // For example:
+  //    Online factor  |  Read factor |  Result 
+  //    ---------------+--------------+-------------------------------
+  //           2       |      3       | Check if signal > 1 * noise
+  //           3       |      3       | Check if signal > 0
+  //           3       |      2       | Check if signal > 0
+  //
+  // In this way, we make sure that we do not suppress away too much
+  // data, and that the read-factor is the most stringent cut. 
+  Float_t nf = TMath::Max(0.F, noiseFactor - (zsEnabled ? zsNoiseFactor : 0));
+  if (counts < noise * nf) counts = 0;
+  if (counts > 0) AliDebugClass(15, "Got a hit strip");
+
+  return counts;
+}
 
 
 //____________________________________________________________________
@@ -439,29 +561,12 @@ AliFMDReconstructor::SubtractPedestal(UShort_t det,
   // Return:
   //    Pedestal subtracted signal or USHRT_MAX in case of problems 
   //
-  AliFMDParameters* param  = AliFMDParameters::Instance();
-  Bool_t            zs     = fZS[det-1];
-  UShort_t          fac    = fZSFactor[det-1];
-  Float_t           ped    = (zs ? 0 : 
-                             param->GetPedestal(det, rng, sec, str));
-  Float_t           noise  = param->GetPedestalWidth(det, rng, sec, str);
-  if(ped < 0 || noise < 0) { 
-    AliWarning(Form("Invalid pedestal (%f) or noise (%f) "
-                   "for FMD%d%c[%02d,%03d]", ped, noise, det, rng, sec, str));
-    return USHRT_MAX;
-  }
-
-  AliFMDDebug(5, ("Subtracting pedestal %f from signal %d", ped, adc));
-  // if (digit->Count3() > 0)      adc = digit->Count3();
-  // else if (digit->Count2() > 0) adc = digit->Count2();
-  // else                          adc = digit->Count1();
-  Int_t counts = adc + Int_t(zs ? fac * noise : - ped);
-  counts       = TMath::Max(Int_t(counts), 0);
-  if (counts < noise * fNoiseFactor) counts = 0;
-  if (counts > 0) AliFMDDebug(15, ("Got a hit strip"));
+  UShort_t counts = SubtractPedestal(det, rng, sec, str, adc, 
+                                    fNoiseFactor, fZS[det-1], 
+                                    fZSFactor[det-1]);
   if (fDiagStep1) fDiagStep1->Fill(adc, counts);
   
-  return  UShort_t(counts);
+  return counts;
 }
 
 //____________________________________________________________________
@@ -470,7 +575,6 @@ AliFMDReconstructor::Adc2Energy(UShort_t det,
                                Char_t   rng, 
                                UShort_t sec, 
                                UShort_t str, 
-                               Float_t  eta, 
                                UShort_t count) const
 {
   // Converts number of ADC counts to energy deposited. 
@@ -504,7 +608,6 @@ AliFMDReconstructor::Adc2Energy(UShort_t det,
   //    rng    Ring ID
   //    sec    Sector ID
   //    rng    Strip ID
-  //    eta     Psuedo-rapidity
   //    counts  Number of ADC counts over pedestal
   // Return 
   //    The energy deposited in a single strip, or -1 in case of problems
@@ -523,6 +626,56 @@ AliFMDReconstructor::Adc2Energy(UShort_t det,
 
   Double_t edep  = ((count * param->GetEdepMip()) 
                    / (gain * param->GetDACPerMIP()));
+  return edep;
+}
+
+//____________________________________________________________________
+Float_t
+AliFMDReconstructor::Adc2Energy(UShort_t det, 
+                               Char_t   rng, 
+                               UShort_t sec, 
+                               UShort_t str, 
+                               Float_t  eta, 
+                               UShort_t count) const
+{
+  // Converts number of ADC counts to energy deposited. 
+  // Note, that this member function can be overloaded by derived
+  // classes to do strip-specific look-ups in databases or the like,
+  // to find the proper gain for a strip. 
+  // 
+  // In the first simple version, we calculate the energy deposited as 
+  // 
+  //    EnergyDeposited = cos(theta) * gain * count
+  // 
+  // where 
+  // 
+  //           Pre_amp_MIP_Range
+  //    gain = ----------------- * Energy_deposited_per_MIP
+  //           ADC_channel_size    
+  // 
+  // is constant and the same for all strips.
+  //
+  // For the production we use the conversion measured in the NBI lab.
+  // The total conversion is then:
+  // 
+  //    gain = ADC / DAC
+  // 
+  //                  EdepMip * count
+  //      => energy = ----------------
+  //                  gain * DACPerADC
+  // 
+  // Parameters: 
+  //    det    Detector ID
+  //    rng    Ring ID
+  //    sec    Sector ID
+  //    rng    Strip ID
+  //    eta     Psuedo-rapidity
+  //    counts  Number of ADC counts over pedestal
+  // Return 
+  //    The energy deposited in a single strip, or -1 in case of problems
+  //
+  Double_t edep = Adc2Energy(det, rng, sec, str, count);
+  
   if (fDiagStep2) fDiagStep2->Fill(count, edep);  
   if (fAngleCorrect) {
     Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
index 446b32e3c524bd1436299bdec9a499f5c1266e79..a04989d0c16086ea7816d6bfab3dc4c07fed568e 100644 (file)
@@ -117,6 +117,16 @@ public:
    */
   virtual void   FillESD(AliRawReader*, TTree* clusterTree, 
                         AliESDEvent* esd) const;
+
+  /** 
+   * Create SDigits from raw data
+   * 
+   * @param reader  The raw reader
+   * @param sdigits Array to fill with AliFMDSDigit objects. 
+   */  
+  virtual void Digitize(AliRawReader* reader, 
+                       TClonesArray* sdigits) const;
+  
   /** 
    * Not used 
    */
@@ -199,6 +209,49 @@ protected:
                             UShort_t sec, 
                             UShort_t str, 
                             Short_t  adc) const;
+  /** 
+   * Process the signal from a single strip. 
+   * 
+   * @param sdigits Array to fill
+   * @param det     Detector number 
+   * @param rng     Ring identifier 
+   * @param sec     Sector number
+   * @param str     Strip number 
+   * @param sam     Sample number 
+   * @param adc     Number of ADC counts for this strip
+   */  
+  virtual void DigitizeSignal(TClonesArray* sdigits, 
+                             UShort_t      det, 
+                             Char_t         rng, 
+                             UShort_t       sec, 
+                             UShort_t       str, 
+                             UShort_t       sam,
+                             Short_t        adc) const;
+  /** 
+   * Subtract the pedestal off the ADC counts. 
+   * 
+   * @param det           Detector number
+   * @param rng           Ring identifier
+   * @param sec           Sector number
+   * @param str           Strip number
+   * @param adc           ADC counts
+   * @param noiseFactor   If pedestal substracted pedestal is less then
+   *        this times the noise, then consider this to be 0. 
+   * @param zsEnabled     Whether zero-suppression is on.
+   * @param zsNoiseFactor Noise factor used in on-line pedestal
+   *        subtraction. 
+   * 
+   * @return The pedestal subtracted ADC counts (possibly 0), or @c
+   *         USHRT_MAX in case of problems.
+   */  
+  virtual UShort_t SubtractPedestal(UShort_t det, 
+                                   Char_t   rng, 
+                                   UShort_t sec, 
+                                   UShort_t str, 
+                                   UShort_t adc, 
+                                   Float_t  noiseFactor,
+                                   Bool_t   zsEnabled, 
+                                   UShort_t zsNoiseFactor) const;
   /** 
    * Substract pedestals from raw ADC in @a digit
    * 
@@ -233,6 +286,29 @@ protected:
    *
    * @return Energy deposited @f$ E_i@f$ 
    */
+  virtual Float_t  Adc2Energy(UShort_t det, 
+                             Char_t   rng, 
+                             UShort_t sec, 
+                             UShort_t str, 
+                             UShort_t count) const;
+  /** 
+   * Converts number of ADC counts to energy deposited.   This is
+   * done by 
+   * @f[
+   * E_i = A_i g_i
+   * @f]
+   * where @f$ A_i@f$ is the pedestal subtracted ADC counts, and @f$
+   * g_i@f$ is the gain for the @f$ i^{\mbox{th}}@f$ strip. 
+   * 
+   * @param det          Detector number  
+   * @param rng   Ring identifier 
+   * @param sec   Sector number
+   * @param str   Strip number 
+   * @param eta   Psuedo-rapidity of digit.
+   * @param count Pedestal subtracted ADC counts
+   *
+   * @return Energy deposited @f$ E_i@f$ 
+   */
   virtual Float_t  Adc2Energy(UShort_t det, 
                              Char_t   rng, 
                              UShort_t sec, 
index 9077fc5cec7ab90366cb13545be1f8178046a236..e4aeabdf2e57b10117db781232c912b330df0155 100644 (file)
@@ -79,7 +79,9 @@ AliFMDSDigit::AliFMDSDigit()
     fCount2(-1),
     fCount3(-1), 
     fCount4(-1), 
-    fLabels(0)
+    fNParticles(0),
+    fNPrimaries(0)
+    // ,     fLabels(0)
 {
   // cTOR 
 }
@@ -104,8 +106,8 @@ AliFMDSDigit::AliFMDSDigit(UShort_t       detector,
     fCount3(count3),
     fCount4(count4),
     fNParticles(npart), 
-    fNPrimaries(nprim)
-    fLabels(refs)
+    fNPrimaries(nprim)
+    // , fLabels(refs)
 {
   //
   // Creates a real data digit object
@@ -120,6 +122,7 @@ AliFMDSDigit::AliFMDSDigit(UShort_t       detector,
   //    count1    ADC count (a 10-bit word)
   //    count2    ADC count (a 10-bit word) -1 if not used
   //    count3    ADC count (a 10-bit word) -1 if not used
+  for (Int_t i = 0; i < refs.fN; i++) AddTrack(refs.fArray[i]);
 }
 
 //____________________________________________________________________
@@ -128,17 +131,28 @@ AliFMDSDigit::Print(Option_t* option) const
 {
   // Print digit to standard out 
   AliFMDBaseDigit::Print();
-  std::cout << "\t" << fEdep << " -> "
-           << fCount1 << " (" << fCount2 << "," << fCount3 << "," 
-           << fCount4 << ") = " << Counts() << std::flush;
+  std::cout << "\t" 
+           << std::setw(10) << fEdep << " -> "
+           << std::setw(4) << fCount1 << " (" 
+           << std::setw(4) << fCount2 << "," 
+           << std::setw(4) << fCount3 << "," 
+           << std::setw(4) << fCount4 << ") = " 
+           << std::setw(4) << Counts() << std::flush;
+
   TString opt(option);
   if (opt.Contains("p", TString::kIgnoreCase)) 
-    std::cout << " [" << fNPrimaries << "/" << fNParticles << "]"
+    std::cout << " [" 
+             << std::setw(2) << fNPrimaries << "/" 
+             << std::setw(2) << fNParticles << "]"
              << std::flush;
   if (opt.Contains("l", TString::kIgnoreCase)) {
     std::cout << " ";
+    for (Int_t i = 0; i < GetNTrack(); i++) 
+      std::cout << (i == 0 ? "" : ",") << std::setw(5) << fTracks[i];
+#if 0
     for (Int_t i = 0; i < fLabels.fN; i++) 
       std::cout << (i == 0 ? "" : ",") << fLabels.fArray[i];
+#endif
   }
   std::cout << std::endl;
 }
index 18598eeeddf9f73b9785142b2f6ad6596a6067bb..c786cc264014b63a559ba3e999202d8b311b46b1 100644 (file)
 class AliFMDSDigit : public AliFMDBaseDigit
 {
 public:
-  /** CTOR */
+  /** 
+   * CTOR 
+   */
   AliFMDSDigit();
-  /** Constrctor 
-      @param detector Detector 
-      @param ring     Ring
-      @param sector   Sector
-      @param strip    Strip 
-      @param edep     Energy deposited 
-      @param count    ADC (first sample)
-      @param count2   ADC (second sample, or -1 if not used)
-      @param count3   ADC (third sample, or -1 if not used) */
+  /** 
+   * Constrctor 
+   *
+   * @param detector Detector 
+   * @param ring     Ring
+   * @param sector   Sector
+   * @param strip    Strip 
+   * @param edep     Energy deposited 
+   * @param count    ADC (first sample)
+   * @param count2   ADC (second sample, or -1 if not used)
+   * @param count3   ADC (third sample, or -1 if not used) 
+   */
   AliFMDSDigit(UShort_t       detector, 
               Char_t         ring='\0', 
               UShort_t       sector=0, 
@@ -50,29 +55,67 @@ public:
               UShort_t       npart=0,
               UShort_t       nprim=0, 
               const TArrayI& refs=TArrayI());
-  /** DTOR */
+  /** 
+   * DTOR 
+   */
   virtual ~AliFMDSDigit() {}
-  /** @return ADC count (first sample) */
-  UShort_t Count1()                const { return fCount1;   }
-  /** @return ADC count (second sample, or -1 if not used) */
-  Short_t  Count2()                const { return fCount2;   }
-  /** @return ADC count (third sample, or -1 if not used) */
-  Short_t  Count3()                const { return fCount3;   }
-  /** @return ADC count (third sample, or -1 if not used) */
-  Short_t  Count4()                const { return fCount4;   }
-  /** @return Canonical ADC counts */
-  UShort_t Counts()                const;
-  /** @return Energy deposited */
-  Float_t  Edep()                  const { return fEdep;     }
-  /** @return Number of particles that hit this strip */
+  /** 
+   * 
+   * @return ADC count (first sample) 
+   */
+  UShort_t Count1() const { return fCount1;   }
+  /** 
+   * 
+   * @return ADC count (second sample, or -1 if not used) 
+   */
+  Short_t  Count2() const { return fCount2;   }
+  /** 
+   * 
+   * @return ADC count (third sample, or -1 if not used) 
+   */
+  Short_t  Count3() const { return fCount3;   }
+  /** 
+   * 
+   * @return ADC count (third sample, or -1 if not used) 
+   */
+  Short_t  Count4() const { return fCount4;   }
+  /** 
+   *
+   * @return Canonical ADC counts 
+   */
+  UShort_t Counts() const;
+  /** 
+   * 
+   * @return Energy deposited 
+   */
+  Float_t  Edep() const { return fEdep;     }
+  /** 
+   * 
+   * @return Number of particles that hit this strip 
+   */
   UShort_t NParticles() const { return fNParticles; }
-  /** @return Number of primary particles that hit this strip */
+  /** 
+   * 
+   * @return Number of primary particles that hit this strip 
+   */
   UShort_t NPrimaries() const { return fNPrimaries; }
+#if 0
   /** @return the track labels */
   const TArrayI& TrackLabels() const { return fLabels; }
-  /** Print info 
-      @param opt Not used */
+#endif
+  /** 
+   * Print info 
+   *
+   * @param opt Not used 
+   */
   void     Print(Option_t* opt="") const;
+  /** 
+   * Set the count value 
+   * 
+   * @param s Sample number 
+   * @param c Counts 
+   */
+  void SetCount(UShort_t s, Short_t c);
 protected:
   Float_t  fEdep;       // Energy deposited 
   UShort_t fCount1;     // Digital signal 
@@ -81,8 +124,10 @@ protected:
   Short_t  fCount4;     // Digital signal (-1 if not used)
   UShort_t fNParticles; // Total number of particles that hit this strip
   UShort_t fNPrimaries; // Number of primary particles that his this strip
+#if 0
   TArrayI  fLabels;     // MC-truth track labels
-  ClassDef(AliFMDSDigit,4)     // Summable FMD digit
+#endif
+  ClassDef(AliFMDSDigit,5)     // Summable FMD digit
 };
   
 inline UShort_t 
@@ -94,6 +139,16 @@ AliFMDSDigit::Counts() const
   return fCount1;
 }
 
+inline void
+AliFMDSDigit::SetCount(UShort_t i, Short_t c)
+{
+  switch (i) {
+  case 0: fCount1 = c; break;
+  case 1: fCount2 = c; break;
+  case 2: fCount3 = c; break;
+  case 3: fCount4 = c; break;
+  }
+}
 
 #endif
 //____________________________________________________________________
index 3ca4fbca7594ace925d6075b5efaebb4b789f6d3..b6d524cff9942f151dc52c480190394d2f1eca93 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-//Struct to encode a strip address into one integer
-//developed by Christian Holm Christensen (cholm@nbi.dk).
-// The functions are static
-// to ensure applicability from anywhere. This is needed to smoothly store
-//strip addresses in track references.
+// Struct to encode a strip address into one integer
+// developed by Christian Holm Christensen (cholm@nbi.dk).
+// 
+// The functions are static to ensure applicability from
+// anywhere. This is needed to smoothly store strip addresses in track
+// references. 
+//
 // Added by Hans H. Dalsgaard (hans.dalsgaard@cern.ch) 
 
 
-struct AliFMDStripIndex
+class AliFMDStripIndex
 {
+public:
+  /** 
+   * Constructor
+   * 
+   */  
+  AliFMDStripIndex() {}
+  /** 
+   * Destructor 
+   * 
+   */
+  virtual ~AliFMDStripIndex() {}
+  /** 
+   * Pack an identifier from detector coordinates
+   * 
+   * @param det  Detector
+   * @param rng  Ring
+   * @param sec  Sector
+   * @param str  Strip
+   * 
+   * @return Packed identifier 
+   */
   static UInt_t Pack(UShort_t det, Char_t rng, UShort_t sec, UShort_t str) 
   {
     UInt_t irg  = (rng == 'I' || rng == 'i' ? 0 : 1);
@@ -34,8 +57,17 @@ struct AliFMDStripIndex
                   ((det & 0x003) << 17));
     return id;
   }
+  /** 
+   * Unpack an identifier to detector coordinates
+   * 
+   * @param id   Identifier to unpack
+   * @param det  On return, the detector
+   * @param rng  On return, the ring
+   * @param sec  On return, the sector
+   * @param str  On return, the strip
+   */  
   static void Unpack(UInt_t id, 
-             UShort_t& det, Char_t& rng, UShort_t& sec, UShort_t& str)
+                    UShort_t& det, Char_t& rng, UShort_t& sec, UShort_t& str)
   {
     str = ((id >>  0) & 0x1FF);
     sec = ((id >>  9) & 0x03F);
@@ -43,6 +75,11 @@ struct AliFMDStripIndex
     det = ((id >> 17) & 0x003);
   }
 
-  ClassDef(AliFMDStripIndex,0)
+  ClassDef(AliFMDStripIndex,1)
 };
 #endif
+//
+// Local Variables:
+//   mode: C++
+// End:
+//
index bac9fd7dc7f69756457fa80e928e7e7cc11205d9..3c43d8a2653359390e74e731052f858df979a48c 100644 (file)
 #include <TGeoPhysicalNode.h>
 #include "AliFMDGeometry.h"
 
+//____________________________________________________________________
+Double_t
+AliFMDSurveyToAlignObjs::GetUnitFactor() const
+{
+  // Returns the conversion factor from the measured values to
+  // centimeters. 
+  TString units(fSurveyObj->GetUnits());
+  if      (units.CompareTo("mm", TString::kIgnoreCase) == 0) return .1;
+  else if (units.CompareTo("cm", TString::kIgnoreCase) == 0) return 1.;
+  else if (units.CompareTo("m",  TString::kIgnoreCase) == 0) return 100.;
+  return 1;
+}
+
 //____________________________________________________________________
 Bool_t
-AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, 
-                                     Double_t* trans)
+AliFMDSurveyToAlignObjs::GetPoint(const char* name, 
+                                 TVector3&   point, 
+                                 TVector3&   error) const
 {
+  // Get named point.   On return, point will contain the point
+  // coordinates in centimeters, and error will contain the
+  // meassurement errors in centimeters too.  If no point is found,
+  // returns false, otherwise true. 
+  Double_t unit = GetUnitFactor();
+  TObject* obj  = fSurveyPoints->FindObject(name);
+  if (!obj) return kFALSE;
+  
+   AliSurveyPoint* p = static_cast<AliSurveyPoint*>(obj);
+   point.SetXYZ(unit * p->GetX(), 
+               unit * p->GetY(),
+               unit * p->GetZ());
+   error.SetXYZ(unit * p->GetPrecisionX(),
+               unit * p->GetPrecisionY(),
+               unit * p->GetPrecisionZ());
+   return kTRUE;
+}
 
-  // The possile survey points 
-  const char*  names[] = { "FMD2_ITOP",  "FMD2_OTOP", 
-                          "FMD2_IBOTM", "FMD2_OBOTM", 
-                          "FMD2_IBOT",  "FMD2_OBOT", 
-                          0 };
-  const char** name    = names;
-  Int_t    i = 0;
-  TGraph2DErrors g;
+//____________________________________________________________________
+Bool_t 
+AliFMDSurveyToAlignObjs::CalculatePlane(const     TVector3& a, 
+                                       const     TVector3& b,
+                                       const     TVector3& c, 
+                                       Double_t* trans,
+                                       Double_t* rot) const
+{
+  // Vector a->b, b->c, and normal to plane defined by these two
+  // vectors. 
+  TVector3 ab(b-a), bc(c-b);
   
-  Double_t unit = 1.;
-  TString units(fSurveyObj->GetUnits());
-  if      (units.CompareTo("mm", TString::kIgnoreCase) == 0) unit = .1;
-  else if (units.CompareTo("cm", TString::kIgnoreCase) == 0) unit = 1.;
-  else if (units.CompareTo("m",  TString::kIgnoreCase) == 0) unit = 100.;
+  // Normal vector to the plane of the fiducial marks obtained
+  // as cross product of the two vectors on the plane d0^d1
+  TVector3 nn(ab.Cross(bc));
+  if (nn.Mag() < 1e-8) { 
+    Info("DoIt", "Normal vector is null vector");
+    return kFALSE;
+  }
   
+  // We express the plane in Hessian normal form.
+  //
+  //   n x = -p,
+  //
+  // where n is the normalised normal vector given by 
+  // 
+  //   n_x = a / l,   n_y = b / l,   n_z = c / l,  p = d / l
+  //
+  // with l = sqrt(a^2+b^2+c^2) and a, b, c, and d are from the
+  // normal plane equation 
+  //
+  //  ax + by + cz + d = 0
+  // 
+  // Normalize
+  TVector3 n(nn.Unit());
+  Double_t p = - (n * a);
+  
+  // The center of the square with the fiducial marks as the
+  // corners.  The mid-point of one diagonal - md.  Used to get the
+  // center of the surveyd box. 
+  TVector3 md(a + c);
+  md *= 1/2.;
+  
+  // The center of the box. 
+  TVector3 orig(md);
+  trans[0] = orig[0];
+  trans[1] = orig[1];
+  trans[2] = orig[2];
+  
+  // Normalize the spanning vectors 
+  TVector3 uab(ab.Unit());
+  TVector3 ubc(bc.Unit());
+  
+  for (size_t i = 0; i < 3; i++) { 
+    rot[i * 3 + 0] = uab[i];
+    rot[i * 3 + 1] = ubc[i];
+    rot[i * 3 + 2] = n[i];
+  }
+  return kTRUE;
+}
+
+
+//____________________________________________________________________
+Bool_t
+AliFMDSurveyToAlignObjs::GetFMD1Plane(Double_t* rot, Double_t* trans) const
+{
+
+  // The possile survey points 
+  TVector3  icb, ict, ocb, oct, dummy;
+  Int_t     missing = 0;
+  if (!GetPoint("V0L_ICB", icb, dummy)) missing++;
+  if (!GetPoint("V0L_ICT", icb, dummy)) missing++;
+  if (!GetPoint("V0L_OCB", icb, dummy)) missing++;
+  if (!GetPoint("V0L_OCT", icb, dummy)) missing++;
+
+  // Check that we have enough points
+  if (missing > 1) { 
+    Error("GetFMD1Plane", "Only got %d survey points - no good for FMD1 plane",
+         4-missing);
+    return kFALSE;
+  }
+
+  if (!CalculatePlane(icb, ict, oct, rot, trans)) return kFALSE;
+
+  return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDSurveyToAlignObjs::DoFMD1()
+{
+  return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, Double_t* trans) const
+{
+
+  // The possile survey points 
+  const char*    names[] = { "FMD2_ITOP",  "FMD2_OTOP", 
+                            "FMD2_IBOTM", "FMD2_OBOTM", 
+                            "FMD2_IBOT",  "FMD2_OBOT", 
+                            0 };
+  const char**   name    = names;
+  Int_t          i       = 0;
+  TGraph2DErrors g;
   // Loop and fill graph 
   while (*name) {
-    TObject* obj = fSurveyPoints->FindObject(*name);
-    name++;
-    if (!obj) continue;
-    
-    AliSurveyPoint* p = static_cast<AliSurveyPoint*>(obj);
-    Double_t x  = unit * p->GetX();
-    Double_t y  = unit * p->GetY();
-    Double_t z  = unit * p->GetZ();
-    Double_t ex = unit * p->GetPrecisionX();
-    Double_t ey = unit * p->GetPrecisionY();
-    Double_t ez = unit * p->GetPrecisionZ();
-    
-    g.SetPoint(i, x, y, z);
-    g.SetPointError(i, ex, ey, ez);
+    TVector3 p, e;
+    if (!GetPoint(*name++, p, e)) continue;
+
+    g.SetPoint(i, p.X(), p.Y(), p.Z());
+    g.SetPointError(i, e.X(), e.Y(), e.Z());
     i++;
   }
   
@@ -83,9 +197,9 @@ AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot,
   TVector3 ua(a.Unit());
   TVector3 ub(b.Unit());
   Double_t angAb = ua.Angle(ub);
-  PrintVector("ua: ", ua);
-  PrintVector("ub: ", ub);
-  std::cout << "Angle: " << angAb * 180 / TMath::Pi() << std::endl;
+  // PrintVector("ua: ", ua);
+  // PrintVector("ub: ", ub);
+  // std::cout << "Angle: " << angAb * 180 / TMath::Pi() << std::endl;
   
   for (size_t i = 0; i < 3; i++) { 
     rot[i * 3 + 0] = ua[i];
index e510dbdfb28bc8cea70bc4c8bb176aa8cb55e00a..02b793cf47a0850c574858009010a6431beff2ee 100644 (file)
@@ -13,9 +13,18 @@ public:
   void Run();
   Bool_t CreateAlignObjs() { return kTRUE; }
 protected:
+  Bool_t DoFMD1();
+  Bool_t GetFMD1Plane(Double_t* rot, Double_t* trans) const;
   Bool_t DoFMD2();
-  Bool_t GetFMD2Plane(Double_t* rot, Double_t* trans);
-
+  Bool_t GetFMD2Plane(Double_t* rot, Double_t* trans) const;
+
+  Double_t GetUnitFactor() const;
+  Bool_t   GetPoint(const char* name, TVector3& p, TVector3& e) const;
+  Bool_t   CalculatePlane(const     TVector3& a, 
+                         const     TVector3& b,
+                         const     TVector3& c, 
+                         Double_t* trans,
+                         Double_t* rot) const;
   static void PrintVector(const char* text, const Double_t* v);
   static void PrintVector(const char* text, const TVector3& v);
   static void PrintRotation(const char* text, const Double_t* rot);
index a5bc57bb915709086b83cb6bf0b24a03e3c208cb..8e4105e328ededd1509c2e4148d63e38fff1d5a3 100644 (file)
@@ -170,7 +170,7 @@ enum MC_t {
 // Functions
 Float_t       EtaToTheta(Float_t eta);
 EG_t          LookupEG(const Char_t* name);
-AliGenerator* GeneratorFactory(EG_t eg, Rad_t rad, TString& comment);
+AliGenerator* GeneratorFactory(EG_t eg, Rad_t rad);
 AliGenHijing* HijingStandard();
 void          ProcessEnvironmentVars(EG_t& eg, Int_t& seed);
 
@@ -192,10 +192,6 @@ Config()
   Int_t seed = 12345; //Set 0 to use the current time
   MC_t  mc   = kGEANT3TGEO;
   
-  //____________________________________________________________________
-  // Comment line 
-  static TString  comment;
-  
   //____________________________________________________________________
   // Get settings from environment variables
   ProcessEnvironmentVars(eg, seed);
@@ -367,7 +363,7 @@ Config()
   
   //__________________________________________________________________
   // Generator Configuration
-  AliGenerator* gener = GeneratorFactory(eg, rad, comment);
+  AliGenerator* gener = GeneratorFactory(eg, rad);
   gener->SetOrigin(0, 0, 0);    // vertex position
   gener->SetSigma(0, 0, 5.3);   // Sigma in (X,Y,Z) (cm) on IP position
   gener->SetCutVertexZ(1.);     // Truncate at 1 sigma
@@ -375,30 +371,6 @@ Config()
   gener->SetTrackingFlag(1);
   gener->Init();
     
-  //__________________________________________________________________
-  // 
-  // Comments 
-  // 
-  switch (mag) {
-  case k2kG: comment = comment.Append(" | L3 field 0.2 T"); break;
-  case k4kG: comment = comment.Append(" | L3 field 0.4 T"); break;
-  case k5kG: comment = comment.Append(" | L3 field 0.5 T"); break;
-  }
-
-  switch (rad) {
-  case kGluonRadiation: 
-    comment = comment.Append(" | Gluon Radiation On");  break;
-  default:
-    comment = comment.Append(" | Gluon Radiation Off"); break;
-  }
-
-  switch(geo) {
-  case kHoles: comment = comment.Append(" | Holes for PHOS/HMPID"); break;
-  default:     comment = comment.Append(" | No holes for PHOS/HMPID"); break;
-  }
-
-  std::cout << "\n\n Comment: " << comment << "\n" << std::endl;
-
   //__________________________________________________________________
   // Field (L3 0.4 T)
   AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., mag);
@@ -488,7 +460,7 @@ LookupEG(const Char_t* name)
 
 //____________________________________________________________________  
 AliGenerator* 
-GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)  
+GeneratorFactory(EG_t eg, Rad_t rad)  
 {
   Int_t isw = 3;
   if (rad == kNoGluonRadiation) isw = 0;
@@ -520,7 +492,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
   switch (eg) {
   case test50:
     {
-      comment = comment.Append(":HIJINGparam test 50 particles");
       AliGenHIJINGpara *gener = new AliGenHIJINGpara(50);
       gener->SetMomentumRange(0, 999999.);
       gener->SetPhiRange(0., 360.);
@@ -533,7 +504,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kParam_8000:
     {
-      comment = comment.Append(":HIJINGparam N=8000");
       AliGenHIJINGpara *gener = new AliGenHIJINGpara(86030);
       gener->SetMomentumRange(0, 999999.);
       gener->SetPhiRange(0., 360.);
@@ -546,7 +516,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kParam_4000:
     {
-      comment = comment.Append("HIJINGparam N=4000");
       AliGenHIJINGpara *gener = new AliGenHIJINGpara(43015);
       gener->SetMomentumRange(0, 999999.);
       gener->SetPhiRange(0., 360.);
@@ -559,7 +528,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kParam_2000:
     {
-      comment = comment.Append("HIJINGparam N=2000");
       AliGenHIJINGpara *gener = new AliGenHIJINGpara(21507);
       gener->SetMomentumRange(0, 999999.);
       gener->SetPhiRange(0., 360.);
@@ -572,7 +540,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kParam_fmd:
     {
-      comment = comment.Append("HIJINGparam N=100");
       AliGenHIJINGpara *gener = new AliGenHIJINGpara(500);
       gener->SetMomentumRange(0, 999999.);
       gener->SetPhiRange(0., 360.);
@@ -588,7 +555,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     //
   case kHijing_cent1:
     {
-      comment = comment.Append("HIJING cent1");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -597,7 +563,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kHijing_cent2:
     {
-      comment = comment.Append("HIJING cent2");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 2.);
@@ -609,7 +574,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     //
   case kHijing_per1:
     {
-      comment = comment.Append("HIJING per1");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(5., 8.6);
@@ -618,7 +582,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kHijing_per2:
     {
-      comment = comment.Append("HIJING per2");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(8.6, 11.2);
@@ -627,7 +590,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kHijing_per3:
     {
-      comment = comment.Append("HIJING per3");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(11.2, 13.2);
@@ -636,7 +598,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kHijing_per4:
     {
-      comment = comment.Append("HIJING per4");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(13.2, 15.);
@@ -645,7 +606,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kHijing_per5:
     {
-      comment = comment.Append("HIJING per5");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(15., 100.);
@@ -657,7 +617,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     //
   case kHijing_jj25:
     {
-      comment = comment.Append("HIJING Jet 25 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -674,7 +633,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
 
   case kHijing_jj50:
     {
-      comment = comment.Append("HIJING Jet 50 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -691,7 +649,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
 
   case kHijing_jj75:
     {
-      comment = comment.Append("HIJING Jet 75 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -708,7 +665,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
 
   case kHijing_jj100:
     {
-      comment = comment.Append("HIJING Jet 100 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -725,7 +681,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
 
   case kHijing_jj200:
     {
-      comment = comment.Append("HIJING Jet 200 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -744,7 +699,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     //
   case kHijing_gj25:
     {
-      comment = comment.Append("HIJING Gamma 25 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -761,7 +715,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
 
   case kHijing_gj50:
     {
-      comment = comment.Append("HIJING Gamma 50 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -778,7 +731,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
 
   case kHijing_gj75:
     {
-      comment = comment.Append("HIJING Gamma 75 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -795,7 +747,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
 
   case kHijing_gj100:
     {
-      comment = comment.Append("HIJING Gamma 100 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -812,7 +763,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
 
   case kHijing_gj200:
     {
-      comment = comment.Append("HIJING Gamma 200 GeV");
       AliGenHijing *gener = HijingStandard();
       // impact parameter range
       gener->SetImpactParameterRange(0., 5.);
@@ -828,7 +778,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kHijing_pA:
     {
-      comment = comment.Append("HIJING pA");
 
       AliGenCocktail *gener  = new AliGenCocktail();
 
@@ -865,7 +814,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6:
     {
-      comment = comment.Append(":Pythia p-p @ 14 TeV");
       AliGenPythia *gener = new AliGenPythia(-1); 
       gener->SetMomentumRange(0,999999);
       gener->SetThetaRange(0., 180.);
@@ -878,7 +826,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets20_24:
     {
-      comment = comment.Append(":Pythia jets 20-24 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -897,7 +844,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets24_29:
     {
-      comment = comment.Append(":Pythia jets 24-29 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -916,7 +862,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets29_35:
     {
-      comment = comment.Append(":Pythia jets 29-35 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -935,7 +880,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets35_42:
     {
-      comment = comment.Append(":Pythia jets 35-42 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -954,7 +898,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets42_50:
     {
-      comment = comment.Append(":Pythia jets 42-50 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -973,7 +916,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets50_60:
     {
-      comment = comment.Append(":Pythia jets 50-60 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -992,7 +934,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets60_72:
     {
-      comment = comment.Append(":Pythia jets 60-72 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -1011,7 +952,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets72_86:
     {
-      comment = comment.Append(":Pythia jets 72-86 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -1030,7 +970,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets86_104:
     {
-      comment = comment.Append(":Pythia jets 86-104 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -1049,7 +988,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets104_125:
     {
-      comment = comment.Append(":Pythia jets 105-125 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -1068,7 +1006,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets125_150:
     {
-      comment = comment.Append(":Pythia jets 125-150 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -1087,7 +1024,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPythia6Jets150_180:
     {
-      comment = comment.Append(":Pythia jets 150-180 GeV @ 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);//        Centre of mass energy
       gener->SetProcess(kPyJets);//        Process type
@@ -1106,7 +1042,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kD0PbPb5500:
     {
-      comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(10);
       gener->SetProcess(kPyD0PbPbMNR);
       gener->SetStrucFunc(kCTEQ4L);
@@ -1123,7 +1058,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kCharmSemiElPbPb5500:
     {
-      comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV");
       AliGenPythia * gener = new AliGenPythia(10);
       gener->SetProcess(kPyCharmPbPbMNR);
       gener->SetStrucFunc(kCTEQ4L);
@@ -1139,7 +1073,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kBeautySemiElPbPb5500:
     {
-      comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV");
       AliGenPythia *gener = new AliGenPythia(10);
       gener->SetProcess(kPyBeautyPbPbMNR);
       gener->SetStrucFunc(kCTEQ4L);
@@ -1155,7 +1088,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kCocktailTRD:
     {
-      comment = comment.Append(" Cocktail for TRD at 5.5 TeV");
       AliGenCocktail *gener  = new AliGenCocktail();
 
       AliGenParam *jpsi = new AliGenParam(10,
@@ -1201,7 +1133,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPyJJ:
     {
-      comment = comment.Append(" Jet-jet at 5.5 TeV");
       AliGenPythia *gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);
       gener->SetProcess(kPyJets);
@@ -1215,7 +1146,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kPyGJ:
     {
-      comment = comment.Append(" Gamma-jet at 5.5 TeV");
       AliGenPythia *gener = new AliGenPythia(-1);
       gener->SetEnergyCMS(5500.);
       gener->SetProcess(kPyDirectGamma);
@@ -1230,7 +1160,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailCent1:
     {
-      comment = comment.Append(" Muon Cocktail Cent1");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1245,7 +1174,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailPer1:
     {
-      comment = comment.Append(" Muon Cocktail Per1");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1260,7 +1188,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailPer4:
     {
-      comment = comment.Append(" Muon Cocktail Per4");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1275,7 +1202,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailCent1HighPt:
     {
-      comment = comment.Append(" Muon Cocktail HighPt Cent1");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1290,7 +1216,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailPer1HighPt :
     {
-      comment = comment.Append(" Muon Cocktail HighPt Per1");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1305,7 +1230,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailPer4HighPt:
     {
-      comment = comment.Append(" Muon Cocktail HighPt Per4");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1320,7 +1244,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailCent1Single:
     {
-      comment = comment.Append(" Muon Cocktail Single Cent1");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1335,7 +1258,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailPer1Single :
     {
-      comment = comment.Append(" Muon Cocktail Single Per1");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1350,7 +1272,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kMuonCocktailPer4Single:
     {
-      comment = comment.Append(" Muon Cocktail Single Per4");
       AliGenMUONCocktail * gener = new AliGenMUONCocktail();
       gener->SetPtRange(1.0,100.);       // Transverse momentum range   
       gener->SetPhiRange(0.,360.);    // Azimuthal angle range  
@@ -1365,7 +1286,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kFMD1Flat: 
     {
-      comment = comment.Append(" Flat in FMD1 range");
       AliGenBox* gener = new AliGenBox(2000);
       gener->SetPart(kPiPlus);
       gener->SetMomentumRange(3,4);
@@ -1376,7 +1296,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kFMD2Flat: 
     {
-      comment = comment.Append(" Flat in FMD2 range");
       AliGenBox* gener = new AliGenBox(10);
       gener->SetPart(kPiPlus);
       gener->SetMomentumRange(3,4);
@@ -1387,7 +1306,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kFMD3Flat: 
     {
-      comment = comment.Append(" Flat in FMD3 range");
       AliGenBox* gener = new AliGenBox(2000);
       gener->SetPart(kPiPlus);
       gener->SetMomentumRange(3,4);
@@ -1398,7 +1316,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
     break;
   case kFMDFlat:
     {
-      comment = comment.Append(" Flat in FMD range");
       AliGenCocktail* gener = new AliGenCocktail();
       gener->SetMomentumRange(3,4);
       gener->SetPhiRange(0, 360);
index a1a051bb4c241d580588f318291e1dc3fe023a96..5a73e1f4155291d9dff3a81af19bf1ff4bdaa0cf 100644 (file)
@@ -26,16 +26,16 @@ void
 Reconstruct()
 {
   // Debug the FMD.
-  // AliLog::SetModuleDebugLevel("FMD", 1);
+  AliLog::SetModuleDebugLevel("FMD", 1);
 
   // To reconstruct raw data from FDR-I, please enable below lines: 
   // AliFMDParameters::Instance()->UseRcuTrailer(false);
   // AliFMDParameters::Instance()->UseCompleteHeader(false);
  
-  TFile* magF = TFile::Open("mag.root", "READ");
-  AliMagF* mag = static_cast<AliMagF*>(magF->Get("mag"));
-  if (!mag) return;
-  AliTracker::SetFieldMap(mag, true);
+  // TFile* magF = TFile::Open("mag.root", "READ");
+  // AliMagF* mag = static_cast<AliMagF*>(magF->Get("mag"));
+  // if (!mag) return;
+  // AliTracker::SetFieldMap(mag, true);
   AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
   AliCDBManager::Instance()->SetRun(0);
 
@@ -45,8 +45,8 @@ Reconstruct()
   rec.SetRunReconstruction("FMD");
   rec.SetRunTracking(""); 
   rec.SetFillESD("FMD"); 
-  rec.SetRunQA("");
-  rec.SetInput(".");
+  rec.SetRunQA(":");
+  rec.SetInput("raw.root");
   // rec.SetRecoParam("TOF", new AliTOFRecoParam());
   
   rec.Run(); 
index da627ef8cdad80a03d72d42f7a18bbb92cf65ec6..7415b8b5986f459a19a6506ca71c65d23b3ab34b 100644 (file)
@@ -28,9 +28,9 @@ Simulate(Int_t n=1)
   sim.SetConfigFile("$(ALICE_ROOT)/FMD/Config.C");
   sim.SetMakeSDigits("FMD");
   sim.SetMakeDigits("FMD"); 
-  sim.SetWriteRawData("FMD"); 
+  sim.SetWriteRawData("FMD", "raw.root"); 
   // sim.SetMakeDigitsFromHits("FMD"); 
-  sim.SetRunQA(":");
+  sim.SetRunQA("FMD:ALL");
   TStopwatch w; 
   w.Start(); 
   sim.Run(n);  
diff --git a/FMD/Survey_943928_FMD.txt b/FMD/Survey_943928_FMD.txt
new file mode 100644 (file)
index 0000000..4cb8b95
--- /dev/null
@@ -0,0 +1,49 @@
+> Title:
+
+Measurement of targets on FMD2
+
+> Date:
+
+12.08.2008
+
+> Subdetector:
+
+FMD1
+
+> Report URL:
+
+https://edms.cern.ch/document/943928
+
+> Version:
+
+1
+
+> General Observations:
+
+The V0A detector has been installed and first measured on 18.07.2008.
+Afterwards it has been moved close to final position and measured on
+25.07.2008.  The final alignment of the V0A detector has been done in
+three iterations (called sessions) on 12.08.2008.
+
+> Coordinate System:
+
+ALICEPH
+
+> Units:
+
+mm
+
+> Nr Columns:
+
+7
+
+> Column Names:
+
+Point Name,XPH,YPH,ZPH,Point Type,Target Used,Precision(mm)
+
+
+> Data:
+V0L_ICB                127.73  -220.41 3246.58 M       Y       1
+V0L_ICT                127.37  219.48  3247.24 M       Y       1
+V0L_OCB                -126.99 -220.40 3247.27 M       Y       1
+V0L_OCT                -126.99 219.51  3247.49 M       Y       1
diff --git a/FMD/Survey_976326_FMD.txt b/FMD/Survey_976326_FMD.txt
new file mode 100644 (file)
index 0000000..618f247
--- /dev/null
@@ -0,0 +1,51 @@
+> Title:
+
+Measurement of targets on FMD2
+
+> Date:
+
+05.06.2008
+
+> Subdetector:
+
+FMD2
+
+> Report URL:
+
+https://edms.cern.ch/document/976326
+
+> Version:
+
+1
+
+> General Observations:
+
+During the installation period of detectors in the Space Frame,
+several targets on the FMD2 have been measured regularly in order to
+detect relative movements of the TPC.
+
+
+> Coordinate System:
+
+ALICEPH
+
+> Units:
+
+mm
+
+> Nr Columns:
+
+7
+
+> Column Names:
+
+Point Name,XPH,YPH,ZPH,Point Type,Target Used,Precision(mm)
+
+
+> Data:
+
+FMD2_IBOT      192.7   -257.4  892.4   M       Y       1
+FMD2_OBOT      -194.5  -257.2  894.3   M       Y       1
+FMD2_OBOTM     -271.7  -122.0  892.1   M       Y       1
+FMD2_ITOP      246.9   100.4   892.5   M       Y       1
+FMD2_OTOP      -232.4  113.6   891.7   M       Y       1
diff --git a/FMD/scripts/DrawBothDigits.C b/FMD/scripts/DrawBothDigits.C
new file mode 100644 (file)
index 0000000..3a018e5
--- /dev/null
@@ -0,0 +1,110 @@
+//____________________________________________________________________
+//
+// $Id: DrawBothDigits.C 20907 2007-09-25 08:44:03Z cholm $
+//
+// Script that contains a class to draw eloss from hits, versus ADC
+// counts from digits, using the AliFMDInputHits class in the util library. 
+//
+// It draws the energy loss versus the p/(mq^2).  It can be overlayed
+// with the Bethe-Bloc curve to show how the simulation behaves
+// relative to the expected. 
+//
+// Use the script `Compile.C' to compile this class using ACLic. 
+//
+#include <TH2F.h>
+#include <AliFMDDigit.h>
+#include <AliFMDSDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDEdepMap.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <AliLog.h>
+
+/** @class DrawBothDigits
+    @brief Draw hit energy loss versus digit ADC
+    @code 
+    Root> .L Compile.C
+    Root> Compile("DrawBothDigits.C")
+    Root> DrawBothDigits c
+    Root> c.Run();
+    @endcode
+    @ingroup FMD_script
+ */
+class DrawBothDigits : public AliFMDInput
+{
+private:
+  TH2F* fTrackNos; // Histogram 
+  AliFMDEdepMap fCache;
+public:
+  //__________________________________________________________________
+  DrawBothDigits(Int_t max=300) 
+    : AliFMDInput("galice.root")
+  { 
+    AddLoad(kDigits);
+    AddLoad(kSDigits);
+    fTrackNos = new TH2F("trackNos", "Track numbers", 
+                        max+1, -1.5, max-.5, max+1, -1.5, max-.5);
+    fTrackNos->SetXTitle("Digit track");
+    fTrackNos->SetYTitle("SDigit track");
+  }
+  //__________________________________________________________________
+  Bool_t Begin(Int_t evno)
+  {
+    fCache.Reset();
+    return AliFMDInput::Begin(evno);
+  }
+  //__________________________________________________________________
+  Bool_t ProcessSDigit(AliFMDSDigit* sdigit)
+  {
+    if (!sdigit) return kTRUE;
+    AliFMDEdepHitPair& entry = fCache(sdigit->Detector(), 
+                                     sdigit->Ring(), 
+                                     sdigit->Sector(), 
+                                     sdigit->Strip());
+    entry.fLabels.Set(sdigit->GetNTrack());
+    Info("ProcessSDigit", "Got %d SDigit tracks", sdigit->GetNTrack());
+    for (size_t i = 0; i < sdigit->GetNTrack(); i++) 
+      entry.fLabels.fArray[i] = sdigit->GetTrack(i);
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t ProcessDigit(AliFMDDigit* digit)
+  {
+    if (!digit) return kTRUE;
+    AliFMDEdepHitPair& entry = fCache(digit->Detector(), 
+                                     digit->Ring(), 
+                                     digit->Sector(), 
+                                     digit->Strip());
+    TArrayI& stracks = entry.fLabels;
+    Info("ProcessDigit", "Got %d SDigit tracks, and %d Digit tracks", 
+        stracks.fN, digit->GetNTrack());
+    for (Int_t i = 0; (i < stracks.fN || i < digit->GetNTrack()); i++) { 
+      Int_t strack = (i < stracks.fN ? stracks.fArray[i] : -1);
+      Int_t dtrack = (i < digit->GetNTrack() ? digit->GetTrack(i) : -1);
+      fTrackNos->Fill(strack, dtrack);
+    }
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t Finish()
+  {
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+    fTrackNos->SetStats(kFALSE);
+    // fTrackNos->Scale(1. / fAdc->GetEntries());
+    fTrackNos->Draw("colz");
+    return kTRUE;
+  }
+
+  ClassDef(DrawBothDigits,0);
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
index 107ba89cf10896ee3979a54b4c6c0d5b13c9380f..62ddea5b4d5b23684bf58ac4a35ed107ee75b4ef 100644 (file)
@@ -38,6 +38,7 @@ private:
 public:
   //__________________________________________________________________
   DrawDigits(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1023.5) 
+    : AliFMDInput("galice.root")
   { 
     AddLoad(kDigits);
     fAdc = new TH1D("adc", "ADC", m, amin, amax);
@@ -49,6 +50,7 @@ public:
   {
     if (!digit) return kTRUE;
     fAdc->Fill(digit->Counts());
+    digit->Print("l");
     return kTRUE;
   }
   //__________________________________________________________________
index e5e2b7cf6b8bd64966d4f1e6af3a280a598b8eb5..94924fca07b6443f78e4bc28963aa52940aee0d5 100644 (file)
@@ -64,20 +64,24 @@ public:
     return AliFMDInput::Begin(ev);
   }
   //__________________________________________________________________
-  Bool_t ProcessESD(UShort_t /* det */
-                   Char_t   /* ring */
-                   UShort_t /* sec */
-                   UShort_t /* strip */
+  Bool_t ProcessESD(UShort_t det
+                   Char_t   rng
+                   UShort_t sec
+                   UShort_t str
                    Float_t  /* eta */, 
                    Float_t  mult)
   {
     // Cache the energy loss 
-    if (mult/fCorr > 0.01) fMult->Fill(mult/fCorr);
+    if (mult > 0) 
+      Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
+    if (mult/fCorr > 0.001) fMult->Fill(mult/fCorr);
     return kTRUE;
   }
   //__________________________________________________________________
   TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
   {
+    if (TMath::Abs(max-min) < .25) return 0;
+    std::cout << "Fit peack in range " << min << " to " << max << std::endl;
     TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
     hist->GetXaxis()->SetRangeUser(0, 4);
     hist->Fit(l, "0Q", "", min, max);
@@ -112,6 +116,7 @@ public:
   //__________________________________________________________________
   Bool_t Finish()
   {
+    Info("Finish", "Will draw results");
     gStyle->SetPalette(1);
     gStyle->SetOptTitle(0);
     gStyle->SetOptFit(1111111);
@@ -120,16 +125,20 @@ public:
     gStyle->SetPadColor(0);
     gStyle->SetPadBorderSize(0);
     
+    if (fMult->GetEntries() <= 0) return kFALSE;
+    
     TCanvas* c = new TCanvas("c", "C");
     c->cd();
-    c->SetLogy();
+    // c->SetLogy();
     fMult->GetXaxis()->SetRangeUser(0,4);
     fMult->Scale(1. / fMult->GetEntries());
     fMult->SetStats(kFALSE);
     fMult->SetFillColor(2);
     fMult->SetFillStyle(3001);
     fMult->Draw("hist e");
-
+    
+    return kTRUE;
+    
     Double_t mean, rms;
     MaxInRange(fMult, 0.2, mean, rms);
     Double_t x1   = mean-rms/2; // .75; // .8;  // .65 / fCorr;
@@ -139,49 +148,61 @@ public:
     MaxInRange(fMult, x2, mean, rms);
     Double_t x3   = mean + rms;
     TF1*     l2   = FitPeak(2, fMult, x2, x3);
-    Double_t diff = l2->GetParameter(1)-l1->GetParameter(1);
-    TF1*     f    = new TF1("user", "landau(0)+landau(3)", x1, x3);
-
+    TF1*     f    = 0;
+    Double_t diff = 0;
+    if (l2) {
+      diff          = l2->GetParameter(1)-l1->GetParameter(1);
+      f             = new TF1("user", "landau(0)+landau(3)", x1, x3);
+      f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
+                    "A_{2}", "Mpv_{2}", "#sigma_{2}",
+                    "A_{3}", "Mpv_{3}", "#sigma_{3}");
+      f->SetParameters(l1->GetParameter(0), 
+                      l1->GetParameter(1), 
+                      l1->GetParameter(2), 
+                      l2 ? l2->GetParameter(0) : 0, 
+                      l2 ? l2->GetParameter(1) : 0, 
+                      l2 ? l2->GetParameter(2) : 0,
+                      l2->GetParameter(0)/10, 
+                      l2->GetParameter(1) + diff, 
+                      l2->GetParameter(2));
+    }
+    else { 
+      x3 = x2;
+      f  = new TF1("usr", "landau", x1, x3);
+    }
+    
+    std::cout << "Range: " << x1 << "-" << x3 << std::endl;
+    
     fMult->GetXaxis()->SetRangeUser(0, 4);
-    f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
-                  "A_{2}", "Mpv_{2}", "#sigma_{2}",
-                  "A_{3}", "Mpv_{3}", "#sigma_{3}");
-    f->SetParameters(l1->GetParameter(0), 
-                    l1->GetParameter(1), 
-                    l1->GetParameter(2), 
-                    l2->GetParameter(0), 
-                    l2->GetParameter(1), 
-                    l2->GetParameter(2),
-                    l2->GetParameter(0)/10, 
-                    l2->GetParameter(1) + diff, 
-                    l2->GetParameter(2));
-    fMult->Fit(f, "0Q", "", x1, x3);
-    fMult->Fit(f, "ME0", "E1", x1, x3);
+    fMult->Fit(f, "0QR", "", x1, x3);
+    fMult->Fit(f, "ME0R", "E1", x1, x3);
     fMult->DrawClone("same hist");
 
     l1->SetLineColor(3);
     l1->SetLineWidth(2);
     l1->SetRange(0, 4);
     l1->Draw("same");
-    l2->SetLineColor(4);
-    l2->SetLineWidth(2);
-    l2->SetRange(0, 4);
-    l2->Draw("same");
+    if (l2) {
+      l2->SetLineColor(4);
+      l2->SetLineWidth(2);
+      l2->SetRange(0, 4);
+      l2->Draw("same");
+    }
     f->SetLineWidth(2);
     f->SetRange(0, 4);
     f->Draw("same");
 
     TLegend* l = new TLegend(0.6, 0.6, .89, .89);
     l->AddEntry(l1, "1 particle Landau", "l");
-    l->AddEntry(l2, "2 particle Landau", "l");
+    if (l2) l->AddEntry(l2, "2 particle Landau", "l");
     l->AddEntry(f,  "1+2 particle Landau", "l");
     l->SetFillColor(0);
     l->Draw("same");
 
 
-#if 1
+#if 0
     c = new TCanvas("c2", "Landaus");
-    c->SetLogy();
+    // c->SetLogy();
     fMult->DrawClone("axis");
     f->Draw("same");
     Double_t* p1 = f->GetParameters();
index 2d6fabb3c3e813eb795c3fe9f7781628d91d9844..c314038cf4e16bf3b9170474a0d3dde9dfc3e02d 100644 (file)
@@ -90,7 +90,7 @@ public:
       std::cout << "No track" << std::endl;
       return kFALSE;
     }
-    if (!p->IsPrimary()) return kTRUE;
+    // if (!p->IsPrimary()) return kTRUE;
     if (hit->IsStop()) return kTRUE;
 
     Float_t x = hit->P();
@@ -103,7 +103,7 @@ public:
     y /= q * q;
     fElossVsPMQ->Fill(x, y);
     fEloss->Fill(y);
-
+    fNHits++;
     return kTRUE;
   }
   //__________________________________________________________________
@@ -155,24 +155,27 @@ public:
 
     c = new TCanvas("eloss", "Energy loss per unit material");
     // c->SetLogx();
-    c->SetLogy();
-    fEloss->Scale(1. / fEloss->GetEntries());
-    fEloss->GetXaxis()->SetRangeUser(1, 10);
-    fEloss->Fit("landau", "", "", 1, 10);
-    TF1* land = fEloss->GetFunction("landau");
-    land->SetLineWidth(2);
-    Double_t max = fEloss->GetMaximum();
-    TGraph* resp  = GetResp();
-    Double_t* x   = resp->GetX();
-    Double_t* y   = resp->GetY();
-    TGraph*   g   = new TGraph(resp->GetN());
-    TGraph*   co  = GetCorr();
-    std::cout << "Correction factor: " << co->Eval(fBetaGammaMip) << std::endl;
-    Double_t  xs  = fRho; // * 1.19; // / 
-    for (Int_t i = 0; i < g->GetN(); i++) 
-      g->SetPoint(i, x[i] * xs, y[i] * max);
-    g->Draw("C same");
-
+    if (fEloss->GetEntries() != 0) { 
+      c->SetLogy();
+      fEloss->Scale(1. / fEloss->GetEntries());
+      fEloss->GetXaxis()->SetRangeUser(1, 10);
+      fEloss->Fit("landau", "", "", 1, 10);
+      TF1* land = fEloss->GetFunction("landau");
+      land->SetLineWidth(2);
+      Double_t max = fEloss->GetMaximum();
+      TGraph* resp  = GetResp();
+      Double_t* x   = resp->GetX();
+      Double_t* y   = resp->GetY();
+      TGraph*   g   = new TGraph(resp->GetN());
+      TGraph*   co  = GetCorr();
+      std::cout << "Correction factor: " << co->Eval(fBetaGammaMip) 
+               << std::endl;
+      Double_t  xs  = fRho; // * 1.19; // / 
+      for (Int_t i = 0; i < g->GetN(); i++) 
+       g->SetPoint(i, x[i] * xs, y[i] * max);
+      g->Draw("C same");
+    }
+    
     l = new TLegend(.6, .6, .89, .89);
     l->AddEntry(fEloss, fEloss->GetTitle(), "lf");
     l->AddEntry(land,   "Landau fit", "l");
index 94c8df276855884e63b44553f3cea51ff3b37eee..e5bd31b88065563e88ff83d50be6140faeca6129 100644 (file)
@@ -65,6 +65,7 @@ public:
   //__________________________________________________________________
   Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
   {
+    // Info("ProcessHit", "Processing hit %p", hit);
     // Cache the energy loss 
     if (!hit) return kFALSE;
     UShort_t det = hit->Detector();
@@ -77,11 +78,13 @@ public:
     }
     fMap(det, rng, sec, str).fEdep += hit->Edep();
     fMap(det, rng, sec, str).fN++;
+    // hit->Print();
     return kTRUE;
   }
   //__________________________________________________________________
   Bool_t ProcessDigit(AliFMDDigit* digit)
   {
+    // Info("ProcessDigit", "Processing digit %p", digit);
     if (!digit) return kFALSE;
     UShort_t det = digit->Detector();
     Char_t   rng = digit->Ring();
@@ -92,6 +95,10 @@ public:
       return kFALSE;
     }
     fElossVsAdc->Fill(fMap(det, rng, sec, str).fEdep, digit->Counts());
+    if (fMap(det, rng, sec, str).fEdep > 0) 
+      Info("ProcessDigit", "FMD%d%c[%2d,%3d] Edep=%8.5f -> ADC=0x%03x",
+          det, rng, sec, str, 
+          fMap(det, rng, sec, str).fEdep, digit->Counts());
     return kTRUE;
   }
   //__________________________________________________________________
diff --git a/FMD/scripts/DrawHitsSDigits.C b/FMD/scripts/DrawHitsSDigits.C
new file mode 100644 (file)
index 0000000..6a51ec0
--- /dev/null
@@ -0,0 +1,121 @@
+//____________________________________________________________________
+//
+// $Id: DrawHitsSDigits.C 20907 2007-09-25 08:44:03Z cholm $
+//
+// Script that contains a class to draw eloss from hits, versus ADC
+// counts from sdigits, using the AliFMDInputHits class in the util library. 
+//
+// It draws the energy loss versus the p/(mq^2).  It can be overlayed
+// with the Bethe-Bloc curve to show how the simulation behaves
+// relative to the expected. 
+//
+// Use the script `Compile.C' to compile this class using ACLic. 
+//
+#include <TH2D.h>
+#include <AliFMDHit.h>
+#include <AliFMDSDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDEdepMap.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <AliLog.h>
+#include <TMath.h>
+
+/** @class DrawHitsSDigits
+    @brief Draw hit energy loss versus sdigit ADC
+    @code 
+    Root> .L Compile.C
+    Root> Compile("DrawHitsSDigits.C")
+    Root> DrawHitsSDigits c
+    Root> c.Run();
+    @endcode
+    @ingroup FMD_script
+ */
+class DrawHitsSDigits : public AliFMDInput
+{
+private:
+  TH2D* fElossVsAdc; // Histogram 
+  AliFMDEdepMap fMap;
+public:
+  //__________________________________________________________________
+  DrawHitsSDigits(Int_t n=900, Double_t emin=1e-3, Double_t emax=10, 
+                Int_t m=1100, Double_t amin=-0.5, Double_t amax=1099.5) 
+  { 
+    AddLoad(kSDigits);
+    AddLoad(kHits);
+    TArrayF eloss(MakeLogScale(n, emin, emax));
+    TArrayF adcs(m+1);
+    adcs[0] = amin;
+    for (Int_t i = 1; i < m+1; i++) adcs[i] = adcs[i-1] + (amax-amin)/m;
+    fElossVsAdc = new TH2D("bad", "#Delta E vs. ADC", 
+                          eloss.fN-1, eloss.fArray, adcs.fN-1, adcs.fArray);
+    fElossVsAdc->SetXTitle("#Delta E/#Delta x [MeV/cm]");
+    fElossVsAdc->SetYTitle("ADC value");
+  }
+  //__________________________________________________________________
+  /** Begining of event
+      @param ev Event number
+      @return @c false on error */
+  Bool_t Begin(Int_t ev) 
+  {
+    fMap.Reset();
+    return AliFMDInput::Begin(ev);
+  }
+  //__________________________________________________________________
+  Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
+  {
+    Info("ProcessHit", "Processing hit %p", hit);
+    // Cache the energy loss 
+    if (!hit) return kFALSE;
+    UShort_t det = hit->Detector();
+    Char_t   rng = hit->Ring();
+    UShort_t sec = hit->Sector();
+    UShort_t str = hit->Strip();
+    if (str > 511) {
+      AliWarning(Form("Bad strip number %d in hit", str));
+      return kTRUE;
+    }
+    fMap(det, rng, sec, str).fEdep += hit->Edep();
+    fMap(det, rng, sec, str).fN++;
+    hit->Print();
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t ProcessSDigit(AliFMDSDigit* sdigit)
+  {
+    Info("ProcessSDigit", "Processing sdigit %p", sdigit);
+    if (!sdigit) return kFALSE;
+    UShort_t det = sdigit->Detector();
+    Char_t   rng = sdigit->Ring();
+    UShort_t sec = sdigit->Sector();
+    UShort_t str = sdigit->Strip();
+    if (str > 511) {
+      AliWarning(Form("Bad strip number %d in sdigit", str));
+      return kFALSE;
+    }
+    fElossVsAdc->Fill(fMap(det, rng, sec, str).fEdep, sdigit->Counts());
+    sdigit->Print();
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t Finish()
+  {
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+    fElossVsAdc->SetStats(kFALSE);
+    fElossVsAdc->Draw("COLZ");
+    return kTRUE;
+  }
+
+  ClassDef(DrawHitsSDigits,0);
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
index 53d404a1e76b568446da7c537fe9875f6ced7665..9ba1a86a4f5cec81a29fe16065c3cccfcacb3d92 100644 (file)
@@ -68,6 +68,7 @@ public:
     
   //__________________________________________________________________
   DrawSDigits(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1023.5) 
+    : AliFMDInput("galice.root")
   { 
     AddLoad(kSDigits);
     AddLoad(kGeometry);
@@ -121,6 +122,7 @@ public:
   {
     if (!digit) return kTRUE;
     fAdc->Fill(digit->Counts());
+    digit->Print("lp");
     if (digit->NParticles() == 0) return kTRUE;
     
 
@@ -143,7 +145,6 @@ public:
     Double_t ratio = digit->NPrimaries() / digit->NParticles();
     if (phi < 0) phi += 360;
     fPrimRatio[primIdx]->Fill(eta, phi, ratio);
-    
     return kTRUE;
   }
   //__________________________________________________________________
diff --git a/FMD/scripts/TestRaw2SDigits.C b/FMD/scripts/TestRaw2SDigits.C
new file mode 100644 (file)
index 0000000..874c70c
--- /dev/null
@@ -0,0 +1,10 @@
+void
+TestRaw2SDigits()
+{
+  AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
+  AliCDBManager::Instance()->SetRun(0);
+  AliLog::SetModuleDebugLevel("FMD", 2);
+  AliSimulation sim;
+  sim.SetConfigFile("$ALICE_ROOT/FMD/Config.C");
+  sim.ConvertRaw2SDigits("raw.root", 0);
+}