]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPC.cxx
closing files properly
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
index 11b1bc1f73ed9803abb5d823c1406a2a892c43b9..f08e3fb4cab032095fcfd3420e5cf6e6681046cf 100644 (file)
 #include "AliTPCTrackHitsV2.h"
 #include "AliTrackReference.h"
 #include "AliMC.h"
+#include "AliStack.h"
 #include "AliTPCDigitizer.h"
 #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"
 
 ClassImp(AliTPC) 
 //_____________________________________________________________________________
-AliTPC::AliTPC()
+  AliTPC::AliTPC():AliDetector(),
+                  fDefaults(0),
+                  fSens(0),
+                  fNsectors(0),
+                  fDigitsArray(0),
+                  fTPCParam(0),
+                  fTrackHits(0),
+                  fHitType(0),
+                  fDigitsSwitch(0),
+                  fSide(0),
+                   fPrimaryIonisation(0),
+                  fNoiseDepth(0),
+                  fNoiseTable(0),
+                  fCurrentNoise(0),
+                  fActiveSectors(0),
+                   fGainFactor(1.)
+                   
+
 {
   //
   // Default constructor
   //
   fIshunt   = 0;
-  fHits     = 0;
-  fDigits   = 0;
-  fNsectors = 0;
-  fDigitsArray = 0;
-  fDefaults = 0;
-  fTrackHits = 0; 
   //  fTrackHitsOld = 0;   
 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
   fHitType = 4; // ROOT containers
 #else
   fHitType = 2; //default CONTAINERS - based on ROOT structure
 #endif 
-  fTPCParam = 0;    
-  fNoiseTable = 0;
-  fActiveSectors =0;
-  fSens = 0;
-
 }
  
 //_____________________________________________________________________________
 AliTPC::AliTPC(const char *name, const char *title)
-      : AliDetector(name,title)
+  : AliDetector(name,title),
+                   fDefaults(0),
+                  fSens(0),
+                  fNsectors(0),
+                  fDigitsArray(0),
+                  fTPCParam(0),
+                  fTrackHits(0),
+                  fHitType(0),
+                  fDigitsSwitch(0),
+                  fSide(0),
+                   fPrimaryIonisation(0),
+                  fNoiseDepth(0),
+                  fNoiseTable(0),
+                  fCurrentNoise(0),
+                   fActiveSectors(0),
+                   fGainFactor(1.)
+                  
 {
   //
   // Standard constructor
@@ -118,8 +146,6 @@ AliTPC::AliTPC(const char *name, const char *title)
   // Initialise arrays of hits and digits 
   fHits     = new TClonesArray("AliTPChit",  176);
   gAlice->GetMCApp()->AddHitList(fHits); 
-  fDigitsArray = 0;
-  fDefaults = 0;
   //
   fTrackHits = new AliTPCTrackHitsV2;  
   fTrackHits->SetHitPrecision(0.002);
@@ -131,23 +157,20 @@ AliTPC::AliTPC(const char *name, const char *title)
   //fTrackHitsOld->SetStepPrecision(0.003);  
   //fTrackHitsOld->SetMaxDistance(100); 
 
-  fNoiseTable =0;
 
 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
   fHitType = 4; // ROOT containers
 #else
   fHitType = 2;
 #endif
-  fActiveSectors = 0;
-  //
-  // Initialise counters
-  fNsectors = 0;
+
+
 
   //
   fIshunt     =  0;
   //
   // Initialise color attributes
-  SetMarkerColor(kYellow);
+  //PH SetMarkerColor(kYellow);
 
   //
   //  Set TPC parameters
@@ -160,16 +183,9 @@ AliTPC::AliTPC(const char *name, const char *title)
     AliWarning("In Config.C you must set non-default parameters.");
     fTPCParam=0;
   }
-  fSens = 0;
-
 }
 
 //_____________________________________________________________________________
-AliTPC::AliTPC(const AliTPC& t):AliDetector(t){
-  //
-  // dummy copy constructor
-  //
-}
 AliTPC::~AliTPC()
 {
   //
@@ -375,7 +391,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
   //----------------------------------------------------------------------
@@ -649,7 +665,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);
   //-----------------------------------------------------------  
   // tracking media for solids
   //-----------------------------------------------------------
@@ -716,8 +732,7 @@ void    AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
 {
   // activate interesting sectors
   SetTreeAddress();//just for security
-  if (fActiveSectors) delete [] fActiveSectors;
-  fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
+  if (!fActiveSectors) fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
   for (Int_t i=0;i<n;i++) 
     if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector())  fActiveSectors[sectors[i]]=kTRUE;
@@ -731,8 +746,7 @@ void    AliTPC::SetActiveSectors(Int_t flag)
   //loop over tracks
   SetTreeAddress();//just for security
   if (fHitType==0) return;  // if Clones hit - not short volume ID information
-  if (fActiveSectors) delete [] fActiveSectors;
-  fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
+  if (!fActiveSectors) fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
   if (flag) {
     for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
     return;
@@ -759,8 +773,17 @@ void    AliTPC::SetActiveSectors(Int_t flag)
       br1->GetEvent(track);
       br2->GetEvent(track);
       Int_t *volumes = fTrackHits->GetVolumes();
-      for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++)
-       fActiveSectors[volumes[j]]=kTRUE;
+      for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++) {
+       if (volumes[j]>-1 && volumes[j]<fTPCParam->GetNSector()) {
+         fActiveSectors[volumes[j]]=kTRUE;
+       }
+       else {
+           AliError(Form("Volume %d -> sector number %d is outside (0..%d)",
+                         j,
+                         volumes[j],
+                         fTPCParam->GetNSector()));
+       }
+      }
     }
     
     //
@@ -897,7 +920,7 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
   fTPCParam->SetZeroSup(0);
 
   // Loop over sectors
-  const Int_t maxTime = fTPCParam->GetMaxTBin();
+  const Int_t kmaxTime = fTPCParam->GetMaxTBin();
   const Int_t kNIS = fTPCParam->GetNInnerSector();
   const Int_t kNOS = fTPCParam->GetNOuterSector();
   const Int_t kNS = kNIS + kNOS;
@@ -926,7 +949,7 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
     allBins = new Short_t*[nRows];
     for (Int_t iRow = 0; iRow < nRows; iRow++) {
       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
-      Int_t maxBin = maxTime*maxPad;
+      Int_t maxBin = kmaxTime*maxPad;
       allBins[iRow] = new Short_t[maxBin];
       memset(allBins[iRow],0,sizeof(Short_t)*maxBin);
     }
@@ -950,19 +973,19 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
                      iPad, 0, maxPad -1));
 
       Int_t iTimeBin = input.GetTime();
-      if ( iTimeBin < 0 || iTimeBin >= maxTime)
+      if ( iTimeBin < 0 || iTimeBin >= kmaxTime)
        AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
-                     iTimeBin, 0, maxTime -1));
+                     iTimeBin, 0, kmaxTime -1));
       
-      Int_t maxBin = maxTime*maxPad;
+      Int_t maxBin = kmaxTime*maxPad;
 
-      if (((iPad*maxTime+iTimeBin) >= maxBin) ||
-         ((iPad*maxTime+iTimeBin) < 0))
+      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*maxTime+iTimeBin] = input.GetSignal();
+      allBins[iRow][iPad*kmaxTime+iTimeBin] = input.GetSignal();
 
     } // End loop over altro data
     
@@ -976,8 +999,8 @@ Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
 
       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
       for(Int_t iPad = 0; iPad < maxPad; iPad++) {
-       for(Int_t iTimeBin = 0; iTimeBin < maxTime; iTimeBin++) {
-         Short_t q = allBins[iRow][iPad*maxTime + iTimeBin];
+       for(Int_t iTimeBin = 0; iTimeBin < kmaxTime; iTimeBin++) {
+         Short_t q = allBins[iRow][iPad*kmaxTime + iTimeBin];
          if (q <= 0) continue;
          q *= 16;
          dig->SetDigitFast((Short_t)q,iTimeBin,iPad);
@@ -1207,8 +1230,7 @@ void AliTPC::Hits2Digits()
   AliRunLoader* runLoader = fLoader->GetRunLoader(); 
 
   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
-    runLoader->GetEvent(iEvent);
-    SetActiveSectors();   
+    //PH    runLoader->GetEvent(iEvent);
     Hits2Digits(iEvent);
   }
 
@@ -1223,6 +1245,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.");
@@ -1744,11 +1767,38 @@ 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();
 
   Float_t gasgain = fTPCParam->GetGasGain();
+  gasgain = gasgain/fGainFactor;
   Int_t i;
   Float_t xyz[5]; 
 
@@ -1809,7 +1859,7 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
       }
 
       // Remove hits which arrive before the TPC opening gate signal
-      if(((fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
+      if(((fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
          /fTPCParam->GetDriftV()+tpcHit->Time())<fTPCParam->GetGateDelay()) {
        tpcHit = (AliTPChit*) NextHit();
        continue;
@@ -1850,7 +1900,7 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
       //---------------------------------------------------
 
 
-      Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
+      Float_t time = 1.e6*(fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
        /fTPCParam->GetDriftV(); 
       // in microseconds!      
       Float_t attProb = fTPCParam->GetAttCoef()*
@@ -1864,19 +1914,49 @@ 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); 
+       //
+       // 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;
+         correction = calib->GetPadTime0()->GetCalROC(isec)->GetValue(TMath::Nint(xyz[0]),TMath::Nint(xyz[1]));
+       }
+       xyz[2]+=correction;
+       //
        // Electron track time (for pileup simulation)
        xyz[4] = tpcHit->Time()/fTPCParam->GetTSample();
        // row 0 - cross talk from the innermost row
@@ -2048,14 +2128,32 @@ 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)
+  //______________________________________________________________________
+  AliTPChit::AliTPChit()
+            :AliHit(),
+            fSector(0),
+            fPadRow(0),
+            fQ(0),
+            fTime(0)
+{
+  //
+  // default
+  //
+
+}
 //_____________________________________________________________________________
-AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
-AliHit(shunt,track)
+AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits)
+          :AliHit(shunt,track),
+            fSector(0),
+            fPadRow(0),
+            fQ(0),
+            fTime(0)
 {
   //
   // Creates a TPC hit object
@@ -2141,14 +2239,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()
@@ -2283,7 +2373,7 @@ void AliTPC::RemapTrackHitIDs(Int_t *map)
     TClonesArray * arr = fTrackHits->GetArray();;
     for (Int_t i=0;i<arr->GetEntriesFast();i++){
       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
-      info->fTrackID = map[info->fTrackID];
+      info->SetTrackID(map[info->GetTrackID()]);
     }
   }
 }
@@ -2393,11 +2483,11 @@ void AliTPC::LoadPoints2(Int_t)
   for(trk=0; trk<tracks; ++trk) {
     if(ntrk[trk]) {
       points = new AliPoints();
-      points->SetMarkerColor(GetMarkerColor());
-      points->SetMarkerSize(GetMarkerSize());
+      points->SetMarkerColor(kYellow); //PH kYellow it the default in TPC
+      points->SetMarkerSize(1);//PH Default size=1
       points->SetDetector(this);
       points->SetParticle(trk);
-      points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
+      points->SetPolyMarker(ntrk[trk],coor[trk],1); // Default style=1
       fPoints->AddAt(points,trk);
       delete [] coor[trk];
       coor[trk]=0;
@@ -2476,7 +2566,7 @@ void AliTPC::LoadPoints3(Int_t)
   for(trk=0; trk<tracks; ++trk) {
     if(ntrk[trk]) {
       points = new AliPoints();
-      points->SetMarkerColor(GetMarkerColor()+1);
+      points->SetMarkerColor(kMagenta); //PH kYellow + 1 is kMagenta
       points->SetMarkerStyle(5);
       points->SetMarkerSize(0.2);
       points->SetDetector(this);