]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPC.cxx
Corrected raw->sdigits method. Using default altro raw-data reader class. Avoid alloc...
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
index 3b9b00677575c0c03674162ead8744fee1a2fde3..af8bc32128b5d6871e480467dcd9767601026f5c 100644 (file)
 #include <Riostream.h>
 #include <stdlib.h>
 
+#include <TF2.h>
 #include <TFile.h>  
+#include <TGeoGlobalMagField.h>
 #include <TInterpreter.h>
 #include <TMath.h>
 #include <TMatrixF.h>
-#include <TVector.h>
 #include <TObjectTable.h>
 #include <TParticle.h>
 #include <TROOT.h>
 #include <TRandom.h>
+#include <TStopwatch.h>
+#include <TString.h>
 #include <TSystem.h>     
 #include <TTree.h>
+#include <TVector.h>
 #include <TVirtualMC.h>
-#include <TString.h>
-#include <TF2.h>
-#include <TStopwatch.h>
+#include <TParameter.h>
 
 #include "AliDigits.h"
 #include "AliMagF.h"
@@ -64,7 +66,6 @@
 #include "AliTPCPRF2D.h"
 #include "AliTPCParamSR.h"
 #include "AliTPCRF1D.h"
-//#include "AliTPCTrackHits.h"
 #include "AliTPCTrackHitsV2.h"
 #include "AliTrackReference.h"
 #include "AliMC.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()
 {
@@ -230,8 +238,8 @@ void AliTPC::CreateMaterials()
   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
   //-----------------------------------------------------------------
 
-   Int_t iSXFLD=gAlice->Field()->Integ();
-  Float_t sXMGMX=gAlice->Field()->Max();
+   Int_t iSXFLD=((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();
+  Float_t sXMGMX=((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();
 
   Float_t amat[5]; // atomic numbers
   Float_t zmat[5]; // z
@@ -259,7 +267,7 @@ void AliTPC::CreateMaterials()
   wmat[0]=0.2729;
   wmat[1]=0.7271;
 
-  density=0.001977;
+  density=0.001754609;
 
 
   AliMixture(10,"CO2",amat,zmat,density,2,wmat);
@@ -283,9 +291,11 @@ void AliTPC::CreateMaterials()
   // drift gases 
   //----------------------------------------------------------------
 
-
+  //
   // Drift gases 1 - nonsensitive, 2 - sensitive
-  // Ne-CO2-N (85-10-5)
+  // Ne-CO2-N2 (90-10-5) (volume) values at 20deg and 1 atm.
+  // rho(Ne) = 0.839 g/cm^3, rho(CO2) = 1.842 g/cm^3, rho(N2) = 1.165 g/cm^3
+  // for the calculation - everything is normalized to 1
 
   amat[0]= 20.18;
   amat[1]=12.011;
@@ -297,12 +307,12 @@ void AliTPC::CreateMaterials()
   zmat[2]=8.;
   zmat[3]=7.;
 
-  wmat[0]=0.7707;
-  wmat[1]=0.0539;
-  wmat[2]=0.1438;
-  wmat[3]=0.0316;
+  wmat[0]=0.756992632;
+  wmat[1]=0.056235789;
+  wmat[2]=0.128469474;
+  wmat[3]=0.058395789;
  
-  density=0.0010252;
+  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);
@@ -443,7 +453,7 @@ void AliTPC::CreateMaterials()
   wmat[0]=0.194;
   wmat[1]=0.023;
   wmat[2]=0.443;
-  wmat[3]=0.340;
+  wmat[3]=0.34;
 
   density=1.82;
 
@@ -461,12 +471,12 @@ void AliTPC::CreateMaterials()
   zmat[2]=8.;
   zmat[3]=14.;
 
-  wmat[0]=0.225;
+  wmat[0]=0.257;
   wmat[1]=0.03;
-  wmat[2]=0.443;
+  wmat[2]=0.412;
   wmat[3]=0.3;
 
-  density=1.163;
+  density=1.725;
 
   AliMixture(21, "Prepreg3",amat,zmat,density,4,wmat);
 
@@ -518,6 +528,24 @@ void AliTPC::CreateMaterials()
 
   AliMaterial(25,"Cu",amat[0],zmat[0],density,999.,999.);
 
+  // brass
+
+  amat[0] = 63.546;
+  zmat[0] = 29.;
+  //
+  amat[1]= 65.409;
+  zmat[1]= 30.;
+  //
+  wmat[0]= 0.6;
+  wmat[1]= 0.4;
+
+  //
+  density = 8.23;
+  
+  //
+  AliMixture(33,"Brass",amat,zmat,density,2,wmat);
+  
   // Epoxy - C14 H20 O3
  
   amat[0]=12.011;
@@ -535,6 +563,28 @@ void AliTPC::CreateMaterials()
   density=1.25;
 
   AliMixture(26,"Epoxy",amat,zmat,density,-3,wmat);
+  //
+  // epoxy film - 90% epoxy, 10% glass fiber 
+  //
+  amat[0]=12.01;
+  amat[1]=1.;
+  amat[2]=15.994;
+  amat[3]=28.086;
+
+  zmat[0]=6.;
+  zmat[1]=1.;
+  zmat[2]=8.;
+  zmat[3]=14.;
+
+  wmat[0]=0.596;
+  wmat[1]=0.071;
+  wmat[2]=0.257;
+  wmat[3]=0.076;
+
+
+  density=1.345;
+
+  AliMixture(34, "Epoxy-film",amat,zmat,density,4,wmat);
 
   // Plexiglas  C5H8O2
 
@@ -622,6 +672,7 @@ void AliTPC::CreateMaterials()
   density=1.;
 
   AliMixture(32,"Water",amat,zmat,density,-2,wmat);  
+
  
   //----------------------------------------------------------
   // tracking media for gases
@@ -655,7 +706,9 @@ void AliTPC::CreateMaterials()
   AliMedium(18,"Steel",29,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001); 
   AliMedium(19,"Peek",30,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
   AliMedium(21,"Alumina",31,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);    
-  AliMedium(22,"Water",32,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);  
+  AliMedium(22,"Water",32,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
+  AliMedium(23,"Brass",33,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
+  AliMedium(24,"Epoxyfm",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);  
 }
 
 void AliTPC::GenerNoise(Int_t tablesize)
@@ -730,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();
@@ -893,8 +946,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);
@@ -908,55 +970,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");
@@ -964,7 +1036,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++) {
@@ -984,12 +1055,6 @@ 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;
-
   }
 
   fLoader->WriteSDigits("OVERWRITE");
@@ -997,6 +1062,11 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
   if(GetDigitsArray()) delete GetDigitsArray();
   SetDigitsArray(0x0);
 
+  // cleanup storage
+  for (Int_t iRow = 0; iRow < nRowsMax; iRow++)
+    delete [] allBins[iRow];
+  delete [] allBins;
+
   return kTRUE;
 }
 
@@ -1264,14 +1334,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"); 
@@ -1331,6 +1419,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);
 
@@ -1356,7 +1463,6 @@ void AliTPC::Hits2SDigits()
   //-----------------------------------------------------------
   //   summable digits - 16 bit "ADC", no noise, no saturation
   //-----------------------------------------------------------
-  if (0) fDebugStreamer = new TTreeSRedirector("TPCSimdebug.root");
 
   if (!fTPCParam->IsGeoRead()){
     //
@@ -1375,7 +1481,7 @@ void AliTPC::Hits2SDigits()
     SetActiveSectors();
     Hits2SDigits2(iEvent);
   }
-
+  
   fLoader->UnloadHits();
   fLoader->UnloadSDigits();
   if (fDebugStreamer) {
@@ -1412,11 +1518,7 @@ void AliTPC::Hits2DigitsSector(Int_t isec)
 
   Stat_t ntracks = tH->GetEntries();
 
-  if( ntracks > 0){
 
-  //------------------------------------------- 
-  //  Only if there are any tracks...
-  //-------------------------------------------
 
     TObjArray **row;
     
@@ -1465,7 +1567,6 @@ void AliTPC::Hits2DigitsSector(Int_t isec)
       
     delete [] row; // delete the array of pointers to TObjArray-s
         
-  } // ntracks >0
 
 } // end of Hits2DigitsSector
 
@@ -1484,6 +1585,11 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
   //-----------------------------------------------------------------
  
   Float_t zerosup = fTPCParam->GetZeroSup();
+  AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor(); 
+  AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); 
+  AliTPCCalROC * gainROC = gainTPC->GetCalROC(isec);  // pad gains per given sector
+  AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(isec);  // noise per given sector
+
 
   fCurrentIndex[1]= isec;
   
@@ -1537,8 +1643,12 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
     for(Int_t ip=0;ip<nofPads;ip++){
       gi++;
       Float_t q=total(ip,it);      
-      if(fDigitsSwitch == 0){
-       q+=GetNoise();
+      if(fDigitsSwitch == 0){  
+       Float_t gain = gainROC->GetValue(irow,ip);  // get gain for given - pad-row pad 
+       Float_t noisePad = noiseROC->GetValue(irow,ip); 
+       //
+       q*=gain;
+       q+=GetNoise()*noisePad;
         if(q <=fzerosup) continue; // do not fill zeros
         q = TMath::Nint(q);
         if(q >= fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat() - 1;  // saturation
@@ -1575,7 +1685,14 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
     
     } // end of loop over time buckets
   }  // end of lop over pads 
+  //
+  // test
+  //
+  //
 
+  // glitch filters if normal simulated digits
+  //
+  if(!fDigitsSwitch) ((AliSimDigits*)dig)->GlitchFilter();
   //
   //  This row has been digitized, delete nonused stuff
   //
@@ -1625,6 +1742,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++){
@@ -1632,7 +1758,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));
@@ -1790,9 +1917,11 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
   //-----------------------------------------------------------------
   AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
   if (gAlice){ // Set correctly the magnetic field in the ExB calculation
-    AliMagF * field = gAlice->Field();
-    if (field) {
-      calib->SetExBField(field->SolenoidField());
+    if (!calib->GetExB()){
+      AliMagF * field = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField()); 
+      if (field) {
+       calib->SetExBField(field);
+      }
     }
   }