]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPC.cxx
Script supersceeded by AliForwarddNdetaTask.C and
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
index a2948ae7139e4c31c9bc5afb01846d2b4121f7fc..04ce0a97156836d75d5b3da22f7cdf70ed70ff4c 100644 (file)
@@ -52,6 +52,7 @@
 #include <TTree.h>
 #include <TVector.h>
 #include <TVirtualMC.h>
+#include <TParameter.h>
 
 #include "AliDigits.h"
 #include "AliMagF.h"
@@ -78,7 +79,7 @@
 #include "AliTPCCalROC.h"
 #include "AliTPCExB.h"
 #include "AliRawReader.h"
-#include "AliTPCRawStream.h"
+#include "AliTPCRawStreamV3.h"
 #include "TTreeStream.h"
 
 ClassImp(AliTPC) 
@@ -99,7 +100,8 @@ ClassImp(AliTPC)
                   fCurrentNoise(0),
                   fActiveSectors(0),
                    fGainFactor(1.),
-    fDebugStreamer(0)
+                   fDebugStreamer(0),
+                   fLHCclockPhaseSw(0)
 
 {
   //
@@ -132,8 +134,9 @@ AliTPC::AliTPC(const char *name, const char *title)
                   fNoiseTable(0),
                   fCurrentNoise(0),
                    fActiveSectors(0),
-    fGainFactor(1.),
-     fDebugStreamer(0)
+                   fGainFactor(1.),
+                   fDebugStreamer(0),
+                   fLHCclockPhaseSw(0)
                   
 {
   //
@@ -184,7 +187,12 @@ AliTPC::AliTPC(const char *name, const char *title)
   }
 
 }
-
+void AliTPC::CreateDebugStremer(){
+  //
+  // Create Debug streamer to check simulation
+  // 
+  fDebugStreamer = new TTreeSRedirector("TPCSimdebug.root");
+}
 //_____________________________________________________________________________
 AliTPC::~AliTPC()
 {
@@ -259,7 +267,7 @@ void AliTPC::CreateMaterials()
   wmat[0]=0.2729;
   wmat[1]=0.7271;
 
-  density=0.0018872;
+  density=0.001754609;
 
 
   AliMixture(10,"CO2",amat,zmat,density,2,wmat);
@@ -304,11 +312,11 @@ void AliTPC::CreateMaterials()
   wmat[2]=0.128469474;
   wmat[3]=0.058395789;
  
-  density=0.0009733;
+  density=0.000904929;
 
   AliMixture(12,"Ne-CO2-N-1",amat,zmat,density,4,wmat);
   AliMixture(13,"Ne-CO2-N-2",amat,zmat,density,4,wmat);
-  AliMixture(30,"Ne-CO2-N-3",amat,zmat,density,4,wmat);
+  AliMixture(35,"Ne-CO2-N-3",amat,zmat,density,4,wmat);
   //----------------------------------------------------------------------
   //               solid materials
   //----------------------------------------------------------------------
@@ -674,7 +682,7 @@ void AliTPC::CreateMaterials()
   AliMedium(1, "Ne-CO2-N-1", 12, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
   AliMedium(2, "Ne-CO2-N-2", 13, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001); 
-  AliMedium(20, "Ne-CO2-N-3", 30, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
+  AliMedium(20, "Ne-CO2-N-3", 35, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
   //-----------------------------------------------------------  
   // tracking media for solids
   //-----------------------------------------------------------
@@ -775,7 +783,7 @@ void    AliTPC::SetActiveSectors(Int_t flag)
   else branch = fLoader->TreeH()->GetBranch("TPC");
   Stat_t ntracks = fLoader->TreeH()->GetEntries();
   // loop over all hits
-  AliDebug(1,Form("Got %d tracks",ntracks));
+  AliDebug(1,Form("Got %d tracks", (Int_t) ntracks));
   
   for(Int_t track=0;track<ntracks;track++) {
     ResetHits();
@@ -828,6 +836,7 @@ void AliTPC::Digits2Raw()
     return;
   }
 
+  //
   AliSimDigits digarr;
   AliSimDigits* digrow = &digarr;
   digits->GetBranch("Segment")->SetAddress(&digrow);
@@ -938,8 +947,17 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
   const Int_t kNOS = fTPCParam->GetNOuterSector();
   const Int_t kNS = kNIS + kNOS;
 
-  Short_t** allBins = NULL; //array which contains the data for one sector
-  
+  // Setup storage
+  AliTPCROC * roc = AliTPCROC::Instance();
+  Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
+  Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
+  Short_t** allBins = new Short_t*[nRowsMax];
+  for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
+    Int_t maxBin = kmaxTime*nPadsMax;
+    allBins[iRow] = new Short_t[maxBin];
+    memset(allBins[iRow],0,sizeof(Short_t)*maxBin);
+  }
+
   for(Int_t iSector = 0; iSector < kNS; iSector++) {
     
     Int_t nRows = fTPCParam->GetNRow(iSector);
@@ -953,55 +971,65 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
       indexDDL = (iSector-kNIS) * 4 + kNIS * 2;
     }
 
-    // Loas the raw data for corresponding DDLs
+    // Load the raw data for corresponding DDLs
     rawReader->Reset();
-    AliTPCRawStream input(rawReader);
+
+    AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
+    AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 
-    // Alocate and init the array with the sector data
-    allBins = new Short_t*[nRows];
-    for (Int_t iRow = 0; iRow < nRows; iRow++) {
-      Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
-      Int_t maxBin = kmaxTime*maxPad;
-      allBins[iRow] = new Short_t[maxBin];
+    // Clean storage
+    for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
+      Int_t maxBin = kmaxTime*nPadsMax;
       memset(allBins[iRow],0,sizeof(Short_t)*maxBin);
     }
 
     // Begin loop over altro data
-    while (input.Next()) {
+    while (input.NextDDL()) {
 
       if (input.GetSector() != iSector)
        AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSector,input.GetSector()));
 
-      Int_t iRow = input.GetRow();
-      if (iRow < 0 || iRow >= nRows)
-       AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
-                     iRow, 0, nRows -1));
-      Int_t iPad = input.GetPad();
-
-      Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
-
-      if (iPad < 0 || iPad >= maxPad)
-       AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
-                     iPad, 0, maxPad -1));
-
-      Int_t iTimeBin = input.GetTime();
-      if ( iTimeBin < 0 || iTimeBin >= kmaxTime)
-       AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
-                     iTimeBin, 0, kmaxTime -1));
-      
-      Int_t maxBin = kmaxTime*maxPad;
-
-      if (((iPad*kmaxTime+iTimeBin) >= maxBin) ||
-         ((iPad*kmaxTime+iTimeBin) < 0))
-       AliFatal(Form("Index outside the allowed range"
-                     " Sector=%d Row=%d Pad=%d Timebin=%d"
-                     " (Max.index=%d)",iSector,iRow,iPad,iTimeBin,maxBin));
-
-      allBins[iRow][iPad*kmaxTime+iTimeBin] = input.GetSignal();
+      //loop over pads
+      while ( input.NextChannel() ) {
+
+        Int_t iRow = input.GetRow();
+        if (iRow < 0 || iRow >= nRows)
+          AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
+                        iRow, 0, nRows -1));
+        Int_t iPad = input.GetPad();
+
+        Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
+
+        if (iPad < 0 || iPad >= maxPad)
+          AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
+                        iPad, 0, maxPad -1));
+
+        //loop over bunches
+        while ( input.NextBunch() ){
+          Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
+          Int_t  bunchlength  = (Int_t)input.GetBunchLength();
+          const UShort_t *sig = input.GetSignals();
+          for (Int_t iTime = 0; iTime<bunchlength; iTime++){
+            Int_t iTimeBin=startTbin-iTime;
+            if ( iTimeBin < 0 || iTimeBin >= kmaxTime) {
+              continue;
+              //AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
+              //               iTimeBin, 0, kmaxTime -1));
+            }
+
+            Int_t maxBin = kmaxTime*maxPad;
+            if (((iPad*kmaxTime+iTimeBin) >= maxBin) ||
+                ((iPad*kmaxTime+iTimeBin) < 0))
+              AliFatal(Form("Index outside the allowed range"
+                            " Sector=%d Row=%d Pad=%d Timebin=%d"
+                            " (Max.index=%d)",iSector,iRow,iPad,iTimeBin,maxBin));
+            allBins[iRow][iPad*kmaxTime+iTimeBin] = sig[iTime];
+          }
+        }
+      } // End loop over altro data
+    }
 
