]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - VZERO/AliVZEROTriggerMask.cxx
DA output for test shuttle.
[u/mrichter/AliRoot.git] / VZERO / AliVZEROTriggerMask.cxx
index 0e9452341b13eb7b791245e838cf5de490959385..ee23c57c017378f191084bd63e3e438bf9d3acc5 100644 (file)
  **************************************************************************/
 
 /* $Id$ */
+//-------------------------------------------
+// Class : AliVZEROTriggerMask
+//
+// Fill up the trigger mask word.
+//
+
 #include <Riostream.h>
 #include <TGeoManager.h>
 #include <TGeoMatrix.h>
 #include <TGeoPhysicalNode.h>
+#include <TF1.h>
+#include <TMath.h>
+
 #include <AliGeomManager.h>
-#include "AliCDBManager.h"
-#include "AliCDBStorage.h"
-#include "AliCDBEntry.h"
 #include "AliLog.h"
 #include "AliVZEROTriggerMask.h"
-#include "AliVZEROdigit.h"
+#include "AliVZEROConst.h"
+#include "AliVZEROCalibData.h"
+#include "AliESDVZERO.h"
+#include "AliVZEROReconstructor.h"
 
 //______________________________________________________________________
 ClassImp(AliVZEROTriggerMask)
@@ -33,18 +42,18 @@ ClassImp(AliVZEROTriggerMask)
 
 AliVZEROTriggerMask::AliVZEROTriggerMask()
   :TObject(),
-   fAdcThresHold(0.0),
-   fTimeWindowWidthBBA(50.0),
-   fTimeWindowWidthBGA(20.0),
-   fTimeWindowWidthBBC(50.0),
-   fTimeWindowWidthBGC(20.0),
-   fBBtriggerV0A(0),
-   fBGtriggerV0A(0),
-   fBBtriggerV0C(0),
-   fBGtriggerV0C(0)
-   
+   fV0ADist(0),
+   fV0CDist(0),
+   fRecoParam(NULL)
 {
-   SetAdcThreshold();
+  // Default constructor
+  //
+  Float_t zV0A = TMath::Abs(GetZPosition("VZERO/V0A"));
+  Float_t zV0C = TMath::Abs(GetZPosition("VZERO/V0C"));
+
+  // distance in time units from nominal vertex to V0
+  fV0ADist = zV0A/TMath::Ccgs()*1e9;
+  fV0CDist = zV0C/TMath::Ccgs()*1e9;
 }
 
 //________________________________________________________________________________
@@ -53,8 +62,10 @@ Double_t AliVZEROTriggerMask::GetZPosition(const char* symname){
 //
   Double_t *tr;
   TGeoPNEntry *pne = gGeoManager->GetAlignableEntry(symname);
-  if (!pne) return 0;
-
+  if (!pne) {
+    AliFatalClass(Form("TGeoPNEntry with symbolic name %s does not exist!",symname));
+    return 0;
+  }
 
   TGeoPhysicalNode *pnode = pne->GetPhysicalNode();
   if(pnode){
@@ -63,7 +74,7 @@ Double_t AliVZEROTriggerMask::GetZPosition(const char* symname){
   }else{
           const char* path = pne->GetTitle();
           if(!gGeoManager->cd(path)){
-                  AliErrorClass(Form("Volume path %s not valid!",path));
+                  AliFatalClass(Form("Volume path %s not valid!",path));
                   return 0;
           }
          tr = gGeoManager->GetCurrentMatrix()->GetTranslation();
@@ -75,55 +86,153 @@ Double_t AliVZEROTriggerMask::GetZPosition(const char* symname){
 //________________________________________________________________________________
 
 
-void AliVZEROTriggerMask::FillMasks(TTree* vzeroDigitsTree,
-                               TClonesArray* vzeroDigits)
+void AliVZEROTriggerMask::FillMasks(AliESDVZERO *esdV0,
+                                   AliVZEROCalibData *cal,
+                                   TF1 *slewing)
 {
-  const Double_t LightSpeed = 2.9979245800; // cm/100 ps
-  Float_t Z_V0A = TMath::Abs(GetZPosition("VZERO/V0A"));
-  Float_t Z_V0C = TMath::Abs(GetZPosition("VZERO/V0C"));
+  // Fill up the trigger mask word
+  // using the TDC data (already corrected for
+  // slewing and misalignment between channels)
 
-  // distance in time units from nominal vertex to V0
-  Float_t v0a_dist = Z_V0A/LightSpeed; // 100 of picoseconds
-  Float_t v0c_dist = Z_V0C/LightSpeed; // 100 of picoseconds
-  Float_t bunch_separation = 1000.0; // 100 of picoseconds
-
-  // mask
-  UInt_t one=1;
-  // loop over vzero entries
-  Int_t nEntries = (Int_t)vzeroDigitsTree->GetEntries();
-  for (Int_t e=0; e<nEntries; e++) {
-    vzeroDigitsTree->GetEvent(e);
-
-    Int_t nDigits = vzeroDigits->GetEntriesFast();
-    
-    for (Int_t d=0; d<nDigits; d++) {
-      //      vzeroDigitsTree->GetEvent(d);
-      AliVZEROdigit* digit = (AliVZEROdigit*)vzeroDigits->At(d);
-      
-      Int_t   PMNumber   = digit->PMNumber();
-      Float_t adc        = digit->ADC();
-      Float_t tdc        = digit->Time()*10.0; // in 100 of picoseconds
-
-      if (adc>fAdcThresHold) {
-       if (PMNumber<32) { // in V0C
-         if (tdc>(v0c_dist-fTimeWindowWidthBBC/2.0) &&
-             tdc<(v0c_dist+fTimeWindowWidthBBC/2.0))
-           fBBtriggerV0C+=(one<<PMNumber);
-         if (tdc>(bunch_separation-v0c_dist-fTimeWindowWidthBGC/2.0) &&
-             tdc<(bunch_separation-v0c_dist+fTimeWindowWidthBGC/2.0))
-          fBGtriggerV0C+=(one<<PMNumber); 
+  esdV0->SetBit(AliESDVZERO::kTriggerBitsFilled,kTRUE);
+  esdV0->SetBit(AliESDVZERO::kDecisionFilled,kTRUE);
+
+  const Float_t p1 = 2.50; // photostatistics term in the time resolution
+  const Float_t p2 = 3.00; // slewing related term in the time resolution
+
+  // loop over vzero channels
+  Float_t timeAW = 0,timeCW = 0;
+  Float_t weightA = 0,weightC = 0;
+  Int_t ntimeA = 0, ntimeC = 0;
+  Double_t timesA[32], timesC[32];
+  Double_t wA[32],wC[32];
+  Int_t indA[32], indC[32];
+  for (Int_t i = 0; i < 64; ++i) {
+    Float_t adc = esdV0->GetAdc(i);
+    if (adc > GetRecoParam()->GetAdcThresHold()) {
+      Float_t tdc = esdV0->GetTime(i);
+      if (tdc > (AliVZEROReconstructor::kInvalidTime + 1e-6)) {
+       Float_t nphe = adc*kChargePerADC/(cal->GetGain(i)*TMath::Qe());
+       Float_t timeErr = TMath::Sqrt(kIntTimeRes*kIntTimeRes+
+                                     p1*p1/nphe+
+                                     p2*p2*(slewing->GetParameter(0)*slewing->GetParameter(1))*(slewing->GetParameter(0)*slewing->GetParameter(1))*
+                                     TMath::Power(adc/cal->GetCalibDiscriThr(i,kTRUE),2.*(slewing->GetParameter(1)-1.))/
+                                     (cal->GetCalibDiscriThr(i,kTRUE)*cal->GetCalibDiscriThr(i,kTRUE)));
+
+       if (i < 32) { // in V0C
+         timesC[ntimeC] = tdc;
+         wC[ntimeC] = 1.0/(timeErr*timeErr);
+         indC[ntimeC] = i;
+         ntimeC++;
+         timeCW += tdc/(timeErr*timeErr);
+         weightC += 1.0/(timeErr*timeErr);
        }
-       if (PMNumber>31) { // in V0A
-         Int_t shift = PMNumber-32;
-         if (tdc>(v0a_dist-fTimeWindowWidthBBA/2.0) &&
-             tdc<(v0a_dist+fTimeWindowWidthBBA/2.0)) 
-           fBBtriggerV0A+=(one<<shift);
-         if (tdc>(bunch_separation-v0a_dist-fTimeWindowWidthBGA/2.0) &&
-             tdc<(bunch_separation-v0a_dist+fTimeWindowWidthBGA/2.0))
-           fBGtriggerV0A+=(one<<shift); 
+       else { // in V0A
+         timesA[ntimeA] = tdc;
+         wA[ntimeA] = 1.0/(timeErr*timeErr);
+         indA[ntimeA] = i - 32;
+         ntimeA++;
+         timeAW += tdc/(timeErr*timeErr);
+         weightA += 1.0/(timeErr*timeErr);
        }
       }
-    } // end of loop over digits
-  } // end of loop over events in digits tree
+    }
+  } // end of loop over channels
+
+  if (weightA > 0) timeAW = timeAW/weightA;
+  else timeAW = AliVZEROReconstructor::kInvalidTime;
+
+  if (weightC > 0) timeCW = timeCW/weightC;
+  else timeCW = AliVZEROReconstructor::kInvalidTime;
+
+  esdV0->SetBit(AliESDVZERO::kRobustMeanTime,kTRUE);
+
+  Double_t medianTimeA = AliVZEROReconstructor::kInvalidTime;
+  if (ntimeA > 0) medianTimeA = TMath::Median(ntimeA,timesA,wA);
+  Double_t medianTimeC = AliVZEROReconstructor::kInvalidTime;
+  if (ntimeC > 0) medianTimeC = TMath::Median(ntimeC,timesC,wC);
+
+  Float_t robTimeAW = 0,robTimeCW = 0;
+  Float_t robWeightA = 0,robWeightC = 0;
+  Int_t nrobTimeA = 0, nrobTimeC = 0;
+  Int_t robIndA[32], robIndC[32];
+  for(Int_t i = 0; i < ntimeA; ++i) {
+    AliDebug(1,Form("ChannelsAResiduals %f %f %d",timesA[i]-medianTimeA,1./TMath::Sqrt(wA[i]),ntimeA));
+    if (TMath::Abs(timesA[i]-medianTimeA) < GetRecoParam()->GetMaxResid()/TMath::Sqrt(wA[i])) {
+      robIndA[nrobTimeA] = indA[i];
+      nrobTimeA++;
+      robTimeAW += timesA[i]*wA[i];
+      robWeightA += wA[i];
+    }
+  }
+  for(Int_t i = 0; i < ntimeC; ++i) {
+    AliDebug(1,Form("ChannelsCResiduals %f %f %d",timesC[i]-medianTimeC,1./TMath::Sqrt(wC[i]),ntimeC));
+    if (TMath::Abs(timesC[i]-medianTimeC) < GetRecoParam()->GetMaxResid()/TMath::Sqrt(wC[i])) {
+      robIndC[nrobTimeC] = indC[i];
+      nrobTimeC++;
+      robTimeCW += timesC[i]*wC[i];
+      robWeightC += wC[i];
+    }
+  }
+
+  if (robWeightA > 0) robTimeAW = robTimeAW/robWeightA;
+  else robTimeAW = AliVZEROReconstructor::kInvalidTime;
+
+  if (robWeightC > 0) robTimeCW = robTimeCW/robWeightC;
+  else robTimeCW = AliVZEROReconstructor::kInvalidTime;
+
+  AliDebug(1,Form("V0timesA %f %f %f %f %d",timeAW,(weightA > 0) ? (1./TMath::Sqrt(weightA)) : 0,
+                 medianTimeA,robTimeAW,ntimeA));
+  AliDebug(1,Form("V0timesC %f %f %f %f %d",timeCW,(weightC > 0) ? (1./TMath::Sqrt(weightC)) : 0,
+                 medianTimeC,robTimeCW,ntimeC));
+
+  esdV0->SetV0ATime(robTimeAW);
+  esdV0->SetV0CTime(robTimeCW);
+  esdV0->SetV0ATimeError((robWeightA > 0) ? (1./TMath::Sqrt(robWeightA)) : 0);
+  esdV0->SetV0CTimeError((robWeightC > 0) ? (1./TMath::Sqrt(robWeightC)) : 0);
+
+  esdV0->SetV0ADecision(AliESDVZERO::kV0Empty);
+  esdV0->SetV0CDecision(AliESDVZERO::kV0Empty);
+
+  if (robTimeAW > (fV0ADist + GetRecoParam()->GetTimeWindowBBALow()) &&
+      robTimeAW < (fV0ADist + GetRecoParam()->GetTimeWindowBBAUp())) 
+    esdV0->SetV0ADecision(AliESDVZERO::kV0BB);
+  else if (robTimeAW > (-fV0ADist + GetRecoParam()->GetTimeWindowBGALow()) &&
+          robTimeAW < (-fV0ADist + GetRecoParam()->GetTimeWindowBGAUp()))
+    esdV0->SetV0ADecision(AliESDVZERO::kV0BG);
+  else if (robTimeAW > (AliVZEROReconstructor::kInvalidTime + 1e-6))
+    esdV0->SetV0ADecision(AliESDVZERO::kV0Fake);
+
+  if (robTimeCW > (fV0CDist + GetRecoParam()->GetTimeWindowBBCLow()) &&
+      robTimeCW < (fV0CDist + GetRecoParam()->GetTimeWindowBBCUp())) 
+    esdV0->SetV0CDecision(AliESDVZERO::kV0BB);
+  else if (robTimeCW > (-fV0CDist + GetRecoParam()->GetTimeWindowBGCLow()) &&
+          robTimeCW < (-fV0CDist + GetRecoParam()->GetTimeWindowBGCUp()))
+    esdV0->SetV0CDecision(AliESDVZERO::kV0BG);
+  else if (robTimeCW > (AliVZEROReconstructor::kInvalidTime + 1e-6))
+    esdV0->SetV0CDecision(AliESDVZERO::kV0Fake);
+
+  UInt_t aBBtriggerV0A = 0; // bit mask for Beam-Beam trigger in V0A
+  UInt_t aBGtriggerV0A = 0; // bit mask for Beam-Gas trigger in V0A
+  UInt_t aBBtriggerV0C = 0; // bit mask for Beam-Beam trigger in V0C
+  UInt_t aBGtriggerV0C = 0; // bit mask for Beam-Gas trigger in V0C
+
+  for(Int_t i = 0; i < nrobTimeA; ++i) {
+    if (esdV0->GetV0ADecision() == AliESDVZERO::kV0BB)
+      aBBtriggerV0A |= (1 << (robIndA[i]));
+    else if (esdV0->GetV0ADecision() == AliESDVZERO::kV0BG)
+      aBGtriggerV0A |= (1 << (robIndA[i]));
+  }
+
+  for(Int_t i = 0; i < nrobTimeC; ++i) {
+    if (esdV0->GetV0CDecision() == AliESDVZERO::kV0BB)
+      aBBtriggerV0C |= (1 << (robIndC[i]));
+    else if (esdV0->GetV0CDecision() == AliESDVZERO::kV0BG)
+      aBGtriggerV0C |= (1 << (robIndC[i]));
+  }
+
+  esdV0->SetBBtriggerV0A(aBBtriggerV0A);
+  esdV0->SetBGtriggerV0A(aBGtriggerV0A);
+  esdV0->SetBBtriggerV0C(aBBtriggerV0C);
+  esdV0->SetBGtriggerV0C(aBGtriggerV0C);
 }