]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPC.cxx
Updated flags for low flux case (A. Dainese)
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
index 372e8fd2130e84ee473937767101d672dd507ec3..1b2ae45453061d36a6b8beb464d1181bc6441955 100644 (file)
 #include "AliTPCBuffer.h"
 #include "AliTPCDDLRawData.h"
 #include "AliLog.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCCalPad.h"
+#include "AliTPCCalROC.h"
+#include "AliTPCExB.h"
 #include "AliRawReader.h"
 #include "AliTPCRawStream.h"
+#include "TTreeStream.h"
 
 ClassImp(AliTPC) 
 //_____________________________________________________________________________
@@ -92,10 +97,13 @@ ClassImp(AliTPC)
                   fHitType(0),
                   fDigitsSwitch(0),
                   fSide(0),
+                   fPrimaryIonisation(0),
                   fNoiseDepth(0),
                   fNoiseTable(0),
                   fCurrentNoise(0),
-                  fActiveSectors(0)
+                  fActiveSectors(0),
+                   fGainFactor(1.),
+    fDebugStreamer(0)
 
 {
   //
@@ -123,10 +131,14 @@ AliTPC::AliTPC(const char *name, const char *title)
                   fHitType(0),
                   fDigitsSwitch(0),
                   fSide(0),
+                   fPrimaryIonisation(0),
                   fNoiseDepth(0),
                   fNoiseTable(0),
                   fCurrentNoise(0),
-                  fActiveSectors(0)
+                   fActiveSectors(0),
+    fGainFactor(1.),
+     fDebugStreamer(0)
+                  
 {
   //
   // Standard constructor
@@ -168,11 +180,13 @@ AliTPC::AliTPC(const char *name, const char *title)
 
 
   if (!strcmp(title,"Default")) {       
-    fTPCParam = new AliTPCParamSR;
+    //fTPCParam = new AliTPCParamSR;
+    fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
   } else {
     AliWarning("In Config.C you must set non-default parameters.");
     fTPCParam=0;
   }
+
 }
 
 //_____________________________________________________________________________
@@ -185,14 +199,14 @@ AliTPC::~AliTPC()
   fIshunt   = 0;
   delete fHits;
   delete fDigits;
-  delete fTPCParam;
+  //delete fTPCParam;
   delete fTrackHits; //MI 15.09.2000
   //  delete fTrackHitsOld; //MI 10.12.2001
   
   fDigitsArray = 0x0;
   delete [] fNoiseTable;
   delete [] fActiveSectors;
-
+  if (fDebugStreamer) delete fDebugStreamer;
 }
 
 //_____________________________________________________________________________
@@ -381,7 +395,7 @@ void AliTPC::CreateMaterials()
 
   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);
   //----------------------------------------------------------------------
   //               solid materials
   //----------------------------------------------------------------------
@@ -646,6 +660,57 @@ void AliTPC::CreateMaterials()
   density=7.87;
 
   AliMaterial(29,"Fe",amat[0],zmat[0],density,999.,999.);
+  //
+  // Peek - (C6H4-O-OC6H4-O-C6H4-CO)n
+  amat[0]=12.011;
+  amat[1]=1.;
+  amat[2]=15.9994;
+
+  zmat[0]=6.;
+  zmat[1]=1.;
+  zmat[2]=8.;
+
+  wmat[0]=19.;
+  wmat[1]=12.;
+  wmat[2]=3.;
+  //
+  density=1.3;
+  //
+  AliMixture(30,"Peek",amat,zmat,density,-3,wmat);  
+  //
+  //  Ceramics - Al2O3
+  //
+  amat[0] = 26.98;
+  amat[1]= 15.9994;
+
+  zmat[0] = 13.;
+  zmat[1]=8.;
+  wmat[0]=2.;
+  wmat[1]=3.;
+  density = 3.97;
+
+  AliMixture(31,"Alumina",amat,zmat,density,-2,wmat);   
+
+  //
+  // liquids
+  //
+
+  // water
+
+  amat[0]=1.;
+  amat[1]=15.9994;
+
+  zmat[0]=1.;
+  zmat[1]=8.;
+
+  wmat[0]=2.;
+  wmat[1]=1.;
+
+  density=1.;
+
+  AliMixture(32,"Water",amat,zmat,density,-2,wmat);  
  
   //----------------------------------------------------------
   // tracking media for gases
@@ -655,7 +720,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", 13, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
+  AliMedium(20, "Ne-CO2-N-3", 30, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
   //-----------------------------------------------------------  
   // tracking media for solids
   //-----------------------------------------------------------
@@ -677,7 +742,9 @@ void AliTPC::CreateMaterials()
   AliMedium(16,"G10",22,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
   AliMedium(17,"Plexiglas",27,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
   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);  
 }
 
 void AliTPC::GenerNoise(Int_t tablesize)
@@ -1143,14 +1210,21 @@ void AliTPC::SetDefaults(){
   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
   rl->CdGAFile();
   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
-  if(param){
-    AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...");
-    delete param;
-    param = new AliTPCParamSR();
-  }
-  else {
-    param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
-  }
+  // if(param){
+//     AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...");
+//     delete param;
+//     param = new AliTPCParamSR();
+//   }
+//   else {
+//     param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
+//   }
+  param = (AliTPCParamSR*)AliTPCcalibDB::Instance()->GetParameters();
+  if (!param->IsGeoRead()){
+      //
+      // read transformation matrices for gGeoManager
+      //
+      param->ReadGeoMatrices();
+    }
   if(!param){
     AliFatal("No TPC parameters found");
   }
@@ -1159,11 +1233,24 @@ void AliTPC::SetDefaults(){
   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;  
-  AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
-  rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
-  rf->SetOffset(3*param->GetZSigma());
-  rf->Update();
+
   
+  //AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
+  //rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
+  //rf->SetOffset(3*param->GetZSigma());
+  //rf->Update();
+  //
+  // Use gamma 4
+  //
+  char  strgamma4[1000];
+  sprintf(strgamma4,"AliTPCRF1D::Gamma4((x-0.135+%f)*%f,55,160)",3*param->GetZSigma(), 1000000000*param->GetTSample()/param->GetZWidth());
+  
+  TF1 * fgamma4 = new TF1("fgamma4",strgamma4, -1,1);
+  AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE,1000);
+  rf->SetParam(fgamma4,param->GetZWidth(), 1,0.2);
+  rf->SetOffset(3*param->GetZSigma()); 
+  rf->Update();
+
   TDirectory *savedir=gDirectory;
   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
   if (!f->IsOpen()) 
@@ -1221,7 +1308,6 @@ void AliTPC::Hits2Digits()
 
   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
     //PH    runLoader->GetEvent(iEvent);
-    SetActiveSectors();   
     Hits2Digits(iEvent);
   }
 
@@ -1236,6 +1322,7 @@ void AliTPC::Hits2Digits(Int_t eventnumber)
  //----------------------------------------------------
   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
   rl->GetEvent(eventnumber);
+  SetActiveSectors();   
   if (fLoader->TreeH() == 0x0) {
     if(fLoader->LoadHits()) {
       AliError("Can not load hits.");
@@ -1358,6 +1445,7 @@ void AliTPC::Hits2SDigits()
   //-----------------------------------------------------------
   //   summable digits - 16 bit "ADC", no noise, no saturation
   //-----------------------------------------------------------
+  if (0) fDebugStreamer = new TTreeSRedirector("TPCSimdebug.root");
 
   if (!fTPCParam->IsGeoRead()){
     //
@@ -1379,6 +1467,10 @@ void AliTPC::Hits2SDigits()
 
   fLoader->UnloadHits();
   fLoader->UnloadSDigits();
+  if (fDebugStreamer) {
+    delete fDebugStreamer;
+    fDebugStreamer=0;
+  }    
 }
 //_____________________________________________________________________________
 
@@ -1757,11 +1849,44 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
   // contains the track label and the position of electrons.
   //-----------------------------------------------------------------
 
+  // 
+  // The trasport of the electrons through TPC drift volume
+  //    Drift (drift velocity + velocity map(not yet implemented)))
+  //    Application of the random processes (diffusion, gas gain)
+  //    Systematic effects (ExB effect in drift volume + ROCs)  
+  //
+  // Algorithm:
+  // Loop over primary electrons:
+  //    Creation of the secondary electrons
+  //    Loop over electrons (primary+ secondaries)
+  //        Global coordinate frame:
+  //          1. Skip electrons if attached  
+  //          2. ExB effect in drift volume
+  //             a.) Simulation   calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
+  //             b.) Reconstruction -  calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
+  //          3. Generation of gas gain (Random - Exponential distribution) 
+  //          4. TransportElectron function (diffusion)
+  //
+  //        5. Conversion to the local coordinate frame  pad-row, pad, timebin
+  //        6. Apply Time0 shift - AliTPCCalPad class 
+  //            a.) Plus sign in simulation
+  //            b.) Minus sign in reconstruction 
+  // end of loop          
+  //
   //-----------------------------------------------------------------
   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
+  // Origin: Marian Ivanov,  marian.ivanov@cern.ch
   //-----------------------------------------------------------------
+  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());
+    }
+  }
 
   Float_t gasgain = fTPCParam->GetGasGain();
+  gasgain = gasgain/fGainFactor;
   Int_t i;
   Float_t xyz[5]; 
 
@@ -1877,28 +2002,81 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
       for(Int_t nel=0;nel<qI;nel++){
        // skip if electron lost due to the attachment
        if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
-       xyz[0]=tpcHit->X();
-       xyz[1]=tpcHit->Y();
-       xyz[2]=tpcHit->Z();     
+       
+       //
+       // ExB effect
+       //
+       Double_t dxyz0[3],dxyz1[3];
+       dxyz0[0]=tpcHit->X();
+       dxyz0[1]=tpcHit->Y();
+       dxyz0[2]=tpcHit->Z();   
+       if (calib->GetExB()){
+         calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
+       }else{
+         AliError("Not valid ExB calibration");
+         dxyz1[0]=tpcHit->X();
+         dxyz1[1]=tpcHit->Y();
+         dxyz1[2]=tpcHit->Z();         
+       }
+       xyz[0]=dxyz1[0];
+       xyz[1]=dxyz1[1];
+       xyz[2]=dxyz1[2];        
+       //
+       //
        //
        // protection for the nonphysical avalanche size (10**6 maximum)
        //  
        Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
        xyz[3]= (Float_t) (-gasgain*TMath::Log(rn)); 
        index[0]=1;
-       
+         
        TransportElectron(xyz,index);    
        Int_t rowNumber;
-       fTPCParam->GetPadRow(xyz,index); 
+       Int_t padrow = fTPCParam->GetPadRow(xyz,index); 
+       //
+       // Add Time0 correction due unisochronity
+       // xyz[0] - pad row coordinate 
+       // xyz[1] - pad coordinate
+       // xyz[2] - is in now time bin coordinate system
+       Float_t correction =0;
+       if (calib->GetPadTime0()){
+         if (!calib->GetPadTime0()->GetCalROC(isec)) continue;   
+         Int_t npads = fTPCParam->GetNPads(isec,padrow);
+         //      Int_t pad  = TMath::Nint(xyz[1]+fTPCParam->GetNPads(isec,TMath::Nint(xyz[0]))*0.5);
+         // pad numbering from -npads/2 .. npads/2-1
+         Int_t pad  = TMath::Nint(xyz[1]+npads/2);
+         if (pad<0) pad=0;
+         if (pad>=npads) pad=npads-1;
+         correction = calib->GetPadTime0()->GetCalROC(isec)->GetValue(padrow,pad);
+         //      printf("%d\t%d\t%d\t%f\n",isec,padrow,pad,correction);
+         if (fDebugStreamer){
+           (*fDebugStreamer)<<"Time0"<<
+             "isec="<<isec<<
+             "padrow="<<padrow<<
+             "pad="<<pad<<
+             "x0="<<xyz[0]<<
+             "x1="<<xyz[1]<<
+             "x2="<<xyz[2]<<
+             "hit.="<<tpcHit<<
+             "cor="<<correction<<
+             "\n";
+         }
+       }
+       xyz[2]+=correction;
+       xyz[2]+=fTPCParam->GetNTBinsL1();    // adding Level 1 time bin offset
+       //
        // Electron track time (for pileup simulation)
-       xyz[4] = tpcHit->Time()/fTPCParam->GetTSample();
+       xyz[2]+=tpcHit->Time()/fTPCParam->GetTSample(); // adding time of flight
+       xyz[4] =0;
+
+       //
        // row 0 - cross talk from the innermost row
        // row fNRow+1 cross talk from the outermost row
        rowNumber = index[2]+1; 
        //transform position to local digit coordinates
        //relative to nearest pad row 
        if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
-       Float_t x1,y1;
+       /*      Float_t x1,y1;
        if (isec <fTPCParam->GetNInnerSector()) {
          x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
          y1 = fTPCParam->GetYInner(rowNumber);
@@ -1910,7 +2088,7 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
        // gain inefficiency at the wires edges - linear
        x1=TMath::Abs(x1);
        y1-=1.;
-       if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.)); 
+       if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.)); */
        
        nofElectrons[rowNumber]++;        
        //----------------------------------
@@ -1968,32 +2146,6 @@ void AliTPC::Init()
   AliDebug(1,"*********************************************");
 }
 
-//_____________________________________________________________________________
-void AliTPC::MakeBranch(Option_t* option)
-{
-  //
-  // Create Tree branches for the TPC.
-  //
-  AliDebug(1,"");
-  Int_t buffersize = 4000;
-  char branchname[10];
-  sprintf(branchname,"%s",GetName());
-  
-  const char *h = strstr(option,"H");
-  
-  if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
-  
-  AliDetector::MakeBranch(option);
-  
-  const char *d = strstr(option,"D");
-  
-  if (fDigits   && fLoader->TreeD() && d) {
-    MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0);
-  }    
-
-  if (fHitType>1) MakeBranch2(option,0); // MI change 14.09.2000
-}
 //_____________________________________________________________________________
 void AliTPC::ResetDigits()
 {
@@ -2043,7 +2195,7 @@ void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
   // xyz and index must be already transformed to system 1
   //
 
-  fTPCParam->Transform1to2(xyz,index);
+  fTPCParam->Transform1to2(xyz,index);  // mis-alignment applied in this step
   
   //add diffusion
   Float_t driftl=xyz[2];
@@ -2061,7 +2213,9 @@ void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
   }
-  //add nonisochronity (not implemented yet)  
+  //add nonisochronity (not implemented yet) 
+  
 }
   
 ClassImp(AliTPChit)
@@ -2170,14 +2324,6 @@ void AliTPC::SetTreeAddress2()
 //     else
 //       AliDebug(1,"fHitType&2 Failed (can not find branch)");
 //   }
-  //set address to TREETR
-  
-  TTree *treeTR = TreeTR();
-  if (treeTR && fTrackReferences) {
-    branch = treeTR->GetBranch(GetName());
-    if (branch) branch->SetAddress(&fTrackReferences);
-  }
-
 }
 
 void AliTPC::FinishPrimary()
@@ -2545,7 +2691,15 @@ AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
     AliDebugClass(1,Form("TPC parameters %s found.",paramName));
   } else {
     AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
-    paramTPC = new AliTPCParamSR;
+    //paramTPC = new AliTPCParamSR;
+    paramTPC = AliTPCcalibDB::Instance()->GetParameters();
+    if (!paramTPC->IsGeoRead()){
+      //
+      // read transformation matrices for gGeoManager
+      //
+      paramTPC->ReadGeoMatrices();
+    }
+  
   }
   return paramTPC;