-    } // End loop over altro data
-    
     // Now fill the digits array
     if (fDigitsArray->GetTree()==0) {
       AliFatal("Tree not set in fDigitsArray");
@@ -1009,7 +1037,6 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
 
     for (Int_t iRow = 0; iRow < nRows; iRow++) {
       AliDigits * dig = fDigitsArray->CreateRow(iSector,iRow);
-
       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
       for(Int_t iPad = 0; iPad < maxPad; iPad++) {
        for(Int_t iTimeBin = 0; iTimeBin < kmaxTime; iTimeBin++) {
@@ -1029,19 +1056,35 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
       fDigitsArray->ClearRow(iSector,iRow);  
 
     } // end of the sector digitization
-
-    for (Int_t iRow = 0; iRow < nRows; iRow++)
-      delete [] allBins[iRow];
-
-    delete [] allBins;
-
   }
+  // get LHC clock phase from the digits tree
 
-  fLoader->WriteSDigits("OVERWRITE");
+  TParameter<float> *ph; 
+  Float_t phase;
+  TTree *digtree = fLoader->TreeD();
+  //
+  if(digtree){ // if TreeD exists
+    ph = (TParameter<float>*)digtree->GetUserInfo()->FindObject("lhcphase0");
+    phase = ph->GetVal();
+  }
+  else{ //TreeD does not exist
+    phase = 0.; 
+  }
+    //
+    // store lhc clock phase in S-digits tree
+    //
+    fLoader->TreeS()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",phase));
+   //
+   fLoader->WriteSDigits("OVERWRITE");
 
   if(GetDigitsArray()) delete GetDigitsArray();
   SetDigitsArray(0x0);
 
+  // cleanup storage
+  for (Int_t iRow = 0; iRow < nRowsMax; iRow++)
+    delete [] allBins[iRow];
+  delete [] allBins;
+
   return kTRUE;
 }
 
@@ -1309,14 +1352,32 @@ void AliTPC::Hits2Digits(Int_t eventnumber)
   SetDigitsArray(arr);
 
   fDigitsSwitch=0; // standard digits
-
+  // here LHC clock phase
+  Float_t lhcph = 0.;
+  switch (fLHCclockPhaseSw){
+  case 0: 
+    // no phase
+    lhcph=0.;
+    break;
+  case 1:
+    // random phase
+    lhcph = (Int_t)(gRandom->Rndm()/0.25);    
+    break;
+  case 2:
+    lhcph=0.;
+    // not implemented yet
+    break;
+  }
+  // adding phase to the TreeD user info 
+  fLoader->TreeD()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",lhcph));
+  //
   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) 
     if (IsSectorActive(isec)) {
-      AliDebug(1,Form("Hits2Digits","Sector %d is active.",isec));
+      AliDebug(1,Form("Hits2DigitsSector %d is active.",isec));
       Hits2DigitsSector(isec);
     }
     else {
-      AliDebug(1,Form("Hits2Digits","Sector %d is NOT active.",isec));
+      AliDebug(1,Form("Hits2DigitsSector %d is NOT active.",isec));
     }
   
   fLoader->WriteDigits("OVERWRITE"); 
@@ -1376,6 +1437,25 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber)
   fDigitsSwitch=1; // summable digits
   
     // set zero suppression to "0"
+  // here LHC clock phase
+  Float_t lhcph = 0.;
+  switch (fLHCclockPhaseSw){
+  case 0: 
+    // no phase
+    lhcph=0.;
+    break;
+  case 1:
+    // random phase
+    lhcph = (Int_t)(gRandom->Rndm()/0.25);    
+    break;
+  case 2:
+    lhcph=0.;
+    // not implemented yet
+    break;
+  }
+  // adding phase to the TreeS user info 
+  
+  fLoader->TreeS()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",lhcph));
 
   fTPCParam->SetZeroSup(0);
 
@@ -1401,7 +1481,6 @@ void AliTPC::Hits2SDigits()
   //-----------------------------------------------------------
   //   summable digits - 16 bit "ADC", no noise, no saturation
   //-----------------------------------------------------------
-  if (0) fDebugStreamer = new TTreeSRedirector("TPCSimdebug.root");
 
   if (!fTPCParam->IsGeoRead()){
     //
@@ -1420,7 +1499,7 @@ void AliTPC::Hits2SDigits()
     SetActiveSectors();
     Hits2SDigits2(iEvent);
   }
-
+  
   fLoader->UnloadHits();
   fLoader->UnloadSDigits();
   if (fDebugStreamer) {
@@ -1681,6 +1760,15 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
   TMatrixF &signal = *m1;
   TMatrixF &total = *m2;
   //
+  // Get LHC clock phase
+  //
+  TParameter<float> *ph;
+  if(fDigitsSwitch){// s-digits
+    ph = (TParameter<float>*)fLoader->TreeS()->GetUserInfo()->FindObject("lhcphase0");  
+  }
+  else{ // normal digits
+    ph = (TParameter<float>*)fLoader->TreeD()->GetUserInfo()->FindObject("lhcphase0");
+  } 
   //  Loop over all electrons
   //
   for(Int_t nel=0; nel<nElectrons; nel++){
@@ -1688,7 +1776,8 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
     Float_t aval =  v(idx+4);
     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac(); 
     Float_t xyz[4]={v(idx+1),v(idx+2),v(idx+3),v(idx+5)};
-    Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
+    Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,
+                                                           fCurrentIndex[3],ph->GetVal());
 
     Int_t *index = fTPCParam->GetResBin(0);  
     Float_t *weight = & (fTPCParam->GetResWeight(0));