]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/Sim/AliTPC.cxx
Merge branch 'master', remote branch 'origin' into TPCdev
[u/mrichter/AliRoot.git] / TPC / Sim / AliTPC.cxx
index 4603f8fbc443620d458d529f31de95b7bbeaf9ae..7fa8acbe00b96e3d95a0c4177af003a6b34db71d 100644 (file)
@@ -55,6 +55,8 @@
 #include <TParameter.h>
 
 #include "AliDigits.h"
+#include "AliHeader.h"
+
 #include "AliMagF.h"
 #include "AliRun.h"
 #include "AliRunLoader.h"
 #include "AliTPCDDLRawData.h"
 #include "AliLog.h"
 #include "AliTPCcalibDB.h"
+#include "AliTPCRecoParam.h"
 #include "AliTPCCalPad.h"
 #include "AliTPCCalROC.h"
 #include "AliTPCExB.h"
+#include "AliTPCCorrection.h"
 #include "AliCTPTimeParams.h"
 #include "AliRawReader.h"
 #include "AliTPCRawStreamV3.h"
@@ -245,9 +249,9 @@ void AliTPC::CreateMaterials()
    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
-  Float_t wmat[5]; // proportions
+  Float_t amat[7]; // atomic numbers
+  Float_t zmat[7]; // z
+  Float_t wmat[7]; // proportions
 
   Float_t density;
  
@@ -271,7 +275,7 @@ void AliTPC::CreateMaterials()
   wmat[0]=0.2729;
   wmat[1]=0.7271;
 
-  density=0.001754609;
+  density=1.842e-3;
 
 
   AliMixture(10,"CO2",amat,zmat,density,2,wmat);
@@ -292,52 +296,107 @@ void AliTPC::CreateMaterials()
   AliMixture(11,"Air",amat,zmat,density,2,wmat);
   
   //----------------------------------------------------------------
-  // drift gases 
+  // drift gases 5 mixtures, 5 materials
   //----------------------------------------------------------------
-
   //
   // Drift gases 1 - nonsensitive, 2 - sensitive, 3 - for Kr
   //  Composition by % of volume, values at 20deg and 1 atm.
   //
   //  get the geometry title - defined in Config.C
   //
-
-  TString title(GetTitle());
-  
+  //--------------------------------------------------------------
+  //  predefined gases, composition taken from param file
+  //--------------------------------------------------------------
+  TString names[6]={"Ne","Ar","CO2","N","CF4","CH4"};
+  TString gname;
+  Float_t *comp = fTPCParam->GetComposition();
+  // indices:
+  // 0-Ne, 1-Ar, 2-CO2, 3-N, 4-CF4, 5-CH4
   //
