#include "AliTPCcalibDB.h"
#include "AliTPCCalPad.h"
#include "AliTPCCalROC.h"
+#include "TTreeStream.h"
using std::cout;
using std::cerr;
delete [] masks;
delete [] digarr;
}
+
+
+
+//------------------------------------------------------------------------
+void AliTPCDigitizer::DigitizeWithTailAndCrossTalk(Option_t* option, TTreeSredirector *pcstream)
+{
+ // Modified version of the digitization function
+ // Modification: adding the ion tail and crosstalk:
+ //
+ // pcstream used in order to visually inspect data
+ //
+ //
+ // Crosstalk simulation:
+ // 1.) Calculate per time bin mean charge (per pad) within anode wire segment
+ // 2.) Subsract for the clusters at given time bin fraction of (mean charge) normalized by add hoc constant
+ // AliTPCRecoParam::GetCrosstalkCorrection() (0 if not crosstalk, 1 if ideal crosstalk)
+ // for simplicity we are assuming that wire segents are related to pad-rows
+ // Wire segmentationn is obtatined from the
+ // AliTPCParam::GetWireSegment(Int_t sector, Int_t row); // to be implemented
+ // AliTPCParam::GetNPadsPerSegment(Int_t segmentID); // to be implemented
+ //
+ // Ion tail simulation:
+ // 1.) Needs signal from pad+-1, taking signal from history
+ // merge input tree's with summable digits
+ // output stored in TreeTPCD
+ //
+ char s[100];
+ char ss[100];
+ TString optionString = option;
+ if (!strcmp(optionString.Data(),"deb")) {
+ cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
+ fDebug = 3;
+ }
+ //get detector and geometry
+
+
+ AliRunLoader *rl, *orl;
+ AliLoader *gime, *ogime;
+
+ if (gAlice == 0x0)
+ {
+ Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
+ rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
+ if (rl == 0x0)
+ {
+ Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
+ return;
+ }
+ rl->LoadgAlice();
+ rl->GetAliRun();
+ }
+ AliTPC *pTPC = (AliTPC *) gAlice->GetModule("TPC");
+ AliTPCParam * param = pTPC->GetParam();
+
+ //sprintf(s,param->GetTitle());
+ snprintf(s,100,"%s",param->GetTitle());
+ //sprintf(ss,"75x40_100x60");
+ snprintf(ss,100,"75x40_100x60");
+ if(strcmp(s,ss)==0){
+ printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
+ delete param;
+ param=new AliTPCParamSR();
+ }
+ else{
+ //sprintf(ss,"75x40_100x60_150x60");
+ snprintf(ss,100,"75x40_100x60_150x60");
+ if(strcmp(s,ss)!=0) {
+ printf("No TPC parameters found...\n");
+ exit(2);
+ }
+ }
+
+ pTPC->GenerNoise(500000); //create table with noise
+ //
+ Int_t nInputs = fDigInput->GetNinputs();
+ Int_t * masks = new Int_t[nInputs];
+ for (Int_t i=0; i<nInputs;i++)
+ masks[i]= fDigInput->GetMask(i);
+ Short_t **pdig= new Short_t*[nInputs]; //pointers to the expanded digits array
+ Int_t **ptr= new Int_t*[nInputs]; //pointers to the expanded tracks array
+ Bool_t *active= new Bool_t[nInputs]; //flag for active input segments
+ Char_t phname[100];
+
+ //create digits array for given sectors
+ // make indexes
+ //
+ //create branch's in TPC treeD
+ orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
+ ogime = orl->GetLoader("TPCLoader");
+ TTree * tree = ogime->TreeD();
+ AliSimDigits * digrow = new AliSimDigits;
+
+ if (tree == 0x0)
+ {
+ ogime->MakeTree("D");
+ tree = ogime->TreeD();
+ }
+ tree->Branch("Segment","AliSimDigits",&digrow);
+ //
+ AliSimDigits ** digarr = new AliSimDigits*[nInputs];
+ for (Int_t i1=0;i1<nInputs; i1++)
+ {
+ digarr[i1]=0;
+ // intree[i1]
+ rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
+ gime = rl->GetLoader("TPCLoader");
+ gime->LoadSDigits("read");
+ TTree * treear = gime->TreeS();
+
+ if (!treear)
+ {
+ cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
+ <<" input "<< i1<<endl;
+ for (Int_t i2=0;i2<i1+1; i2++){
+
+ if(digarr[i2]) delete digarr[i2];
+ }
+ delete [] digarr;
+ delete [] active;
+ delete []masks;
+ delete []pdig;
+ delete []ptr;
+ return;
+ }
+
+ //sprintf(phname,"lhcphase%d",i1);
+ snprintf(phname,100,"lhcphase%d",i1);
+ TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
+ ->FindObject("lhcphase0");
+ if(!ph){
+ cerr<<"AliTPCDigitizer: LHC phase not found in"
+ <<" input "<< i1<<endl;
+ for (Int_t i2=0;i2<i1+1; i2++){
+ if(digarr[i2]) delete digarr[i2];
+ }
+ delete [] digarr;
+ delete [] active;
+ delete []masks;
+ delete []pdig;
+ delete []ptr;
+ return;
+ }
+ tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
+ //
+ if (treear->GetIndex()==0)
+ treear->BuildIndex("fSegmentID","fSegmentID");
+ treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
+ }
+
+
+
+
+ //
+
+ param->SetZeroSup(2);
+ Int_t zerosup = param->GetZeroSup();
+ AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
+ AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
+ //
+ // 1.) Make first loop to calculate mean amplitude per pad per segment
+ //
+ // TObjArray * crossTalkSignalArray(nsectors);
+ // TMatrixD crossTalkSignal(nWireSegments, timeBin); // structure with mean signal per pad to be filled in first loop
+ //
+ //
+ //2.) Loop over segments (padrows) of the TPC
+ //
+ for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
+ {
+ Int_t sector, padRow;
+ if (!param->AdjustSectorRow(segmentID,sector,padRow))
+ {
+ cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
+ continue;
+ }
+ AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector); // pad gains per given sector
+ AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector); // noise per given sector
+ digrow->SetID(segmentID);
+
+ Int_t nTimeBins = 0;
+ Int_t nPads = 0;
+
+ Bool_t digitize = kFALSE;
+ for (Int_t i=0;i<nInputs; i++)
+ {
+
+ rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
+ gime = rl->GetLoader("TPCLoader");
+
+ if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
+ digarr[i]->ExpandBuffer();
+ digarr[i]->ExpandTrackBuffer();
+ nTimeBins = digarr[i]->GetNRows();
+ nPads = digarr[i]->GetNCols();
+ active[i] = kTRUE;
+ if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
+ } else {
+ active[i] = kFALSE;
+ }
+ if (GetRegionOfInterest() && !digitize) break;
+ }
+ if (!digitize) continue;
+
+ digrow->Allocate(nTimeBins,nPads);
+ digrow->AllocateTrack(3);
+
+ Float_t q=0;
+ Int_t label[1000]; //stack for 300 events
+ Int_t labptr = 0;
+
+ Int_t nElems = nTimeBins*nPads;
+
+ for (Int_t i=0;i<nInputs; i++)
+ if (active[i]) {
+ pdig[i] = digarr[i]->GetDigits();
+ ptr[i] = digarr[i]->GetTracks();
+ }
+
+ Short_t *pdig1= digrow->GetDigits();
+ Int_t *ptr1= digrow->GetTracks() ;
+
+
+
+ for (Int_t elem=0;elem<nElems; elem++)
+ {
+ // padNumber=elem/nTimeBins
+ // timeBin=elem%nTimeBins
+ q=0;
+ labptr=0;
+ // looop over digits
+ for (Int_t i=0;i<nInputs; i++) if (active[i])
+ {
+ // q += digarr[i]->GetDigitFast(rows,col);
+ q += *(pdig[i]);
+
+ for (Int_t tr=0;tr<3;tr++)
+ {
+ // Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
+ Int_t lab = ptr[i][tr*nElems];
+ if ( (lab > 1) && *(pdig[i])>zerosup)
+ {
+ label[labptr]=lab+masks[i];
+ labptr++;
+ }
+ }
+ pdig[i]++;
+ ptr[i]++;
+ }
+ Int_t padNumber=elem/nTimeBins;
+
+ q/=16.; //conversion factor
+ Float_t gain = gainROC->GetValue(padRow,padNumber); // get gain for given - pad-row pad
+ //if (gain<0.5){
+ //printf("problem\n");
+ //}
+ q*= gain;
+ // Crosstalk correction:
+ // Double_t qCrossTalk=crossTalkSignal(wireSegment,timeBin); // signal matrix from per sector array
+ //
+ //
+ // Ion tail correction:
+ // padNumber=elem/nTimeBins
+ // timeBin=elem%nTimeBins
+ // elem=padNumber*nTimeBins+timeBin;
+ // lowerElem=elem-nIonTailBins;
+ // if (lowerElem<0) lowerElem=0;
+ // if (lowerElem in previospad) lowerElem = padNumber*nTimeBins;
+ //
+ // for (Int_t celem=elem-1; celem>lowerElem; celem--){
+ // Int_t deltaT=elem-celem
+ //
+ // }
+ //
+ Float_t noisePad = noiseROC->GetValue(padRow,padNumber);
+ // Float_t noise = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
+ Float_t noise = pTPC->GetNoise();
+ q+=noise*noisePad;
+ q=TMath::Nint(q); // round to the nearest integer
+ /*
+ fill infor to check consistency of the data
+ */
+ if (q > zerosup)
+ {
+ if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
+ //digrow->SetDigitFast((Short_t)q,rows,col);
+ *pdig1 =Short_t(q);
+ for (Int_t tr=0;tr<3;tr++)
+ {
+ if (tr<labptr)
+ ptr1[tr*nElems] = label[tr];
+ }
+ }
+ pdig1++;
+ ptr1++;
+ }
+ //
+ // glitch filter
+ //
+ digrow->GlitchFilter();
+ //
+ digrow->CompresBuffer(1,zerosup);
+ digrow->CompresTrackBuffer(1);
+ tree->Fill();
+ if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
+ } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
+
+
+ orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
+ ogime = orl->GetLoader("TPCLoader");
+ ogime->WriteDigits("OVERWRITE");
+
+ //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
+
+ delete digrow;
+ for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
+ delete []masks;
+ delete []pdig;
+ delete []ptr;
+ delete []active;
+ delete []digarr;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff