]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDmcmSim.cxx
Changes for PbPb RAA analysis (Zaida)
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmSim.cxx
index 21dcd667135f455ae5dddf19297028006a475550..247dee10703391062841b998820bb6d0aa6eaead 100644 (file)
@@ -35,6 +35,7 @@
 #include "TRandom.h"
 #include "TClonesArray.h"
 #include "TMath.h"
+#include <TTree.h>
 
 #include "AliLog.h"
 #include "AliRunLoader.h"
@@ -88,6 +89,14 @@ AliTRDmcmSim::AliTRDmcmSim() :
   // AliTRDmcmSim default constructor
   // By default, nothing is initialized.
   // It is necessary to issue Init before use.
+
+  for (Int_t iDict = 0; iDict < 3; iDict++)
+    fDict[iDict] = 0x0;
+
+  fFitPtr[0] = 0;
+  fFitPtr[1] = 0;
+  fFitPtr[2] = 0;
+  fFitPtr[3] = 0;
 }
 
 AliTRDmcmSim::~AliTRDmcmSim() 
@@ -540,17 +549,18 @@ void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *
           continue;
 
         fDict[iDict] = newDict;
-        
-        if (fDict[iDict]->GetDim() == 0) {
-          AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
-          continue;
-        }
         fDict[iDict]->Expand(); 
       }
       else {
         fDict[iDict] = newDict;
         if (fDict[iDict])
           fDict[iDict]->Expand();
+        }
+      
+      // If there is no data, set dictionary to zero to avoid crashes  
+      if (fDict[iDict]->GetDim() == 0)  {
+        AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
+        fDict[iDict] = 0x0;
       }
     }
   }
@@ -594,11 +604,6 @@ void AliTRDmcmSim::SetDataByPad(AliTRDarrayADC* const adcArray, AliTRDdigitsMana
           continue;
 
         fDict[iDict] = newDict;
-        
-        if (fDict[iDict]->GetDim() == 0) {
-          AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
-          continue;
-        }
         fDict[iDict]->Expand(); 
       }
       else {
@@ -606,6 +611,12 @@ void AliTRDmcmSim::SetDataByPad(AliTRDarrayADC* const adcArray, AliTRDdigitsMana
         if (fDict[iDict])
           fDict[iDict]->Expand();
       }
+      
+      // If there is no data, set dictionary to zero to avoid crashes  
+      if (fDict[iDict]->GetDim() == 0)  {
+        AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
+        fDict[iDict] = 0x0;
+      }
     }
   }
 
@@ -668,7 +679,7 @@ Bool_t AliTRDmcmSim::GetHit(Int_t index, Int_t &channel, Int_t &timebin, Int_t &
                         (channel << 8) - ypos) 
     * (0.635 + 0.03 * (fDetector % 6))
     / 256.0;
-  label   = fHits[index].fLabel;
+  label   = fHits[index].fLabel[0];
 
   return kTRUE;
 }
@@ -1147,7 +1158,7 @@ void AliTRDmcmSim::ZSMapping()
   }
 }
 
-void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label
+void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label[])
 {
   // Add the given hit to the fit register which is lateron used for 
   // the tracklet calculation. 
@@ -1178,7 +1189,9 @@ void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Sh
   fHits[fNHits].fQtot = qtot;
   fHits[fNHits].fYpos = ypos;
   fHits[fNHits].fTimebin = timebin;
-  fHits[fNHits].fLabel = label;
+  fHits[fNHits].fLabel[0] = label[0];
+  fHits[fNHits].fLabel[1] = label[1];
+  fHits[fNHits].fLabel[2] = label[2];
   fNHits++;
 }
 
@@ -1348,16 +1361,16 @@ void AliTRDmcmSim::CalcFitreg()
         ypos = 128*(adcLeft - adcRight) / adcCentral;
         if (ypos < 0) ypos = -ypos;
         // make the correction using the position LUT
-        ypos = ypos + fTrapConfig->GetTrapReg((AliTRDtrapConfig::TrapReg_t) (AliTRDtrapConfig::kTPL00 + (ypos & 0x7F))); 
+        ypos = ypos + fTrapConfig->GetTrapReg((AliTRDtrapConfig::TrapReg_t) (AliTRDtrapConfig::kTPL00 + (ypos & 0x7F)),
+                                             fDetector, fRobPos, fMcmPos);
         if (adcLeft > adcRight) ypos = -ypos;
 
-        // label calculation
-        Int_t mcLabel = -1;
+        // label calculation (up to 3)
+        Int_t mcLabel[] = {-1, -1, -1};
         if (fDigitsManager) {
-          Int_t label[9] = { 0 }; // up to 9 different labels possible
-          Int_t count[9] = { 0 };
-          Int_t maxIdx = -1;
-          Int_t maxCount = 0;
+          const Int_t maxLabels = 9;
+          Int_t label[maxLabels] = { 0 }; // up to 9 different labels possible
+          Int_t count[maxLabels] = { 0 };
           Int_t nLabels = 0;
           Int_t padcol[3]; 
           padcol[0] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch);
@@ -1370,26 +1383,29 @@ void AliTRDmcmSim::CalcFitreg()
             for (Int_t iPad = 0; iPad < 3; iPad++) {
               if (padcol[iPad] < 0) 
                 continue;
-              Int_t currLabel = fDict[iDict]->GetData(padrow, padcol[iPad], timebin); //fDigitsManager->GetTrack(iDict, padrow, padcol, timebin, fDetector);
+              Int_t currLabel = fDict[iDict]->GetData(padrow, padcol[iPad], timebin);
              AliDebug(10, Form("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin));
               for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
                 if (currLabel == label[iLabel]) {
                   count[iLabel]++;
-                  if (count[iLabel] > maxCount) {
-                    maxCount = count[iLabel];
-                    maxIdx = iLabel;
-                  }
                   currLabel = -1;
                   break;
                 }
               } 
               if (currLabel >= 0) {
-                label[nLabels++] = currLabel;
+                label[nLabels] = currLabel;
+               count[nLabels] = 1;
+               nLabels++;
               }
             }
           }
-          if (maxIdx >= 0)
-            mcLabel = label[maxIdx];
+         Int_t index[2*maxLabels];
+         TMath::Sort(maxLabels, count, index);
+         for (Int_t i = 0; i < 3; i++) {
+           if (count[index[i]] <= 0)
+             break;
+           mcLabel[i] = label[index[i]];
+         }
         }
 
         // add the hit to the fitregister
@@ -1492,12 +1508,13 @@ void AliTRDmcmSim::FitTracklet()
   // which have been filled in the fit registers. 
 
   // parameters in fitred.asm (fit program)
-  Int_t decPlaces = 5;
   Int_t rndAdd = 0;
-  if (decPlaces >  1) 
+  Int_t decPlaces = 5; // must be larger than 1 or change the following code
+  // if (decPlaces >  1)
     rndAdd = (1 << (decPlaces-1)) + 1;
-  else if (decPlaces == 1)
-    rndAdd = 1;
+  // else if (decPlaces == 1)
+  //   rndAdd = 1;
+
   Int_t ndriftDp = 5;  // decimal places for drift time
   Long64_t shift = ((Long64_t) 1 << 32);
 
@@ -1547,7 +1564,7 @@ void AliTRDmcmSim::FitTracklet()
       nHits   = fit0->fNhits + fit1->fNhits; // number of hits
       sumX    = fit0->fSumX  + fit1->fSumX;
       sumX2   = fit0->fSumX2 + fit1->fSumX2;
-      denom   = nHits*sumX2 - sumX*sumX;
+      denom   = ((Long64_t) nHits)*((Long64_t) sumX2) - ((Long64_t) sumX)*((Long64_t) sumX);
 
       mult    = mult / denom; // exactly like in the TRAP program
       q0      = fit0->fQ0    + fit1->fQ0;
@@ -1614,7 +1631,7 @@ void AliTRDmcmSim::FitTracklet()
       else
       {
         if (slope > 63 || slope < -64) { // wrapping in TRAP!
-          AliError(Form("Overflow in slope: %i, tracklet discarded!", slope));
+          AliDebug(1,Form("Overflow in slope: %i, tracklet discarded!", slope));
           fMCMT[cpu] = 0x10001000;
           continue;
         }
@@ -1635,14 +1652,13 @@ void AliTRDmcmSim::FitTracklet()
         fMCMT[cpu] = (pid << 24) | (padrow << 20) | (slope << 13) | offset;
 
         // calculate MC label
-        Int_t mcLabel = -1;
+        Int_t mcLabel[] = { -1, -1, -1};
        Int_t nHits0 = 0;
        Int_t nHits1 = 0;
         if (fDigitsManager) {
-          Int_t label[30] = {0}; // up to 30 different labels possible
-          Int_t count[30] = {0};
-          Int_t maxIdx = -1;
-          Int_t maxCount = 0;
+         const Int_t maxLabels = 30;
+          Int_t label[maxLabels] = {0}; // up to 30 different labels possible
+          Int_t count[maxLabels] = {0};
           Int_t nLabels = 0;
           for (Int_t iHit = 0; iHit < fNHits; iHit++) {
             if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) ||
@@ -1657,24 +1673,29 @@ void AliTRDmcmSim::FitTracklet()
                fHits[iHit].fTimebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))
              nHits1++;
 
-            Int_t currLabel = fHits[iHit].fLabel;
-            for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
-              if (currLabel == label[iLabel]) {
-                count[iLabel]++;
-                if (count[iLabel] > maxCount) {
-                  maxCount = count[iLabel];
-                  maxIdx = iLabel;
-                }
-                currLabel = -1;
-                break;
-              }
-            }
-            if (currLabel >= 0) {
-              label[nLabels++] = currLabel;
-            }
-          }
-          if (maxIdx >= 0)
-            mcLabel = label[maxIdx];
+           for (Int_t i = 0; i < 3; i++) {
+             Int_t currLabel = fHits[iHit].fLabel[i];
+             for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
+               if (currLabel == label[iLabel]) {
+                 count[iLabel]++;
+                 currLabel = -1;
+                 break;
+               }
+             }
+             if (currLabel >= 0 && nLabels < maxLabels) {
+               label[nLabels] = currLabel;
+               count[nLabels]++;
+               nLabels++;
+             }
+           }
+         }
+         Int_t index[2*maxLabels];
+         TMath::Sort(maxLabels, count, index);
+         for (Int_t i = 0; i < 3; i++) {
+           if (count[index[i]] <= 0)
+             break;
+           mcLabel[i] = label[index[i]];
+         }
         }
         new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
@@ -1683,8 +1704,8 @@ void AliTRDmcmSim::FitTracklet()
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits(fit0->fNhits + fit1->fNhits);
        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits0(nHits0);
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits1(nHits1);
-        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0);
-        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1);
+        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0 >> fgkAddDigits);
+        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1 >> fgkAddDigits);
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetSlope(fitSlope);
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetOffset(fitOffset);
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetError(TMath::Sqrt(TMath::Abs(fitError)/nHits));
@@ -1860,12 +1881,12 @@ Int_t AliTRDmcmSim::GetPID(Int_t q0, Int_t q1)
 
    UInt_t nBinsQ0 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTnbins);  // number of bins in q0 / 4 !!
    UInt_t pidTotalSize = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTLength);
-   if(nBinsQ0==0 || pidTotalSize==0)  // make sure we don't run into trouble if one of the values is not configured
-      return 0;
+   if(nBinsQ0==0 || pidTotalSize==0)  // make sure we don't run into trouble if the value for Q0 is not configured
+     return 0;                        // Q1 not configured is ok for 1D LUT
 
    ULong_t corrQ0 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTcor0, fDetector, fRobPos, fMcmPos);
    ULong_t corrQ1 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTcor1, fDetector, fRobPos, fMcmPos);
-   if(corrQ0==0 || corrQ1==0)  // make sure we don't run into trouble if one of the values is not configured
+   if(corrQ0==0)  // make sure we don't run into trouble if one of the values is not configured
       return 0;
 
    addrQ0 = corrQ0;
@@ -2166,7 +2187,7 @@ ostream& operator<<(ostream& os, const AliTRDmcmSim& mcm)
     Int_t bufLength   = mcm.ProduceRawStream(&buf[0], bufSize);
     
     for (Int_t i = 0; i < bufLength; i++) 
-      std::cout << "0x" << std::hex << buf[i] << std::endl;
+      std::cout << "0x" << std::hex << buf[i] << std::dec << std::endl;
     
     delete [] buf;
   }