#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"
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;
wmat[0]=0.2729;
wmat[1]=0.7271;
- density=0.001754609;
+ density=1.842e-3;
AliMixture(10,"CO2",amat,zmat,density,2,wmat);
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
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
//-----------------------------------------------------------
// 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)
//
// 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
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];
+ }
+ }
//
//
//
// 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)