-  amat[0]= 20.18;
-  amat[1]=12.011;
-  amat[2]=15.9994;
-
-
-  zmat[0]= 10.; 
-  zmat[1]=6.;
-  zmat[2]=8.;
-
-  if(title == TString("Ne-CO2") || title == TString("Default")){
-    wmat[0]=0.8038965;
-    wmat[1]= 0.053519;
-    wmat[2]= 0.1425743;
-    density=0.0009393;
-    //
-    AliMixture(12,"Ne-CO2-1",amat,zmat,density,3,wmat);
-    AliMixture(13,"Ne-CO2-2",amat,zmat,density,3,wmat);
-    AliMixture(35,"Ne-CO2-3",amat,zmat,density,3,wmat);
+  // elements' masses
+  // 
+  amat[0]=20.18; //Ne
+  amat[1]=39.95; //Ar
+  amat[2]=12.011; //C
+  amat[3]=15.9994; //O
+  amat[4]=14.007; //N
+  amat[5]=18.998; //F
+  amat[6]=1.; //H
+  //
+  // elements' atomic numbers
+  //
+  // 
+  zmat[0]=10.; //Ne
+  zmat[1]=18.; //Ar
+  zmat[2]=6.;  //C
+  zmat[3]=8.;  //O
+  zmat[4]=7.;  //N
+  zmat[5]=9.;  //F
+  zmat[6]=1.;  //H
+  //
+  // Mol masses
+  //
+  Float_t wmol[6];
+  wmol[0]=20.18; //Ne
+  wmol[1]=39.948; //Ar
+  wmol[2]=44.0098; //CO2
+  wmol[3]=2.*14.0067; //N2
+  wmol[4]=88.0046; //CF4
+  wmol[5]=16.011; //CH4
+  //
+  Float_t wtot=0.; //total mass of the mixture
+  for(Int_t i =0;i<6;i++){
+    wtot += *(comp+i)*wmol[i]; 
   }
-  else if (title == TString("Ne-CO2-N")){
-     amat[3]=14.007;
-     zmat[3]=7.; 
-     wmat[0]=0.756992632;
-     wmat[1]=0.056235789; 
-     wmat[2]=0.128469474; 
-     wmat[3]=0.058395789;
-     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);
-  
+  wmat[0]=comp[0]*amat[0]/wtot; //Ne
+  wmat[1]=comp[1]*amat[1]/wtot; //Ar
+  wmat[2]=(comp[2]*amat[2]+comp[4]*amat[2]+comp[5]*amat[2])/wtot; //C
+  wmat[3]=comp[2]*amat[3]*2./wtot; //O
+  wmat[4]=comp[3]*amat[4]*2./wtot; //N
+  wmat[5]=comp[4]*amat[5]*4./wtot; //F
+  wmat[6]=comp[5]*amat[6]*4./wtot; //H
+  //
+  // densities (NTP)
+  //
+  Float_t dens[6]={0.839e-3,1.661e-3,1.842e-3,1.251e-3,3.466e-3,0.668e-3};
+ //
+  density=0.;
+  for(Int_t i=0;i<6;i++){
+    density += comp[i]*dens[i];
   }
+  //
+  // names
+  //
+  Int_t cnt=0;
+  for(Int_t i =0;i<6;i++){
+    if(comp[i]){
+     if(cnt)gname+="-"; 
+      gname+=names[i];
+      cnt++;   
+    } 
+  }
+TString gname1,gname2,gname3;
+gname1 = gname + "-1";
+gname2 = gname + "-2";
+gname3 = gname + "-3";
+//
+// take only elements with nonzero weights
+//
+ Float_t amat1[6],zmat1[6],wmat1[6];
+ cnt=0;
+ for(Int_t i=0;i<7;i++){
+   if(wmat[i]){
+     zmat1[cnt]=zmat[i];
+     amat1[cnt]=amat[i];
+     wmat1[cnt]=wmat[i];
+     cnt++;
+   }
+ }
+   
+//
+AliMixture(12,gname1.Data(),amat1,zmat1,density,cnt,wmat1); // nonsensitive
+AliMixture(13,gname2.Data(),amat1,zmat1,density,cnt,wmat1); // sensitive
+AliMixture(40,gname3.Data(),amat1,zmat1,density,cnt,wmat1); //sensitive Kr
 
 
 
@@ -742,7 +801,7 @@ void AliTPC::CreateMaterials()
   AliMedium(1, "DriftGas1", 12, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
   AliMedium(2, "DriftGas2", 13, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001); 
-  AliMedium(20, "DriftGas3", 35, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
+  AliMedium(20, "DriftGas3", 40, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
   //-----------------------------------------------------------  
   // tracking media for solids
   //-----------------------------------------------------------
@@ -2003,8 +2062,11 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
   //        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);
+  //             a.) Simulation   calib->GetExB()->CorrectInverse(dxyz0,dxyz1); (applied on the el. level)
+  //             b.) Reconstruction -  calib->GetExB()->Correct(dxyz0,dxyz1);   (applied on the space point)
+  //             changed to the full distrotion (not only due to the ExB)
+  //             aNew.)  correction->DistortPoint(distPoint, sector);
+  //             bNew.)  correction->CorrectPoint(distPoint, sector);
   //          3. Generation of gas gain (Random - Exponential distribution) 
   //          4. TransportElectron function (diffusion)
   //
@@ -2019,19 +2081,42 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
   // Origin: Marian Ivanov,  marian.ivanov@cern.ch
   //-----------------------------------------------------------------
   AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
-  if (gAlice){ // Set correctly the magnetic field in the ExB calculation
-    if (!calib->GetExB()){
-      AliMagF * field = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField()); 
-      if (field) {
-       calib->SetExBField(field);
+
+  AliTPCCorrection * correctionDist = calib->GetTPCComposedCorrection();  
+
+  AliTPCRecoParam *tpcrecoparam = calib->GetRecoParam(0); //FIXME: event specie should not be set by hand, However the parameters read here are the same for al species
+
+//  AliWarning(Form("Flag for ExB correction \t%d",tpcrecoparam->GetUseExBCorrection())); 
+//  AliWarning(Form("Flag for Composed correction \t%d",calib->GetRecoParam()->GetUseComposedCorrection()));
+
+  if (tpcrecoparam->GetUseExBCorrection()) {
+    if (gAlice){ // Set correctly the magnetic field in the ExB calculation
+      if (!calib->GetExB()){
+        AliMagF * field = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField()); 
+        if (field) {
+         calib->SetExBField(field);
+        }
       }
     }
+  } else if (tpcrecoparam->GetUseComposedCorrection()) {
+    AliMagF * field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); 
+    Double_t bzpos[3]={0,0,0};
+    if (!correctionDist) correctionDist = calib->GetTPCComposedCorrection(field->GetBz(bzpos));
+
+    if (!correctionDist){
+      AliFatal("Correction map does not exist. Check the OCDB or your setup");
+    }
   }
 
   Float_t gasgain = fTPCParam->GetGasGain();
   gasgain = gasgain/fGainFactor;
+  //  const Int_t timeStamp = 1; //where to get it? runloader->GetHeader()->GetTimeStamp(). https://savannah.cern.ch/bugs/?53025
+  const Int_t timeStamp = fLoader->GetRunLoader()->GetHeader()->GetTimeStamp(); //?
+  Double_t correctionHVandPT = calib->GetGainCorrectionHVandPT(timeStamp, calib->GetRun(), isec, 5 ,tpcrecoparam->GetGainCorrectionHVandPTMode());
+  gasgain*=correctionHVandPT;
+
   Int_t i;
-  Float_t xyz[5]
+  Float_t xyz[5]={0,0,0,0,0};
 
   AliTPChit *tpcHit; // pointer to a sigle TPC hit    
   //MI change
@@ -2145,25 +2230,39 @@ 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!
-       
+       // use default hit position
+       xyz[0]=tpcHit->X();
+       xyz[1]=tpcHit->Y();
+       xyz[2]=tpcHit->Z(); 
        //
-       // ExB effect
+       // ExB effect - distort hig if specifiend in the RecoParam
        //
-       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];        
+        if (tpcrecoparam->GetUseExBCorrection()) {
+         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];      
+        } else if (tpcrecoparam->GetUseComposedCorrection()) {
+        //      Use combined correction/distortion  class AliTPCCorrection
+          if (correctionDist){
+            Float_t distPoint[3]={tpcHit->X(),tpcHit->Y(), tpcHit->Z()};
+            correctionDist->DistortPoint(distPoint, isec);
+            xyz[0]=distPoint[0];
+            xyz[1]=distPoint[1];
+            xyz[2]=distPoint[2];
+           }      
+        }
        //
        //
        //
@@ -2182,41 +2281,43 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
        // 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";
+       if (tpcrecoparam->GetUseExBCorrection()) {
+          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";
+           }
          }
-       }
+        }
         if (AliTPCcalibDB::Instance()->IsTrgL0()){  
-         // Modification 14.03
-         // distinguish between the L0 and L1 trigger as it is done in the reconstruction
-         // by defualt we assume L1 trigger is used - make a correction in case of  L0
-         AliCTPTimeParams* ctp = AliTPCcalibDB::Instance()->GetCTPTimeParams();
-         if (ctp){
-           //for TPC standalone runs no ctp info
-           Double_t delay = ctp->GetDelayL1L0()*0.000000025;
-           xyz[2]+=delay/fTPCParam->GetTSample();  // adding the delay (in the AliTPCTramsform opposite sign)
-         }
-       }
-       xyz[2]+=correction;
+          // Modification 14.03
+          // distinguish between the L0 and L1 trigger as it is done in the reconstruction
+          // by defualt we assume L1 trigger is used - make a correction in case of  L0
+          AliCTPTimeParams* ctp = AliTPCcalibDB::Instance()->GetCTPTimeParams();
+          if (ctp){
+            //for TPC standalone runs no ctp info
+            Double_t delay = ctp->GetDelayL1L0()*0.000000025;
+            xyz[2]+=delay/fTPCParam->GetTSample();  // adding the delay (in the AliTPCTramsform opposite sign)
+          }
+        }
+       if (tpcrecoparam->GetUseExBCorrection()) xyz[2]+=correction; // In Correction there is already a corretion for the time 0 offset so not needed
        xyz[2]+=fTPCParam->GetNTBinsL1();    // adding Level 1 time bin offset
        //
        // Electron track time (for pileup simulation)