/* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */
//-----------------------------------------------------------------
-// Implementation of the TPC cluster class
+// Implementation of the TPC Kr cluster class
//
// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl
//-----------------------------------------------------------------
+/*
+Instruction - how to use that:
+There are two macros prepared. One is for preparing clusters from MC
+samples: FindKrClusters.C. The output is kept in TPC.RecPoints.root.
+The other macro is prepared for data analysis: FindKrClustersRaw.C.
+The output is created for each processed file in root file named adc.root.
+For each data subsample the same named file is created. So be careful
+do not overwrite them.
+
+Additional selection criteria to select the GOLD cluster
+Example:
+// open file with clusters
+TFile f("Krypton.root");
+TTree * tree = (TTree*)f.Get("Kr")
+TCut cutR0("cutR0","fADCcluster/fSize<100"); // adjust it according v seetings -
+TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal
+TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.2"); // digital noise removal
+TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal
+TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
+TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;
+This values are typical values to be applied in selectors
+
+
+*
+**** MC ****
+*
+
+To run clusterizaton for MC type:
+.x FindKrClusters.C
+
+If you don't want to use the standard selection criteria then you
+have to do following:
+
+// load the standard setup
+AliRunLoader* rl = AliRunLoader::Open("galice.root");
+AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
+tpcl->LoadDigits();
+rl->LoadgAlice();
+gAlice=rl->GetAliRun();
+TDirectory *cwd = gDirectory;
+AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");
+Int_t ver = tpc->IsVersion();
+rl->CdGAFile();
+AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
+AliTPCDigitsArray *digarr=new AliTPCDigitsArray;
+digarr->Setup(param);
+cwd->cd();
+
+//loop over events
+Int_t nevmax=rl->GetNumberOfEvents();
+for(Int_t nev=0;nev<nevmax ;nev++){
+ rl->GetEvent(nev);
+ TTree* input_tree= tpcl->TreeD();//load tree with digits
+ digarr->ConnectTree(input_tree);
+ TTree *output_tree =tpcl->TreeR();//load output tree
+
+ AliTPCclustererKr *clusters = new AliTPCclustererKr();
+ clusters->SetParam(param);
+ clusters->SetInput(input_tree);
+ clusters->SetOutput(output_tree);
+ clusters->SetDigArr(digarr);
+
+//If you want to change the cluster finder parameters for MC there are
+//several of them:
+
+//1. signal threshold (everything below the given number is treated as 0)
+ clusters->SetMinAdc(3);
+
+//2. number of neighbouring timebins to be considered
+ clusters->SetMinTimeBins(2);
+
+//3. distance of the cluster center to the center of a pad in pad-padrow plane
+//(in cm). Remenber that this is still quantified by pad size.
+ clusters->SetMaxPadRangeCm(2.5);
+
+//4. distance of the cluster center to the center of a padrow in pad-padrow
+//plane (in cm). Remenber that this is still quantified by pad size.
+ clusters->SetMaxRowRangeCm(3.5);
+
+//5. distance of the cluster center to the max time bin on a pad (in tackts)
+//ie. fabs(centerT - time)<7
+ clusters->SetMaxTimeRange(7);
+
+//6. cut reduce peak at 0. There are noises which appear mostly as two
+//timebins on one pad.
+ clusters->SetValueToSize(3.5);
+
+
+ clusters->finderIO();
+ tpcl->WriteRecPoints("OVERWRITE");
+}
+delete rl;//cleans everything
+
+*
+********* DATA *********
+*
+
+To run clusterizaton for DATA for file named raw_data.root type:
+.x FindKrClustersRaw.C("raw_data.root")
+
+If you want to change some criteria do the following:
+
+//
+// remove Altro warnings
+//
+AliLog::SetClassDebugLevel("AliTPCRawStream",-5);
+AliLog::SetClassDebugLevel("AliRawReaderDate",-5);
+AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);
+AliLog::SetModuleDebugLevel("RAW",-5);
+
+//
+// Get database with noises
+//
+// char *ocdbpath = gSystem->Getenv("OCDB_PATH");
+char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
+if (ocdbpath==0){
+ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
+}
+AliCDBManager * man = AliCDBManager::Instance();
+man->SetDefaultStorage(ocdbpath);
+man->SetRun(0);
+AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
+AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
+
+//define tree
+TFile *hfile=new TFile("adc.root","RECREATE","ADC file");
+// Create a ROOT Tree
+TTree *mytree = new TTree("Kr","Krypton cluster tree");
+
+//define infput file
+const char *fileName="data.root";
+AliRawReader *reader = new AliRawReaderRoot(fileName);
+//AliRawReader *reader = new AliRawReaderDate(fileName);
+reader->Reset();
+AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);
+stream->SelectRawData("TPC");
+
+//one general output
+AliTPCclustererKr *clusters = new AliTPCclustererKr();
+clusters->SetOutput(mytree);
+clusters->SetRecoParam(0);//standard reco parameters
+AliTPCParamSR *param=new AliTPCParamSR();
+clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.)
+
+//set cluster finder parameters (from data):
+//1. zero suppression parameter
+ clusters->SetZeroSup(param->GetZeroSup());
+
+//2. first bin
+ clusters->SetFirstBin(60);
+
+//3. last bin
+ clusters->SetLastBin(950);
+
+//4. maximal noise
+ clusters->SetMaxNoiseAbs(2);
+
+//5. maximal amount of sigma of noise
+ clusters->SetMaxNoiseSigma(3);
+
+//The remaining parameters are the same paramters as for MC (see MC section
+//points 1-6)
+ clusters->SetMinAdc(3);
+ clusters->SetMinTimeBins(2);
+ clusters->SetMaxPadRangeCm(2.5);
+ clusters->SetMaxRowRangeCm(3.5);
+ clusters->SetMaxTimeRange(7);
+ clusters->SetValueToSize(3.5);
+
+while (reader->NextEvent()) {
+ clusters->FinderIO(reader);
+}
+
+hfile->Write();
+hfile->Close();
+delete stream;
+
+
+*/
+
#include "AliTPCclustererKr.h"
#include "AliTPCclusterKr.h"
-#include <vector>
+//#include <vector>
+#include <list>
#include "TObject.h"
#include "AliPadMax.h"
#include "AliSimDigits.h"
#include "AliTPCvtpr.h"
#include "AliTPCClustersRow.h"
#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TTreeStream.h"
+
+#include "AliTPCTransform.h"
//used in raw data finder
#include "AliTPCROC.h"
AliTPCclustererKr::AliTPCclustererKr()
:TObject(),
fRawData(kFALSE),
- fRowCl(0),
fInput(0),
fOutput(0),
fParam(0),
fDigarr(0),
fRecoParam(0),
- fIsOldRCUFormat(kFALSE)
+ fZeroSup(2),
+ fFirstBin(60),
+ fLastBin(950),
+ fMaxNoiseAbs(2),
+ fMaxNoiseSigma(3),
+ fMinAdc(3),
+ fMinTimeBins(2),
+// fMaxPadRange(4),
+// fMaxRowRange(3),
+ fMaxTimeRange(7),
+ fValueToSize(3.5),
+ fMaxPadRangeCm(2.5),
+ fMaxRowRangeCm(3.5),
+ fDebugLevel(-1),
+ fHistoRow(0),
+ fHistoPad(0),
+ fHistoTime(0),
+ fHistoRowPad(0)
{
//
// default constructor
AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m)
:TObject(),
fRawData(kFALSE),
- fRowCl(0),
fInput(0),
fOutput(0),
fParam(0),
fDigarr(0),
fRecoParam(0),
- fIsOldRCUFormat(kFALSE)
+ fZeroSup(2),
+ fFirstBin(60),
+ fLastBin(950),
+ fMaxNoiseAbs(2),
+ fMaxNoiseSigma(3),
+ fMinAdc(3),
+ fMinTimeBins(2),
+// fMaxPadRange(4),
+// fMaxRowRange(3),
+ fMaxTimeRange(7),
+ fValueToSize(3.5),
+ fMaxPadRangeCm(2.5),
+ fMaxRowRangeCm(3.5),
+ fDebugLevel(-1),
+ fHistoRow(0),
+ fHistoPad(0),
+ fHistoTime(0),
+ fHistoRowPad(0)
{
//
// copy constructor
//
fParam = param.fParam;
- fRecoParam=param.fRecoParam;
- fIsOldRCUFormat=param.fIsOldRCUFormat;
- fRawData=param.fRawData;
- fRowCl =param.fRowCl ;
- fInput =param.fInput ;
- fOutput=param.fOutput;
- fDigarr=param.fDigarr;
+ fRecoParam = param.fRecoParam;
+ fRawData = param.fRawData;
+ fInput = param.fInput ;
+ fOutput = param.fOutput;
+ fDigarr = param.fDigarr;
+ fZeroSup = param.fZeroSup ;
+ fFirstBin = param.fFirstBin ;
+ fLastBin = param.fLastBin ;
+ fMaxNoiseAbs = param.fMaxNoiseAbs ;
+ fMaxNoiseSigma = param.fMaxNoiseSigma ;
+ fMinAdc = param.fMinAdc;
+ fMinTimeBins = param.fMinTimeBins;
+// fMaxPadRange = param.fMaxPadRange ;
+// fMaxRowRange = param.fMaxRowRange ;
+ fMaxTimeRange = param.fMaxTimeRange;
+ fValueToSize = param.fValueToSize;
+ fMaxPadRangeCm = param.fMaxPadRangeCm;
+ fMaxRowRangeCm = param.fMaxRowRangeCm;
+ fDebugLevel = param.fDebugLevel;
+ fHistoRow = param.fHistoRow ;
+ fHistoPad = param.fHistoPad ;
+ fHistoTime = param.fHistoTime;
+ fHistoRowPad = param.fHistoRowPad;
+
}
AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
{
+ //
+ // assignment operator
+ //
fParam = param.fParam;
- fRecoParam=param.fRecoParam;
- fIsOldRCUFormat=param.fIsOldRCUFormat;
- fRawData=param.fRawData;
- fRowCl =param.fRowCl ;
- fInput =param.fInput ;
- fOutput=param.fOutput;
- fDigarr=param.fDigarr;
+ fRecoParam = param.fRecoParam;
+ fRawData = param.fRawData;
+ fInput = param.fInput ;
+ fOutput = param.fOutput;
+ fDigarr = param.fDigarr;
+ fZeroSup = param.fZeroSup ;
+ fFirstBin = param.fFirstBin ;
+ fLastBin = param.fLastBin ;
+ fMaxNoiseAbs = param.fMaxNoiseAbs ;
+ fMaxNoiseSigma = param.fMaxNoiseSigma ;
+ fMinAdc = param.fMinAdc;
+ fMinTimeBins = param.fMinTimeBins;
+// fMaxPadRange = param.fMaxPadRange ;
+// fMaxRowRange = param.fMaxRowRange ;
+ fMaxTimeRange = param.fMaxTimeRange;
+ fValueToSize = param.fValueToSize;
+ fMaxPadRangeCm = param.fMaxPadRangeCm;
+ fMaxRowRangeCm = param.fMaxRowRangeCm;
+ fDebugLevel = param.fDebugLevel;
+ fHistoRow = param.fHistoRow ;
+ fHistoPad = param.fHistoPad ;
+ fHistoTime = param.fHistoTime;
+ fHistoRowPad = param.fHistoRowPad;
return (*this);
}
//
// destructor
//
+ delete fOutput;
}
void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
{
+ //
+ // set reconstruction parameters
+ //
if (recoParam) {
fRecoParam = recoParam;
}else{
//
fInput = tree;
if (!fInput->GetBranch("Segment")){
- cerr<<"AliTPCclusterKr::findClusterKr(): no proper input tree !\n";
+ cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
fInput=0;
return;
}
}
-void AliTPCclustererKr::SetOutput(TTree * tree)
+void AliTPCclustererKr::SetOutput(TTree * /*tree*/)
{
//
- // set output
+ //dummy
//
- fOutput= tree;
- AliTPCClustersRow clrow;
- AliTPCClustersRow *pclrow=&clrow;
- clrow.SetClass("AliTPCclusterKr");
- clrow.SetArray(1); // to make Clones array
- fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
+ fOutput = new TTreeSRedirector("Krypton.root");
}
////____________________________________________________________________________
//// with new I/O
-Int_t AliTPCclustererKr::finderIO()
+Int_t AliTPCclustererKr::FinderIO()
{
// Krypton cluster finder for simulated events from MC
return 11;
}
- findClusterKrIO();
+ FindClusterKrIO();
return 0;
}
-Int_t AliTPCclustererKr::finderIO(AliRawReader* rawReader)
+Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
{
// Krypton cluster finder for the TPC raw data
//
// fParam must be defined before
- // consider noiceROC or not
-
if(rawReader)fRawData=kTRUE; //set flag to data
if (!fOutput) {
Bool_t isAltro=kFALSE;
AliTPCROC * roc = AliTPCROC::Instance();
- //AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
+ AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
//
AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors
const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors
const Int_t kNS = kNIS + kNOS;//all sectors
- Int_t zeroSup = fParam->GetZeroSup();//zero suppression parameter
+
//crate TPC view
AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
//
// Loop over sectors
//
- for(Int_t sec = 0; sec < kNS; sec++) {
- //AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sec); // noise per given sector
-
+ for(Int_t iSec = 0; iSec < kNS; iSec++) {
+ AliTPCCalROC * noiseROC;
+ if(noiseTPC==0x0){
+ noiseROC =new AliTPCCalROC(iSec);//noise=0
+ }
+ else{
+ noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector
+ }
Int_t nRows = 0; //number of rows in sector
Int_t nDDLs = 0; //number of DDLs
Int_t indexDDL = 0; //DDL index
- if (sec < kNIS) {
+ if (iSec < kNIS) {
nRows = fParam->GetNRowLow();
nDDLs = 2;
- indexDDL = sec * 2;
+ indexDDL = iSec * 2;
}else {
nRows = fParam->GetNRowUp();
nDDLs = 4;
- indexDDL = (sec-kNIS) * 4 + kNIS * 2;
+ indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
}
//
// Load the raw data for corresponding DDLs
//
rawReader->Reset();
- input.SetOldRCUFormat(fIsOldRCUFormat);
rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
if(input.Next()) {
isAltro=kTRUE;
// Allocate memory for rows in sector (pads(depends on row) x timebins)
- for(Int_t row = 0; row < nRows; row++) {
- digarr->CreateRow(sec,row);
+ for(Int_t iRow = 0; iRow < nRows; iRow++) {
+ digarr->CreateRow(iSec,iRow);
}//end loop over rows
}
+ rawReader->Reset();
rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
//
// Begin loop over altro data
//
while (input.Next()) {
-
+
//check sector consistency
- if (input.GetSector() != sec)
- AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",sec,input.GetSector()));
+ if (input.GetSector() != iSec)
+ AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector()));
+ Int_t iRow = input.GetRow();
+ Int_t iPad = input.GetPad();
+ Int_t iTimeBin = input.GetTime();
+
+ //
+ if(fDebugLevel==72){
+ fHistoRow->Fill(iRow);
+ fHistoPad->Fill(iPad);
+ fHistoTime->Fill(iTimeBin);
+ fHistoRowPad->Fill(iPad,iRow);
+ }else if(fDebugLevel>=0&&fDebugLevel<72){
+ if(iSec==fDebugLevel){
+ fHistoRow->Fill(iRow);
+ fHistoPad->Fill(iPad);
+ fHistoTime->Fill(iTimeBin);
+ fHistoRowPad->Fill(iPad,iRow);
+ }
+ }else if(fDebugLevel==73){
+ if(iSec<36){
+ fHistoRow->Fill(iRow);
+ fHistoPad->Fill(iPad);
+ fHistoTime->Fill(iTimeBin);
+ fHistoRowPad->Fill(iPad,iRow);
+ }
+ }else if(fDebugLevel==74){
+ if(iSec>=36){
+ fHistoRow->Fill(iRow);
+ fHistoPad->Fill(iPad);
+ fHistoTime->Fill(iTimeBin);
+ fHistoRowPad->Fill(iPad,iRow);
+ }
+ }
+
//check row consistency
- Short_t iRow = input.GetRow();
+ if (iRow < 0 ) continue;
if (iRow < 0 || iRow >= nRows){
AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
iRow, 0, nRows -1));
}
//check pad consistency
- Short_t iPad = input.GetPad();
- if (iPad < 0 || iPad >= (Short_t)(roc->GetNPads(sec,iRow))) {
+ if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
- iPad, 0, roc->GetNPads(sec,iRow) ));
+ iPad, 0, roc->GetNPads(iSec,iRow) ));
continue;
}
//check time consistency
- Short_t iTimeBin = input.GetTime();
if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
//cout<<iTimeBin<<endl;
continue;
//signal
Int_t signal = input.GetSignal();
- if (signal <= zeroSup) {
- //continue;
- digarr->GetRow(sec,iRow)->SetDigitFast(0,iTimeBin,iPad);
+ if (signal <= fZeroSup ||
+ iTimeBin < fFirstBin ||
+ iTimeBin > fLastBin
+ ) {
+ digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
+ continue;
+ }
+ if (!noiseROC) continue;
+ Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector
+ if (noiseOnPad > fMaxNoiseAbs){
+ digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
+ continue; // consider noisy pad as dead
}
- (//(AliSimDigits*)
- digarr->GetRow(sec,iRow))->SetDigitFast(signal,iTimeBin,iPad);
+ if(signal <= fMaxNoiseSigma * noiseOnPad){
+ digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
+ continue;
+ }
+ digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad);
}//end of loop over altro data
}//end of loop over sectors
SetDigArr(digarr);
- if(isAltro) findClusterKrIO();
+ if(isAltro) FindClusterKrIO();
delete digarr;
return 0;
}
////____________________________________________________________________________
-Int_t AliTPCclustererKr::findClusterKrIO()
+Int_t AliTPCclustererKr::FindClusterKrIO()
{
+
+ //
//fParam and fDigarr must be set to run this method
+ //
- Int_t cluster_counter=0;
- const Short_t NTotalSector=fParam->GetNSector();//number of sectors
- for(Short_t sec=0; sec<NTotalSector; sec++){
+ Int_t clusterCounter=0;
+ const Int_t nTotalSector=fParam->GetNSector();//number of sectors
+ for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
//vector of maxima for each sector
- std::vector<AliPadMax*> maxima_in_sector;
-
- const Short_t Nrows=fParam->GetNRow(sec);//number of rows in sector
- for(Short_t row=0; row<Nrows; row++){
+ //std::vector<AliPadMax*> maximaInSector;
+ TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*
+
+ //
+ // looking for the maxima on the pad
+ //
+
+ const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector
+ for(Int_t iRow=0; iRow<kNRows; ++iRow){
AliSimDigits *digrow;
if(fRawData){
- digrow = (AliSimDigits*)fDigarr->GetRow(sec,row);//real data
+ digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data
}else{
- digrow = (AliSimDigits*)fDigarr->LoadRow(sec,row);//MC
+ digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
}
if(digrow){//if pointer exist
digrow->ExpandBuffer(); //decrunch
- const Short_t npads = digrow->GetNCols(); // number of pads
- const Short_t ntime = digrow->GetNRows(); // number of timebins
- for(Short_t np=0;np<npads;np++){
+ const Int_t kNPads = digrow->GetNCols(); // number of pads
+ const Int_t kNTime = digrow->GetNRows(); // number of timebins
+ for(Int_t iPad=0;iPad<kNPads;iPad++){
- Short_t tb_max=-1;//timebin of maximum
- Short_t value_maximum=-1;//value of maximum in adc
- Short_t increase_begin=-1;//timebin when increase starts
- Short_t sum_adc=0;//sum of adc on the pad in maximum surrounding
- bool if_increase_begin=true;//flag - check if increasing start
- bool if_maximum=false;//flag - check if it could be maximum
-
- for(Short_t nt=0;nt<ntime;nt++){
- Short_t adc = digrow->GetDigitFast(nt,np);
- if(adc<3){
- if(if_maximum){
- if(nt-1-increase_begin<1){//at least 2 time bins
- tb_max=-1;
- value_maximum=-1;
- increase_begin=-1;
- sum_adc=0;
- if_increase_begin=true;
- if_maximum=false;
+ Int_t timeBinMax=-1;//timebin of maximum
+ Int_t valueMaximum=-1;//value of maximum in adc
+ Int_t increaseBegin=-1;//timebin when increase starts
+ Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding
+ bool ifIncreaseBegin=true;//flag - check if increasing started
+ bool ifMaximum=false;//flag - check if it could be maximum
+ Short_t* val = digrow->GetDigitsColumn(iPad);
+ for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
+ if (!ifMaximum) {
+ if (val[iTimeBin]==-1) break; // 0 until the end
+ for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++);
+ }
+ //
+ Short_t adc = val[iTimeBin];
+ if(adc<fMinAdc){//standard was 3
+ if(ifMaximum){
+ if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins
+ timeBinMax=-1;
+ valueMaximum=-1;
+ increaseBegin=-1;
+ sumAdc=0;
+ ifIncreaseBegin=true;
+ ifMaximum=false;
continue;
}
//insert maximum, default values and set flags
- AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
- tb_max,
- np,
- row),
- increase_begin,
- nt-1,
- sum_adc);
- maxima_in_sector.push_back(one_maximum);
+ //Double_t xCord,yCord;
+ //GetXY(iSec,iRow,iPad,xCord,yCord);
+ Double_t x[]={iRow,iPad,iTimeBin};
+ Int_t i[]={iSec};
+ AliTPCTransform trafo;
+ trafo.Transform(x,i,0,1);
+
+ AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
+ timeBinMax,
+ iPad,
+ iRow,
+ x[0],//xCord,
+ x[1],//yCord,
+ x[2]/*timeBinMax*/),
+ increaseBegin,
+ iTimeBin-1,
+ sumAdc);
+ maximaInSector->AddLast(oneMaximum);
- tb_max=-1;
- value_maximum=-1;
- increase_begin=-1;
- sum_adc=0;
- if_increase_begin=true;
- if_maximum=false;
+ timeBinMax=-1;
+ valueMaximum=-1;
+ increaseBegin=-1;
+ sumAdc=0;
+ ifIncreaseBegin=true;
+ ifMaximum=false;
}
continue;
}
- if(if_increase_begin){
- if_increase_begin=false;
- increase_begin=nt;
+ if(ifIncreaseBegin){
+ ifIncreaseBegin=false;
+ increaseBegin=iTimeBin;
}
- if(adc>value_maximum){
- tb_max=nt;
- value_maximum=adc;
- if_maximum=true;
+ if(adc>valueMaximum){
+ timeBinMax=iTimeBin;
+ valueMaximum=adc;
+ ifMaximum=true;
}
- sum_adc+=adc;
- if(nt==ntime-1 && if_maximum && ntime-1-increase_begin>2){//on the edge
+ sumAdc+=adc;
+ if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
+ //at least 3 timebins
//insert maximum, default values and set flags
- AliPadMax *one_maximum = new AliPadMax(AliTPCvtpr(value_maximum,
- tb_max,
- np,
- row),
- increase_begin,
- nt-1,
- sum_adc);
- maxima_in_sector.push_back(one_maximum);
-
- tb_max=-1;
- value_maximum=-1;
- increase_begin=-1;
- sum_adc=0;
- if_increase_begin=true;
- if_maximum=false;
+ //Double_t xCord,yCord;
+ //GetXY(iSec,iRow,iPad,xCord,yCord);
+ Double_t x[]={iRow,iPad,iTimeBin};
+ Int_t i[]={iSec};
+ AliTPCTransform trafo;
+ trafo.Transform(x,i,0,1);
+ AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,
+ timeBinMax,
+ iPad,
+ iRow,
+ x[0],//xCord,
+ x[1],//yCord,
+ x[2]/*timeBinMax*/),
+ increaseBegin,
+ iTimeBin-1,
+ sumAdc);
+ maximaInSector->AddLast(oneMaximum);
+
+ timeBinMax=-1;
+ valueMaximum=-1;
+ increaseBegin=-1;
+ sumAdc=0;
+ ifIncreaseBegin=true;
+ ifMaximum=false;
continue;
}
// cout<<"Pointer does not exist!!"<<endl;
}//end if poiner exists
}//end loop over rows
-
- Short_t max_dig=0;
- Short_t max_sum_adc=0;
- Short_t max_nt=0;
- Short_t max_np=0;
- Short_t max_row=0;
-
- for( std::vector<AliPadMax*>::iterator mp1 = maxima_in_sector.begin();
- mp1 != maxima_in_sector.end(); ++mp1 ) {
-
- AliTPCclusterKr *tmp=new AliTPCclusterKr();
-
- Short_t n_used_pads=1;
- Short_t cluster_value=0;
- cluster_value+=(*mp1)->GetSum();
-
- max_dig =(*mp1)->GetAdc() ;
- max_sum_adc =(*mp1)->GetSum() ;
- max_nt =(*mp1)->GetTime();
- max_np =(*mp1)->GetPad() ;
- max_row =(*mp1)->GetRow() ;
-
- AliSimDigits *digrow_tmp;
- if(fRawData){
- digrow_tmp = (AliSimDigits*)fDigarr->GetRow(sec,(*mp1)->GetRow());
- }else{
- digrow_tmp = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp1)->GetRow());
- }
+ MakeClusters(maximaInSector,iSec,clusterCounter);
+ //
+ maximaInSector->SetOwner(kTRUE);
+ maximaInSector->Delete();
+ delete maximaInSector;
+ }//end sector for
+ cout<<"Number of clusters in event: "<<clusterCounter<<endl;
+ return 0;
+}
- digrow_tmp->ExpandBuffer(); //decrunch
- for(Short_t tb=(*mp1)->GetBegin(); tb<((*mp1)->GetEnd())+1; tb++){
- Short_t adc_tmp = digrow_tmp->GetDigitFast(tb,(*mp1)->GetPad());
- AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp1)->GetPad(),(*mp1)->GetRow());
- tmp->fCluster.push_back(vtpr);
- }
-
- maxima_in_sector.erase(mp1);
- mp1--;
-
- for( std::vector<AliPadMax*>::iterator mp2 = maxima_in_sector.begin();
- mp2 != maxima_in_sector.end(); ++mp2 ) {
+void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){
+ //
+ // Make clusters
+ //
+
+ Int_t maxDig=0;
+ Int_t maxSumAdc=0;
+ Int_t maxTimeBin=0;
+ Int_t maxPad=0;
+ Int_t maxRow=0;
+ Double_t maxX=0;
+ Double_t maxY=0;
+ Double_t maxT=0;
+ Int_t entriesArr = maximaInSector->GetEntriesFast();
+ for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {
+
+ AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);
+ if (!mp1) continue;
+ AliTPCclusterKr clusterKr;
+
+ Int_t nUsedPads=1;
+ Int_t clusterValue=0;
+ clusterValue+=(mp1)->GetSum();
+ list<Int_t> nUsedRows;
+ nUsedRows.push_back((mp1)->GetRow());
+
+ maxDig =(mp1)->GetAdc() ;
+ maxSumAdc =(mp1)->GetSum() ;
+ maxTimeBin =(mp1)->GetTime();
+ maxPad =(mp1)->GetPad() ;
+ maxRow =(mp1)->GetRow() ;
+ maxX =(mp1)->GetX();
+ maxY =(mp1)->GetY();
+ maxT =(mp1)->GetT();
+
+ AliSimDigits *digrowTmp;
+ if(fRawData){
+ digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());
+ }else{
+ digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());
+ }
+
+ digrowTmp->ExpandBuffer(); //decrunch
+
+ for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){
+ Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());
+ AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());
+ clusterKr.AddDigitToCluster(vtpr);
+ }
+ clusterKr.SetCenter();//set centr of the cluster
+
+ for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {
+ AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);
+ if (!mp2) continue;
+ if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX())>fMaxPadRangeCm) continue;
+ if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue;
+ if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) >fMaxTimeRange) continue;
+
+ {
+ clusterValue+=(mp2)->GetSum();
- Short_t two_pad_max=4;
- Short_t two_row_max=3;
+ nUsedPads++;
+ nUsedRows.push_back((mp2)->GetRow());
- if(abs(max_row - (*mp2)->GetRow())<two_row_max &&
- abs(max_np - (*mp2)->GetPad())<two_pad_max &&
- abs(max_nt - (*mp2)->GetTime())<7){
-
- cluster_value+=(*mp2)->GetSum();
- n_used_pads++;
-
- AliSimDigits *digrow_tmp1;
- if(fRawData){
- digrow_tmp1 = (AliSimDigits*)fDigarr->GetRow(sec,(*mp2)->GetRow());
- }else{
- digrow_tmp1 = (AliSimDigits*)fDigarr->LoadRow(sec,(*mp2)->GetRow());
- }
- digrow_tmp1->ExpandBuffer(); //decrunch
-
- for(Short_t tb=(*mp2)->GetBegin(); tb<(*mp2)->GetEnd()+1; tb++){
- Short_t adc_tmp = digrow_tmp1->GetDigitFast(tb,(*mp2)->GetPad());
- AliTPCvtpr *vtpr=new AliTPCvtpr(adc_tmp,tb,(*mp2)->GetPad(),(*mp2)->GetRow());
- tmp->fCluster.push_back(vtpr);
- }
-
- //which one is bigger
- if( (*mp2)->GetAdc() > max_dig ){
- max_dig =(*mp2)->GetAdc() ;
- max_sum_adc =(*mp2)->GetSum() ;
- max_nt =(*mp2)->GetTime();
- max_np =(*mp2)->GetPad() ;
- max_row =(*mp2)->GetRow() ;
- } else if ( (*mp2)->GetAdc() == max_dig ){
- if( (*mp2)->GetSum() > max_sum_adc){
- max_dig =(*mp2)->GetAdc() ;
- max_sum_adc =(*mp2)->GetSum() ;
- max_nt =(*mp2)->GetTime();
- max_np =(*mp2)->GetPad() ;
- max_row =(*mp2)->GetRow() ;
- }
+ AliSimDigits *digrowTmp1;
+ if(fRawData){
+ digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());
+ }else{
+ digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());
+ }
+ digrowTmp1->ExpandBuffer(); //decrunch
+
+ for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){
+ Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());
+ AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());
+ clusterKr.AddDigitToCluster(vtpr);
+ }
+
+ clusterKr.SetCenter();//set center of the cluster
+
+ //which one is bigger
+ if( (mp2)->GetAdc() > maxDig ){
+ maxDig =(mp2)->GetAdc() ;
+ maxSumAdc =(mp2)->GetSum() ;
+ maxTimeBin =(mp2)->GetTime();
+ maxPad =(mp2)->GetPad() ;
+ maxRow =(mp2)->GetRow() ;
+ maxX =(mp2)->GetX() ;
+ maxY =(mp2)->GetY() ;
+ maxT =(mp2)->GetT() ;
+ } else if ( (mp2)->GetAdc() == maxDig ){
+ if( (mp2)->GetSum() > maxSumAdc){
+ maxDig =(mp2)->GetAdc() ;
+ maxSumAdc =(mp2)->GetSum() ;
+ maxTimeBin =(mp2)->GetTime();
+ maxPad =(mp2)->GetPad() ;
+ maxRow =(mp2)->GetRow() ;
+ maxX =(mp2)->GetX() ;
+ maxY =(mp2)->GetY() ;
+ maxT =(mp2)->GetT() ;
}
- maxima_in_sector.erase(mp2);
- mp2--;
}
- }//inside loop
-
- //through out ADC=6,7 on 1 pad, 2 tb and ADC=12 on 2 pads,2 tb
- if(n_used_pads==1 && cluster_value/tmp->fCluster.size()<3.6)continue;
- if(n_used_pads==2 && cluster_value/tmp->fCluster.size()<3.1)continue;
-
- tmp->SetADCcluster(cluster_value);
- tmp->SetNpads(n_used_pads);
- tmp->SetMax(AliTPCvtpr(max_dig,max_nt,max_np,max_row));
- tmp->SetSec(sec);
- tmp->SetSize();
-
- cluster_counter++;
-
- //save each cluster into file
+ delete maximaInSector->RemoveAt(it2);
+ }
+ }//inside loop
+ delete maximaInSector->RemoveAt(it1);
+ clusterKr.SetSize();
+ //through out clusters on the edge and noise
+ //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;
+ if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;
+
+ clusterKr.SetADCcluster(clusterValue);
+ clusterKr.SetNPads(nUsedPads);
+ clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));
+ clusterKr.SetSec(iSec);
+ clusterKr.SetSize();
+
+ nUsedRows.sort();
+ nUsedRows.unique();
+ clusterKr.SetNRows(nUsedRows.size());
+ clusterKr.SetCenter();
+
+ clusterCounter++;
+
+
+ //save each cluster into file
+ if (fOutput){
+ (*fOutput)<<"Kr"<<
+ "Cl.="<<&clusterKr<<
+ "\n";
+ }
+ //end of save each cluster into file adc.root
+ }//outer loop
+}
- AliTPCClustersRow *clrow= new AliTPCClustersRow();
- fRowCl=clrow;
- clrow->SetClass("AliTPCclusterKr");
- clrow->SetArray(1);
- fOutput->GetBranch("Segment")->SetAddress(&clrow);
- Int_t tmp_cluster=0;
- TClonesArray * arr = fRowCl->GetArray();
- AliTPCclusterKr * cl = new ((*arr)[tmp_cluster]) AliTPCclusterKr(*tmp);
-
- fOutput->Fill();
- delete clrow;
- //end of save each cluster into file adc.root
- }//outer loop
+////____________________________________________________________________________
- }//end sector for
- cout<<"Number of clusters in event: "<<cluster_counter<<endl;
- return 0;
-}
+void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){
+ //
+ //gives global XY coordinate of the pad
+ //
-////____________________________________________________________________________
+ Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm
+
+ Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row
+ Float_t padXSize;
+ if(sec<fParam->GetNInnerSector())padXSize=0.4;
+ else padXSize=0.6;
+ Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad
+
+ Float_t sin,cos;
+ fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector
+
+ Double_t xGlob1 = xLocal * cos + yLocal * sin;
+ Double_t yGlob1 = -xLocal * sin + yLocal * cos;
+
+
+ Double_t rot=0;
+ rot=TMath::Pi()/2.;
+
+ xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);
+ yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);
+
+ yGlob=-1*yGlob;
+ if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;
+
+
+ return;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff