In AliFMDReconstructor and AliFMDRevertexer, use AliFMDGeometry
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Sep 2010 14:33:49 +0000 (14:33 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Sep 2010 14:33:49 +0000 (14:33 +0000)
service function to calculate eta, phi, theta, and r.

In AliFMDReconstructor, AliFMD, AliFMDQADataMakerSim,
fix needless TClonesArray allocations, deletions, clears, etc.

FMD/AliFMD.cxx
FMD/AliFMDBoolMap.h
FMD/AliFMDESDRevertexer.cxx
FMD/AliFMDGeometry.cxx
FMD/AliFMDGeometry.h
FMD/AliFMDQADataMakerSim.cxx
FMD/AliFMDReconstructor.cxx
FMD/AliFMDReconstructor.h

index 437fcf5..6544446 100644 (file)
@@ -145,7 +145,7 @@ AliFMD::AliFMD()
   fHits        = 0;
   fDigits      = 0;
   fIshunt      = 0;
-  fBad         = new TClonesArray("AliFMDHit");
+  // fBad         = new TClonesArray("AliFMDHit");
 }
 
 //____________________________________________________________________
@@ -162,15 +162,15 @@ AliFMD::AliFMD(const char *name, const char *title)
   // Standard constructor for Forward Multiplicity Detector
   //
   AliFMDDebug(10, ("\tStandard CTOR"));
-  fBad         = new TClonesArray("AliFMDHit");
+  // fBad         = new TClonesArray("AliFMDHit");
   
   // Initialise Hit array
-  HitsArray();
-  gAlice->GetMCApp()->AddHitList(fHits);
+  // HitsArray();
+  // gAlice->GetMCApp()->AddHitList(fHits);
 
   // (S)Digits for the detectors disk
-  DigitsArray();
-  SDigitsArray();
+  // DigitsArray();
+  // SDigitsArray();
   
   // CHC: What is this?
   fIshunt = 0;
@@ -736,7 +736,9 @@ AliFMD::AddHitByFields(Int_t    track,
 
   AliMC *mcApplication = (AliMC*)gAlice->GetMCApp();
   
-  AliTrackReference* trackRef = AddTrackReference(mcApplication->GetCurrentTrackNumber(), AliTrackReference::kFMD); 
+  AliTrackReference* trackRef = 
+    AddTrackReference(mcApplication->GetCurrentTrackNumber(), 
+                     AliTrackReference::kFMD); 
   UInt_t stripId = AliFMDStripIndex::Pack(detector,ring,sector,strip);
   trackRef->SetUserId(stripId);
   
@@ -898,6 +900,8 @@ AliFMD::HitsArray()
   if (!fHits) { 
     fHits = new TClonesArray("AliFMDHit", 1000);
     fNhits = 0;
+    if (gAlice && gAlice->GetMCApp() && gAlice->GetMCApp()->GetHitLists()) 
+      gAlice->GetMCApp()->AddHitList(fHits);
   }
   return fHits;
 }
index abb9858..c2ebada 100644 (file)
@@ -74,6 +74,8 @@ public:
                                   Char_t   ring,
                                   UShort_t sec,
                                   UShort_t str) const;
+  Bool_t* Data() const { return fData; }
+  Int_t   Total() const { return fTotal; }
 protected:
   Int_t     MaxIndex()            const { return fTotal; }
   Bool_t    AtAsBool(Int_t idx)   const { return fData[idx]; }
index e993cbd..59f14d0 100644 (file)
@@ -21,6 +21,17 @@ AliFMDESDRevertexer::AliFMDESDRevertexer()
 Bool_t
 AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz) const
 {
+  // Recalculate the various quantities based on updated 
+  // primary vertex position. 
+  // 
+  // Parameters: 
+  //    fmdEsd    FMD ESD object 
+  //    vz        New vertex location (along the z-axis)
+  //
+  // Return:
+  //    true on success, false if there was an error during the 
+  //    recalculations.   Please inspect log output for details. 
+  // 
   if (!fmdEsd) return kFALSE;
   
   Bool_t ret = kTRUE;
@@ -94,17 +105,7 @@ AliFMDESDRevertexer::PhysicalCoordinates(UShort_t  det,
   Double_t x=0, y=0, z=0;
   geom->Detector2XYZ(det, rng, sec, str, x, y, z);
 
-  // Check that the conversion succeeded
-  if (x == 0 && y == 0 && z == 0) return kFALSE;
-  
-  // Correct for vertex offset. 
-  z     -= vz;
-  phi   =  TMath::ATan2(y, x);
-  r     =  TMath::Sqrt(y * y + x * x);
-  theta =  TMath::ATan2(r, z);
-  eta   = -TMath::Log(TMath::Tan(theta / 2));
-
-  return kTRUE;
+  return AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z-vz, r, eta, phi, theta);
 }
 
 
index 1201bd3..d7a2e62 100644 (file)
@@ -76,6 +76,7 @@
 #include <TGeoManager.h>
 #include <TGeoVolume.h>
 #include <TGeoNode.h>
+#include <TMath.h>
 static Int_t FindNodeDepth(const char* name, const char* volname);
 
 
@@ -359,6 +360,41 @@ AliFMDGeometry::XYZ2Detector(Double_t  x,
   return kFALSE;
 }
 
+//____________________________________________________________________
+Bool_t
+AliFMDGeometry::XYZ2REtaPhiTheta(Double_t  x,   Double_t y, 
+                                Double_t  z, 
+                                Double_t& r,   Double_t& eta, 
+                                Double_t& phi, Double_t& theta)
+{
+  
+  // Service function to convert Cartisean XYZ to r, eta, phi, and theta.   
+  // 
+  // Note, that the z input should be corrected for the vertex location 
+  // if needed.
+  //
+  //     x      Cartisean X coordinate
+  //     y      Cartisean Y coordinate 
+  //     z      Cartisean Z coordinate 
+  //     r      On return, the radius
+  //     eta    On return, the pseudo-rapidity
+  //     phi    On return, the azimuthal angle
+  //     theta  On return, the polar angle;
+  // 
+  // Return:
+  //     kFALSE in case of problems. 
+
+  if (x == 0 && y == 0 && z == 0) return kFALSE;
+  
+  // Correct for vertex offset. 
+  phi   =  TMath::ATan2(y, x);
+  r     =  TMath::Sqrt(y * y + x * x);
+  theta =  TMath::ATan2(r, z);
+  eta   = -TMath::Log(TMath::Tan(theta / 2));
+
+  return kTRUE;
+}
+
 
 //____________________________________________________________________
 void
index 35760ce..754667a 100644 (file)
@@ -210,6 +210,29 @@ public:
   virtual Bool_t Impact(const TParticle* particle) const;
   /** Declare alignable volumes */
   virtual void SetAlignableVolumes() const;
+
+
+  /** 
+   * Service function to convert Cartisean XYZ to r, eta, phi, and theta.   
+   *
+   * Note, that the z input should be corrected for the vertex location 
+   * if needed.
+   * 
+   * @param x      Cartisean X coordinate
+   * @param y      Cartisean Y coordinate 
+   * @param z      Cartisean Z coordinate 
+   * @param r      On return, the radius
+   * @param eta    On return, the pseudo-rapidity
+   * @param phi    On return, the azimuthal angle
+   * @param theta  On return, the polar angle;
+   *
+   * @return kTRUE on success, kFALSE in case of problems
+   */     
+  static Bool_t XYZ2REtaPhiTheta(Double_t  x,   Double_t y, 
+                                Double_t  z, 
+                                Double_t& r,   Double_t& eta, 
+                                Double_t& phi, Double_t& theta);
+
 protected:
   Bool_t        fIsInitialized; // Whether singleton is initalized
   AliFMDRing*  fInner;         // Inner ring geometry information
index 4e422b4..250e3f1 100644 (file)
@@ -57,14 +57,15 @@ AliFMDQADataMakerSim::AliFMDQADataMakerSim()
 AliFMDQADataMakerSim::AliFMDQADataMakerSim(const AliFMDQADataMakerSim& /*qadm*/) 
   : AliQADataMakerSim()
 {
-  //copy ctor 
-  
+  // copy ctor 
+  // 
   // Parameters: 
   //    qadm    Object to copy from
   
 }
 //_____________________________________________________________________
-AliFMDQADataMakerSim& AliFMDQADataMakerSim::operator = (const AliFMDQADataMakerSim& ) 
+AliFMDQADataMakerSim& 
+AliFMDQADataMakerSim::operator = (const AliFMDQADataMakerSim& ) 
 {
   
   return *this;
@@ -95,7 +96,9 @@ void AliFMDQADataMakerSim::InitSDigits()
   const Bool_t expert   = kTRUE ; 
   const Bool_t image    = kTRUE ; 
   
-  TH1I* hADCCounts = new TH1I("hADCCounts","Dist of ADC counts;ADC counts;Entries",1024,0,1024);
+  TH1I* hADCCounts = new TH1I("hADCCounts",
+                             "Dist of ADC counts;ADC counts;Entries",
+                             1024,0,1024);
   hADCCounts->SetXTitle("ADC counts");
   Add2SDigitsList(hADCCounts, 0, !expert, image);
 }
@@ -107,7 +110,9 @@ void AliFMDQADataMakerSim::InitHits()
   const Bool_t expert   = kTRUE ; 
   const Bool_t image    = kTRUE ; 
   
-  TH1F* hEnergyOfHits = new TH1F("hEnergyOfHits","Energy distribution;Energy [MeV];Counts",100,0,3);
+  TH1F* hEnergyOfHits = new TH1F("hEnergyOfHits",
+                                "Energy distribution;Energy [MeV];Counts",
+                                100,0,3);
   hEnergyOfHits->SetXTitle("Edep");
   hEnergyOfHits->SetYTitle("Counts");
   Add2HitsList(hEnergyOfHits, 0, !expert, image);
@@ -120,7 +125,9 @@ void AliFMDQADataMakerSim::InitDigits()
   const Bool_t expert   = kTRUE ; 
   const Bool_t image    = kTRUE ; 
   
-  TH1I* hADCCounts = new TH1I("hADCCounts","Dist of ADC counts; ADC counts;Entries",1024,0,1024);
+  TH1I* hADCCounts = new TH1I("hADCCounts",
+                             "Dist of ADC counts; ADC counts;Entries",
+                             1024,0,1024);
   hADCCounts->SetXTitle("ADC counts");
   Add2DigitsList(hADCCounts, 0, !expert, image);
 }
@@ -142,11 +149,13 @@ void AliFMDQADataMakerSim::MakeHits()
 void AliFMDQADataMakerSim::MakeHits(TTree * hitTree)
 {
   // make QA data from Hit Tree
-  
-  if (fHitsArray) 
-    fHitsArray->Clear() ; 
-  else
+  // 
+  // Parameters: 
+  //   hitTree    Hits container 
+  //
+  if (!fHitsArray) 
     fHitsArray = new TClonesArray("AliFMDHit", 1000) ; 
+  fHitsArray->Clear() ; 
   
   TBranch * branch = hitTree->GetBranch("FMD") ;
   if (!branch) {
@@ -167,6 +176,9 @@ void AliFMDQADataMakerSim::MakeHits(TTree * hitTree)
 void AliFMDQADataMakerSim::MakeDigits()
 {
   // makes data from Digits
+  // 
+  // Parameters: 
+  //    none
   if(!fDigitsArray) return;
   
   for(Int_t i = 0 ; i < fDigitsArray->GetEntriesFast() ; i++) {
@@ -179,11 +191,14 @@ void AliFMDQADataMakerSim::MakeDigits()
 //_____________________________________________________________________
 void AliFMDQADataMakerSim::MakeDigits(TTree * digitTree)
 {
+  // Make data from digits. 
+  // 
+  // Parameters: 
+  //    digitTree    Tree holding digits. 
   
-  if (fDigitsArray) 
-    fDigitsArray->Clear();
-  else 
+  if (!fDigitsArray) 
     fDigitsArray = new TClonesArray("AliFMDDigit", 1000) ; 
+  fDigitsArray->Clear();
   
   TBranch * branch = digitTree->GetBranch("FMD") ;
   if (!branch)    {
@@ -191,6 +206,9 @@ void AliFMDQADataMakerSim::MakeDigits(TTree * digitTree)
       return;
   } 
   branch->SetAddress(&fDigitsArray) ;
+
+  if (fDigitsArray) fDigitsArray->Clear();
+
   branch->GetEntry(0) ; 
   MakeDigits() ; 
 }
@@ -199,6 +217,9 @@ void AliFMDQADataMakerSim::MakeDigits(TTree * digitTree)
 void AliFMDQADataMakerSim::MakeSDigits()
 {
   // makes data from Digits
+  // 
+  // Parameters: 
+  //   none 
   if(!fSDigitsArray) return;
   
    for(Int_t i = 0 ; i < fSDigitsArray->GetEntriesFast() ; i++) {
@@ -211,11 +232,15 @@ void AliFMDQADataMakerSim::MakeSDigits()
 //_____________________________________________________________________
 void AliFMDQADataMakerSim::MakeSDigits(TTree * sdigitTree)
 {
+  // Make data from digits. 
+  // 
+  // Parameters: 
+  //    digitTree    Tree holding digits. 
   
-  if (fSDigitsArray) 
-    fSDigitsArray->Clear() ;
-  else 
+  if (!fSDigitsArray) 
     fSDigitsArray = new TClonesArray("AliFMDSDigit", 1000) ; 
+  fSDigitsArray->Clear() ;
+
   TBranch * branch = sdigitTree->GetBranch("FMD") ;
   if (!branch)    {
     AliWarning("FMD branch in SDigit Tree not found") ; 
@@ -229,7 +254,10 @@ void AliFMDQADataMakerSim::MakeSDigits(TTree * sdigitTree)
 //_____________________________________________________________________ 
 void AliFMDQADataMakerSim::StartOfDetectorCycle()
 {
-   
+  // Does 
+  // not 
+  // do 
+  // anything 
 }
 //_____________________________________________________________________ 
 //
index 9a50e42..3befa0d 100644 (file)
@@ -120,7 +120,7 @@ AliFMDReconstructor::Init()
   // Current vertex position
   fCurrentVertex = 0;
   // Create array of reconstructed strip multiplicities 
-  fMult = new TClonesArray("AliFMDRecPoint", 51200);
+  // fMult = new TClonesArray("AliFMDRecPoint", 51200);
   // Create ESD output object 
   fESDObj = new AliESDFMD;
   
@@ -162,9 +162,23 @@ AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
 {
   // Convert Raw digits to AliFMDDigit's in a tree 
   AliFMDDebug(1, ("Reading raw data into digits tree"));
+  if (!digitsTree) { 
+    AliError("No digits tree passed");
+    return;
+  }
+  static TClonesArray* array = new TClonesArray("AliFMDDigit");
+  digitsTree->Branch("FMD", &array);
+  
   AliFMDRawReader rawRead(reader, digitsTree);
   // rawRead.SetSampleRate(fFMD->GetSampleRate());
-  rawRead.Exec();
+  // rawRead.Exec();
+  rawRead.ReadAdcs(array);
+
+  Int_t nWrite = digitsTree->Fill();
+  AliFMDDebug(1, ("Got a grand total of %d digits, wrote %d bytes to tree", 
+                  array->GetEntriesFast(), nWrite));
+
+  
   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
   for (size_t i = 1; i <= 3; i++) { 
     fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
@@ -240,6 +254,34 @@ AliFMDReconstructor::UseRecoParam(Bool_t set) const
 }
   
   
+//____________________________________________________________________
+void 
+AliFMDReconstructor::MarkDeadChannels(AliESDFMD* esd) const
+{
+  // Loop over all entries of the ESD and mark 
+  // those that are dead as such 
+  // - otherwise put in the zero signal. 
+  AliFMDParameters* param = AliFMDParameters::Instance();
+
+  for (UShort_t d = 1; d <= 3; d++) { 
+    UShort_t nR = (d == 1 ? 1 : 2);
+
+    for (UShort_t q = 0; q < nR; q++) {
+      Char_t   r  = (q == 0 ? 'I' : 'O');
+      UShort_t nS = (q == 0 ?  20 :  40);
+      UShort_t nT = (q == 0 ? 512 : 256);
+
+      for (UShort_t s = 0; s < nS; s++) { 
+       for (UShort_t t = 0; t < nT; t++) {
+         if (param->IsDead(d, r, s, t)) 
+           esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
+         else if (esd->Multiplicity(d, r, s, t) == AliESDFMD::kInvalidMult) 
+           esd->SetMultiplicity(d, r, s, t, 0);
+       }
+      }
+    }
+  }
+}
 
 //____________________________________________________________________
 void 
@@ -250,7 +292,7 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
   // 
   // Parameters: 
   //   reader  Raw event reader 
-  //   ctree    Not used. 
+  //   ctree    Not used - 'cluster tree' to store rec-points in. 
   AliFMDDebug(1, ("Reconstructing from raw reader"));
   AliFMDRawReader rawReader(reader, 0);
 
@@ -316,18 +358,22 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   //   digitsTree      Pointer to a tree containing digits 
   //   clusterTree     Pointer to output tree 
   // 
+  if (!fMult) fMult = new TClonesArray("AliFMDRecPoint");
+
   AliFMDDebug(1, ("Reconstructing from digits in a tree"));
   GetVertex(fESD);
 
+
+
   static TClonesArray* digits = new TClonesArray("AliFMDDigit");
-  TBranch *digitBranch = digitsTree->GetBranch("FMD");
+  TBranch*      digitBranch   = digitsTree->GetBranch("FMD");
   if (!digitBranch) {
     Error("Exec", "No digit branch for the FMD found");
     return;
   }
-  // TClonesArray* digits = new TClonesArray("AliFMDDigit");
   digitBranch->SetAddress(&digits);
 
+  if (digits)  digits->Clear();
   if (fMult)   fMult->Clear();
   if (fESDObj) fESDObj->Clear();
   
@@ -345,7 +391,7 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 
   Int_t written = clusterTree->Fill();
   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
-  digits->Delete();
+  // digits->Delete();
   // delete digits;
 }
  
@@ -758,14 +804,14 @@ AliFMDReconstructor::PhysicalCoordinates(UShort_t det,
   //    phi     On return, contains the azimuthal angle of the strip
   // 
   AliFMDGeometry* geom = AliFMDGeometry::Instance();
-  Double_t x, y, z, r, theta;
+  Double_t x, y, z, r, theta, deta, dphi;
   geom->Detector2XYZ(det, rng, sec, str, x, y, z);
+
   // Correct for vertex offset. 
   z     -= fCurrentVertex;
-  phi   =  TMath::ATan2(y, x);
-  r     =  TMath::Sqrt(y * y + x * x);
-  theta =  TMath::ATan2(r, z);
-  eta   = -TMath::Log(TMath::Tan(theta / 2));
+  AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, r, deta, dphi, theta);
+  eta = deta;
+  phi = dphi;
 }
 
       
@@ -782,6 +828,9 @@ AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
   // fESDObj->Print();
 
+  // Fix up ESD so that only truely dead channels get the kInvalidMult flag. 
+  MarkDeadChannels(fESDObj);
+
   Double_t oldVz = fCurrentVertex;
   GetVertex(esd);
   if (fVertexType != kNoVertex) { 
index 3b064ba..374a245 100644 (file)
@@ -365,6 +365,14 @@ protected:
                                       Float_t& eta, 
                                       Float_t& phi) const;
   /** 
+   * Mark dead channels as invalid, and those that are marked as invalid 
+   * but are not dead, get the zero signal. 
+   * 
+   * @param esd ESD object to modify. 
+   */
+  void MarkDeadChannels(AliESDFMD* esd) const;
+
+  /** 
    * Set-up reconstructor to use values from reconstruction
    * parameters, if present, for this event.   If the argument @a set
    * is @c false, then restore preset values.