-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- * *\r
- * Author: The ALICE Off-line Project. *\r
- * Contributors are mentioned in the code where appropriate. *\r
- * *\r
- * Permission to use, copy, modify and distribute this software and its *\r
- * documentation strictly for non-commercial purposes is hereby granted *\r
- * without fee, provided that the above copyright notice appears in all *\r
- * copies and that both the copyright notice and this permission notice *\r
- * appear in the supporting documentation. The authors make no claims *\r
- * about the suitability of this software for any purpose. It is *\r
- * provided "as is" without express or implied warranty. *\r
- **************************************************************************/\r
-\r
-/* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */\r
-\r
-//-----------------------------------------------------------------\r
-// Implementation of the TPC Kr cluster class\r
-//\r
-// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl\r
-//-----------------------------------------------------------------\r
-\r
-/*\r
-Instruction - how to use that:\r
-There are two macros prepared. One is for preparing clusters from MC \r
-samples: FindKrClusters.C. The output is kept in TPC.RecPoints.root.\r
-The other macro is prepared for data analysis: FindKrClustersRaw.C. \r
-The output is created for each processed file in root file named adc.root. \r
-For each data subsample the same named file is created. So be careful \r
-do not overwrite them. \r
-\r
-Additional selection criteria to select the GOLD cluster\r
-Example:\r
-// open file with clusters\r
-TFile f("Krypton.root");\r
-TTree * tree = (TTree*)f.Get("Kr")\r
-TCut cutR0("cutR0","fADCcluster/fSize<100"); // adjust it according v seetings - \r
-TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal\r
-TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.2"); // digital noise removal\r
-TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal\r
-TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks\r
-TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1;\r
-This values are typical values to be applied in selectors\r
-\r
-\r
-*\r
-**** MC ****\r
-*\r
-\r
-To run clusterizaton for MC type:\r
-.x FindKrClusters.C\r
-\r
-If you don't want to use the standard selection criteria then you \r
-have to do following:\r
-\r
-// load the standard setup\r
-AliRunLoader* rl = AliRunLoader::Open("galice.root");\r
-AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");\r
-tpcl->LoadDigits();\r
-rl->LoadgAlice();\r
-gAlice=rl->GetAliRun();\r
-TDirectory *cwd = gDirectory;\r
-AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");\r
-Int_t ver = tpc->IsVersion();\r
-rl->CdGAFile();\r
-AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");\r
-AliTPCDigitsArray *digarr=new AliTPCDigitsArray;\r
-digarr->Setup(param);\r
-cwd->cd();\r
-\r
-//loop over events\r
-Int_t nevmax=rl->GetNumberOfEvents();\r
-for(Int_t nev=0;nev<nevmax ;nev++){\r
- rl->GetEvent(nev);\r
- TTree* input_tree= tpcl->TreeD();//load tree with digits\r
- digarr->ConnectTree(input_tree);\r
- TTree *output_tree =tpcl->TreeR();//load output tree\r
-\r
- AliTPCclustererKr *clusters = new AliTPCclustererKr();\r
- clusters->SetParam(param);\r
- clusters->SetInput(input_tree);\r
- clusters->SetOutput(output_tree);\r
- clusters->SetDigArr(digarr);\r
- \r
-//If you want to change the cluster finder parameters for MC there are \r
-//several of them:\r
-\r
-//1. signal threshold (everything below the given number is treated as 0)\r
- clusters->SetMinAdc(3);\r
-\r
-//2. number of neighbouring timebins to be considered\r
- clusters->SetMinTimeBins(2);\r
-\r
-//3. distance of the cluster center to the center of a pad in pad-padrow plane \r
-//(in cm). Remenber that this is still quantified by pad size.\r
- clusters->SetMaxPadRangeCm(2.5);\r
-\r
-//4. distance of the cluster center to the center of a padrow in pad-padrow \r
-//plane (in cm). Remenber that this is still quantified by pad size.\r
- clusters->SetMaxRowRangeCm(3.5);\r
-\r
-//5. distance of the cluster center to the max time bin on a pad (in tackts)\r
-//ie. fabs(centerT - time)<7\r
- clusters->SetMaxTimeRange(7);\r
-\r
-//6. cut reduce peak at 0. There are noises which appear mostly as two \r
-//timebins on one pad.\r
- clusters->SetValueToSize(3.5);\r
-\r
-\r
- clusters->finderIO();\r
- tpcl->WriteRecPoints("OVERWRITE");\r
-}\r
-delete rl;//cleans everything\r
-\r
-*\r
-********* DATA *********\r
-*\r
-\r
-To run clusterizaton for DATA for file named raw_data.root type:\r
-.x FindKrClustersRaw.C("raw_data.root")\r
-\r
-If you want to change some criteria do the following:\r
-\r
-//\r
-// remove Altro warnings\r
-//\r
-AliLog::SetClassDebugLevel("AliRawReaderDate",-5);\r
-AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5);\r
-AliLog::SetModuleDebugLevel("RAW",-5);\r
-\r
-//\r
-// Get database with noises\r
-//\r
-// char *ocdbpath = gSystem->Getenv("OCDB_PATH");\r
-char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";\r
-if (ocdbpath==0){\r
-ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";\r
-}\r
-AliCDBManager * man = AliCDBManager::Instance();\r
-man->SetDefaultStorage(ocdbpath);\r
-man->SetRun(0);\r
-AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();\r
-AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();\r
-\r
-//define tree\r
-TFile *hfile=new TFile("adc.root","RECREATE","ADC file");\r
-// Create a ROOT Tree\r
-TTree *mytree = new TTree("Kr","Krypton cluster tree");\r
-\r
-//define infput file\r
-const char *fileName="data.root";\r
-AliRawReader *reader = new AliRawReaderRoot(fileName);\r
-//AliRawReader *reader = new AliRawReaderDate(fileName);\r
-reader->Reset();\r
-AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader);\r
-stream->SelectRawData("TPC");\r
-\r
-//one general output\r
-AliTPCclustererKr *clusters = new AliTPCclustererKr();\r
-clusters->SetOutput(mytree);\r
-clusters->SetRecoParam(0);//standard reco parameters\r
-AliTPCParamSR *param=new AliTPCParamSR();\r
-clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.)\r
-\r
-//set cluster finder parameters (from data):\r
-//1. zero suppression parameter\r
- clusters->SetZeroSup(param->GetZeroSup());\r
-\r
-//2. first bin\r
- clusters->SetFirstBin(60);\r
-\r
-//3. last bin\r
- clusters->SetLastBin(950);\r
-\r
-//4. maximal noise\r
- clusters->SetMaxNoiseAbs(2);\r
-\r
-//5. maximal amount of sigma of noise\r
- clusters->SetMaxNoiseSigma(3);\r
-\r
-//The remaining parameters are the same paramters as for MC (see MC section \r
-//points 1-6)\r
- clusters->SetMinAdc(3);\r
- clusters->SetMinTimeBins(2);\r
- clusters->SetMaxPadRangeCm(2.5);\r
- clusters->SetMaxRowRangeCm(3.5);\r
- clusters->SetMaxTimeRange(7);\r
- clusters->SetValueToSize(3.5);\r
-\r
-while (reader->NextEvent()) {\r
- clusters->FinderIO(reader);\r
-}\r
-\r
-hfile->Write();\r
-hfile->Close();\r
-delete stream;\r
-\r
-\r
-*/\r
-\r
-#include "AliTPCclustererKr.h"\r
-#include "AliTPCclusterKr.h"\r
-//#include <vector>\r
-#include <list>\r
-#include "TObject.h"\r
-#include "AliPadMax.h"\r
-#include "AliSimDigits.h"\r
-#include "AliTPCv4.h"\r
-#include "AliTPCParam.h"\r
-#include "AliTPCDigitsArray.h"\r
-#include "AliTPCvtpr.h"\r
-#include "AliTPCClustersRow.h"\r
-#include "TTree.h"\r
-#include "TH1F.h"\r
-#include "TH2F.h"\r
-#include "TTreeStream.h"\r
-\r
-#include "AliTPCTransform.h"\r
-\r
-//used in raw data finder\r
-#include "AliTPCROC.h"\r
-#include "AliTPCCalPad.h"\r
-#include "AliTPCAltroMapping.h"\r
-#include "AliTPCcalibDB.h"\r
-#include "AliTPCRawStreamV3.h"\r
-#include "AliTPCRecoParam.h"\r
-#include "AliTPCReconstructor.h"\r
-#include "AliRawReader.h"\r
-#include "AliTPCCalROC.h"\r
-#include "AliRawEventHeaderBase.h"\r
-\r
-using std::cerr;\r
-using std::cout;\r
-using std::endl;\r
-using std::list;\r
-ClassImp(AliTPCclustererKr)\r
-\r
-\r
-AliTPCclustererKr::AliTPCclustererKr()\r
- :TObject(),\r
- fRawData(kFALSE),\r
- fInput(0),\r
- fOutput(0),\r
- fParam(0),\r
- fDigarr(0),\r
- fRecoParam(0),\r
- fZeroSup(2),\r
- fFirstBin(60),\r
- fLastBin(950),\r
- fMaxNoiseAbs(2),\r
- fMaxNoiseSigma(3),\r
- fMinAdc(3),\r
- fMinTimeBins(2),\r
-// fMaxPadRange(4),\r
-// fMaxRowRange(3),\r
- fMaxTimeRange(7),\r
- fValueToSize(3.5),\r
- fMaxPadRangeCm(2.5),\r
- fMaxRowRangeCm(3.5),\r
- fIsolCut(3),\r
- fDebugLevel(-1),\r
- fHistoRow(0),\r
- fHistoPad(0),\r
- fHistoTime(0),\r
- fHistoRowPad(0),\r
- fTimeStamp(0),\r
- fRun(0)\r
-{\r
-//\r
-// default constructor\r
-//\r
-}\r
-\r
-AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m)\r
- :TObject(),\r
- fRawData(kFALSE),\r
- fInput(0),\r
- fOutput(0),\r
- fParam(0),\r
- fDigarr(0),\r
- fRecoParam(0),\r
- fZeroSup(2),\r
- fFirstBin(60),\r
- fLastBin(950),\r
- fMaxNoiseAbs(2),\r
- fMaxNoiseSigma(3),\r
- fMinAdc(3),\r
- fMinTimeBins(2),\r
-// fMaxPadRange(4),\r
-// fMaxRowRange(3),\r
- fMaxTimeRange(7),\r
- fValueToSize(3.5),\r
- fMaxPadRangeCm(2.5),\r
- fMaxRowRangeCm(3.5),\r
- fIsolCut(3),\r
- fDebugLevel(-1),\r
- fHistoRow(0),\r
- fHistoPad(0),\r
- fHistoTime(0),\r
- fHistoRowPad(0),\r
- fTimeStamp(0),\r
- fRun(0)\r
-{\r
-//\r
-// copy constructor\r
-//\r
- fParam = param.fParam;\r
- fRecoParam = param.fRecoParam;\r
- fRawData = param.fRawData;\r
- fInput = param.fInput ;\r
- fOutput = param.fOutput;\r
- fDigarr = param.fDigarr;\r
- fZeroSup = param.fZeroSup ;\r
- fFirstBin = param.fFirstBin ;\r
- fLastBin = param.fLastBin ;\r
- fMaxNoiseAbs = param.fMaxNoiseAbs ;\r
- fMaxNoiseSigma = param.fMaxNoiseSigma ;\r
- fMinAdc = param.fMinAdc;\r
- fMinTimeBins = param.fMinTimeBins;\r
-// fMaxPadRange = param.fMaxPadRange ;\r
-// fMaxRowRange = param.fMaxRowRange ;\r
- fMaxTimeRange = param.fMaxTimeRange;\r
- fValueToSize = param.fValueToSize;\r
- fMaxPadRangeCm = param.fMaxPadRangeCm;\r
- fMaxRowRangeCm = param.fMaxRowRangeCm;\r
- fIsolCut = param.fIsolCut;\r
- fDebugLevel = param.fDebugLevel;\r
- fHistoRow = param.fHistoRow ;\r
- fHistoPad = param.fHistoPad ;\r
- fHistoTime = param.fHistoTime;\r
- fHistoRowPad = param.fHistoRowPad;\r
- fTimeStamp = param.fTimeStamp;\r
- fRun = param.fRun;\r
-\r
-} \r
-\r
-AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)\r
-{\r
- //\r
- // assignment operator\r
- //\r
- if (this == ¶m) return (*this);\r
- \r
- fParam = param.fParam;\r
- fRecoParam = param.fRecoParam;\r
- fRawData = param.fRawData;\r
- fInput = param.fInput ;\r
- fOutput = param.fOutput;\r
- fDigarr = param.fDigarr;\r
- fZeroSup = param.fZeroSup ;\r
- fFirstBin = param.fFirstBin ;\r
- fLastBin = param.fLastBin ;\r
- fMaxNoiseAbs = param.fMaxNoiseAbs ;\r
- fMaxNoiseSigma = param.fMaxNoiseSigma ;\r
- fMinAdc = param.fMinAdc;\r
- fMinTimeBins = param.fMinTimeBins;\r
-// fMaxPadRange = param.fMaxPadRange ;\r
-// fMaxRowRange = param.fMaxRowRange ;\r
- fMaxTimeRange = param.fMaxTimeRange;\r
- fValueToSize = param.fValueToSize;\r
- fMaxPadRangeCm = param.fMaxPadRangeCm;\r
- fMaxRowRangeCm = param.fMaxRowRangeCm;\r
- fIsolCut = param.fIsolCut;\r
- fDebugLevel = param.fDebugLevel;\r
- fHistoRow = param.fHistoRow ;\r
- fHistoPad = param.fHistoPad ;\r
- fHistoTime = param.fHistoTime;\r
- fHistoRowPad = param.fHistoRowPad;\r
- fTimeStamp = param.fTimeStamp;\r
- fRun = param.fRun;\r
- return (*this);\r
-}\r
-\r
-AliTPCclustererKr::~AliTPCclustererKr()\r
-{\r
- //\r
- // destructor\r
- //\r
- delete fOutput;\r
-}\r
-\r
-void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)\r
-{\r
- //\r
- // set reconstruction parameters\r
- //\r
- if (recoParam) {\r
- fRecoParam = recoParam;\r
- }else{\r
- //set default parameters if not specified\r
- fRecoParam = AliTPCReconstructor::GetRecoParam();\r
- if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();\r
- }\r
- return;\r
-}\r
-\r
-\r
-////____________________________________________________________________________\r
-//// I/O\r
-void AliTPCclustererKr::SetInput(TTree * tree)\r
-{\r
- //\r
- // set input tree with digits\r
- //\r
- fInput = tree; \r
- if (!fInput->GetBranch("Segment")){\r
- cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";\r
- fInput=0;\r
- return;\r
- }\r
-}\r
-\r
-void AliTPCclustererKr::SetOutput(TTree * /*tree*/) \r
-{\r
- //\r
- //dummy\r
- //\r
- fOutput = new TTreeSRedirector("Krypton.root");\r
-}\r
-\r
-////____________________________________________________________________________\r
-//// with new I/O\r
-Int_t AliTPCclustererKr::FinderIO()\r
-{\r
- // Krypton cluster finder for simulated events from MC\r
-\r
- if (!fInput) { \r
- Error("Digits2Clusters", "input tree not initialised");\r
- return 10;\r
- }\r
- \r
- if (!fOutput) {\r
- Error("Digits2Clusters", "output tree not initialised");\r
- return 11;\r
- }\r
-\r
- FindClusterKrIO();\r
- return 0;\r
-}\r
-\r
-\r
-\r
-Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)\r
-{\r
- // Krypton cluster finder for the TPC raw data\r
- // this method is unsing AliAltroRawStreamV3\r
- // fParam must be defined before\r
- if (!rawReader) return 1;\r
- //\r
- fRawData=kTRUE; //set flag to data\r
- \r
- if (!fOutput) {\r
- Error("Digits2Clusters", "output tree not initialised");\r
- return 11;\r
- }\r
- \r
- fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param\r
- // used later for memory allocation\r
-\r
- AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();\r
- if (eventHeader){\r
- fTimeStamp = eventHeader->Get("Timestamp");\r
- fRun = rawReader->GetRunNumber();\r
- }\r
-\r
-\r
- Bool_t isAltro=kFALSE;\r
- \r
- AliTPCROC * roc = AliTPCROC::Instance();\r
- AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();\r
- AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();\r
- //\r
- AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);\r
- \r
- const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors\r
- const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors\r
- const Int_t kNS = kNIS + kNOS;//all sectors\r
- \r
- \r
- //crate TPC view\r
- AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim\r
- digarr->Setup(fParam);//as usually parameters\r
- \r
- for(Int_t iSec = 0; iSec < kNS; iSec++) {\r
- AliTPCCalROC * noiseROC;\r
- AliTPCCalROC noiseDummy(iSec);\r
- if(noiseTPC==0x0){\r
- noiseROC = &noiseDummy;//noise=0\r
- }else{\r
- noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector\r
- }\r
- Int_t nRows = 0; //number of rows in sector\r
- Int_t nDDLs = 0; //number of DDLs\r
- Int_t indexDDL = 0; //DDL index\r
- if (iSec < kNIS) {\r
- nRows = fParam->GetNRowLow();\r
- nDDLs = 2;\r
- indexDDL = iSec * 2;\r
- }else {\r
- nRows = fParam->GetNRowUp();\r
- nDDLs = 4;\r
- indexDDL = (iSec-kNIS) * 4 + kNIS * 2;\r
- }\r
- \r
- //\r
- // Load the raw data for corresponding DDLs\r
- //\r
- rawReader->Reset();\r
- rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);\r
- \r
- \r
- while (input.NextDDL()){\r
- // Allocate memory for rows in sector (pads(depends on row) x timebins)\r
- if (!digarr->GetRow(iSec,0)){\r
- for(Int_t iRow = 0; iRow < nRows; iRow++) {\r
- digarr->CreateRow(iSec,iRow);\r
- }//end loop over rows\r
- }\r
- //loop over pads\r
- while ( input.NextChannel() ) {\r
- Int_t iRow = input.GetRow();\r
- Int_t iPad = input.GetPad();\r
- //check row consistency\r
- if (iRow < 0 ) continue;\r
- if (iRow < 0 || iRow >= nRows){\r
- AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",\r
- iRow, 0, nRows -1));\r
- continue;\r
- }\r
- \r
- //check pad consistency\r
- if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {\r
- AliError(Form("Pad index (%d) outside the range (%d -> %d) !",\r
- iPad, 0, roc->GetNPads(iSec,iRow) ));\r
- continue;\r
- }\r
- \r
- //loop over bunches\r
- while ( input.NextBunch() ){\r
- Int_t startTbin = (Int_t)input.GetStartTimeBin();\r
- Int_t bunchlength = (Int_t)input.GetBunchLength();\r
- const UShort_t *sig = input.GetSignals();\r
- isAltro=kTRUE;\r
- for (Int_t iTime = 0; iTime<bunchlength; iTime++){\r
- Int_t iTimeBin=startTbin-iTime;\r
- //\r
- if(fDebugLevel==72){\r
- fHistoRow->Fill(iRow);\r
- fHistoPad->Fill(iPad);\r
- fHistoTime->Fill(iTimeBin);\r
- fHistoRowPad->Fill(iPad,iRow);\r
- }else if(fDebugLevel>=0&&fDebugLevel<72){\r
- if(iSec==fDebugLevel){\r
- fHistoRow->Fill(iRow);\r
- fHistoPad->Fill(iPad);\r
- fHistoTime->Fill(iTimeBin);\r
- fHistoRowPad->Fill(iPad,iRow);\r
- }\r
- }else if(fDebugLevel==73){\r
- if(iSec<36){\r
- fHistoRow->Fill(iRow);\r
- fHistoPad->Fill(iPad);\r
- fHistoTime->Fill(iTimeBin);\r
- fHistoRowPad->Fill(iPad,iRow);\r
- }\r
- }else if(fDebugLevel==74){\r
- if(iSec>=36){\r
- fHistoRow->Fill(iRow);\r
- fHistoPad->Fill(iPad);\r
- fHistoTime->Fill(iTimeBin);\r
- fHistoRowPad->Fill(iPad,iRow);\r
- }\r
- }\r
- \r
- //check time consistency\r
- if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){\r
- //cout<<iTimeBin<<endl;\r
- continue;\r
- AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",\r
- iTimeBin, 0, fRecoParam->GetLastBin() -1));\r
- }\r
- //signal\r
- Float_t signal=(Float_t)sig[iTime];\r
- if (signal <= fZeroSup ||\r
- iTimeBin < fFirstBin ||\r
- iTimeBin > fLastBin\r
- ) {\r
- digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);\r
- continue;\r
- }\r
- if (!noiseROC) continue;\r
- Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector\r
- if (noiseOnPad > fMaxNoiseAbs){\r
- digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);\r
- continue; // consider noisy pad as dead\r
- }\r
- if(signal <= fMaxNoiseSigma * noiseOnPad){\r
- digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);\r
- continue;\r
- }\r
- digarr->GetRow(iSec,iRow)->SetDigitFast(TMath::Nint(signal),iTimeBin,iPad);\r
- }// end loop signals in bunch\r
- }// end loop bunches\r
- } // end loop pads\r
- }// end ddl loop\r
- }// end sector loop\r
- SetDigArr(digarr);\r
- if(isAltro) FindClusterKrIO();\r
- delete digarr;\r
- \r
- return 0;\r
-}\r
-\r
-void AliTPCclustererKr::CleanSector(Int_t sector){\r
- //\r
- // clean isolated digits\r
- // \r
- const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector\r
- for(Int_t iRow=0; iRow<kNRows; ++iRow){\r
- AliSimDigits *digrow;\r
- if(fRawData){\r
- digrow = (AliSimDigits*)fDigarr->GetRow(sector,iRow);//real data\r
- }else{\r
- digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC\r
- }\r
- if(!digrow) continue;\r
- digrow->ExpandBuffer(); //decrunch\r
- const Int_t kNPads = digrow->GetNCols(); // number of pads\r
- const Int_t kNTime = digrow->GetNRows(); // number of timebins\r
- for(Int_t iPad=1;iPad<kNPads-1;iPad++){\r
- Short_t* val = digrow->GetDigitsColumn(iPad);\r
-\r
- for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){\r
- if (val[iTimeBin]<=0) continue;\r
- if (val[iTimeBin-1]+val[iTimeBin+1]<fIsolCut) {val[iTimeBin]=0; continue;}\r
- if (val[iTimeBin-kNTime]+val[iTimeBin+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}\r
- //\r
- if (val[iTimeBin-1-kNTime]+val[iTimeBin+1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}\r
- if (val[iTimeBin+1-kNTime]+val[iTimeBin-1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}\r
-\r
- }\r
- }\r
- }\r
-}\r
-\r
-\r
-////____________________________________________________________________________\r
-Int_t AliTPCclustererKr::FindClusterKrIO()\r
-{\r
-\r
- //\r
- //fParam and fDigarr must be set to run this method\r
- //\r
-\r
- Int_t clusterCounter=0;\r
- const Int_t nTotalSector=fParam->GetNSector();//number of sectors\r
- for(Int_t iSec=0; iSec<nTotalSector; ++iSec){\r
- CleanSector(iSec);\r
-\r
- //vector of maxima for each sector\r
- //std::vector<AliPadMax*> maximaInSector;\r
- TObjArray *maximaInSector=new TObjArray();//to store AliPadMax*\r
-\r
- //\r
- // looking for the maxima on the pad\r
- //\r
-\r
- const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector\r
- for(Int_t iRow=0; iRow<kNRows; ++iRow){\r
- AliSimDigits *digrow;\r
- if(fRawData){\r
- digrow = (AliSimDigits*)fDigarr->GetRow(iSec,iRow);//real data\r
- }else{\r
- digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC\r
- }\r
- if(digrow){//if pointer exist\r
- digrow->ExpandBuffer(); //decrunch\r
- const Int_t kNPads = digrow->GetNCols(); // number of pads\r
- const Int_t kNTime = digrow->GetNRows(); // number of timebins\r
- for(Int_t iPad=0;iPad<kNPads;iPad++){\r
- \r
- Int_t timeBinMax=-1;//timebin of maximum \r
- Int_t valueMaximum=-1;//value of maximum in adc\r
- Int_t increaseBegin=-1;//timebin when increase starts\r
- Int_t sumAdc=0;//sum of adc on the pad in maximum surrounding\r
- bool ifIncreaseBegin=true;//flag - check if increasing started\r
- bool ifMaximum=false;//flag - check if it could be maximum\r
- Short_t* val = digrow->GetDigitsColumn(iPad);\r
- for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){\r
- if (!ifMaximum) {\r
- if (val[iTimeBin]==-1) break; // 0 until the end\r
- for( ; iTimeBin<kNTime-2&&val[iTimeBin]<fMinAdc ;iTimeBin++) {}\r
- }\r
- //\r
- Short_t adc = val[iTimeBin];\r
-\r
- if(adc<fMinAdc){//standard was 3 for fMinAdc\r
- if(ifMaximum){\r
- if(iTimeBin-increaseBegin<fMinTimeBins){//at least 2 time bins\r
- timeBinMax=-1;\r
- valueMaximum=-1;\r
- increaseBegin=-1;\r
- sumAdc=0;\r
- ifIncreaseBegin=true;\r
- ifMaximum=false;\r
- continue;\r
- }\r
- //insert maximum, default values and set flags\r
- //Double_t xCord,yCord;\r
- //GetXY(iSec,iRow,iPad,xCord,yCord);\r
- Double_t x[]={iRow,iPad,iTimeBin};\r
- Int_t i[]={iSec};\r
- AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;\r
-\r
- transform->Transform(x,i,0,1);\r
- \r
- AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,\r
- timeBinMax,\r
- iPad,\r
- iRow,\r
- x[0],//xCord,\r
- x[1],//yCord,\r
- x[2]/*timeBinMax*/),\r
- increaseBegin,\r
- iTimeBin-1,\r
- sumAdc);\r
- maximaInSector->AddLast(oneMaximum);\r
- \r
- timeBinMax=-1;\r
- valueMaximum=-1;\r
- increaseBegin=-1;\r
- sumAdc=0;\r
- ifIncreaseBegin=true;\r
- ifMaximum=false;\r
- }\r
- continue;\r
- }\r
-\r
-\r
-\r
-\r
-\r
-\r
- if(ifIncreaseBegin){\r
- ifIncreaseBegin=false;\r
- increaseBegin=iTimeBin;\r
- }\r
- \r
- if(adc>valueMaximum){\r
- timeBinMax=iTimeBin;\r
- valueMaximum=adc;\r
- ifMaximum=true;\r
- }\r
- sumAdc+=adc;\r
- if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge\r
- //at least 3 timebins\r
- //insert maximum, default values and set flags\r
- //Double_t xCord,yCord;\r
- //GetXY(iSec,iRow,iPad,xCord,yCord);\r
- Double_t x[]={iRow,iPad,iTimeBin};\r
- Int_t i[]={iSec};\r
- //AliTPCTransform trafo;\r
- //trafo.Transform(x,i,0,1);\r
-\r
- AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;\r
-\r
- transform->Transform(x,i,0,1);\r
-\r
- AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum,\r
- timeBinMax,\r
- iPad,\r
- iRow,\r
- x[0],//xCord,\r
- x[1],//yCord,\r
- x[2]/*timeBinMax*/),\r
- increaseBegin,\r
- iTimeBin-1,\r
- sumAdc);\r
- maximaInSector->AddLast(oneMaximum);\r
- \r
- timeBinMax=-1;\r
- valueMaximum=-1;\r
- increaseBegin=-1;\r
- sumAdc=0;\r
- ifIncreaseBegin=true;\r
- ifMaximum=false;\r
- continue;\r
- }\r
- \r
- }//end loop over timebins\r
- }//end loop over pads\r
-// }else{\r
-// cout<<"Pointer does not exist!!"<<endl;\r
- }//end if poiner exists\r
- }//end loop over rows\r
-\r
- MakeClusters(maximaInSector,iSec,clusterCounter);\r
- //\r
- maximaInSector->SetOwner(kTRUE);\r
- maximaInSector->Delete();\r
- delete maximaInSector;\r
- }//end sector for\r
- cout<<"Number of clusters in event: "<<clusterCounter<<endl;\r
- return 0;\r
-}\r
-\r
-void AliTPCclustererKr::MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter){\r
- //\r
- // Make clusters\r
- //\r
-\r
- Int_t maxDig=0;\r
- Int_t maxSumAdc=0;\r
- Int_t maxTimeBin=0;\r
- Int_t maxPad=0;\r
- Int_t maxRow=0;\r
- Double_t maxX=0;\r
- Double_t maxY=0;\r
- Double_t maxT=0;\r
- Int_t entriesArr = maximaInSector->GetEntriesFast();\r
- for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) {\r
- \r
- AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1);\r
- if (!mp1) continue;\r
- AliTPCclusterKr clusterKr;\r
- \r
- Int_t nUsedPads=1;\r
- Int_t clusterValue=0;\r
- clusterValue+=(mp1)->GetSum();\r
- list<Int_t> nUsedRows;\r
- nUsedRows.push_back((mp1)->GetRow());\r
- \r
- maxDig =(mp1)->GetAdc() ;\r
- maxSumAdc =(mp1)->GetSum() ;\r
- maxTimeBin =(mp1)->GetTime();\r
- maxPad =(mp1)->GetPad() ;\r
- maxRow =(mp1)->GetRow() ;\r
- maxX =(mp1)->GetX();\r
- maxY =(mp1)->GetY();\r
- maxT =(mp1)->GetT();\r
- \r
- AliSimDigits *digrowTmp;\r
- if(fRawData){\r
- digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow());\r
- }else{\r
- digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow());\r
- }\r
- \r
- digrowTmp->ExpandBuffer(); //decrunch\r
- \r
- for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){\r
- Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad());\r
- AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT());\r
- clusterKr.AddDigitToCluster(vtpr);\r
- }\r
- clusterKr.SetCenter();//set centr of the cluster\r
- \r
- for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) {\r
- AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2);\r
- if (!mp2) continue;\r
- if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX()) > fMaxPadRangeCm) continue;\r
- if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue; \r
- if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) > fMaxTimeRange) continue;\r
-\r
- {\r
- clusterValue+=(mp2)->GetSum();\r
- \r
- nUsedPads++;\r
- nUsedRows.push_back((mp2)->GetRow());\r
- \r
- AliSimDigits *digrowTmp1;\r
- if(fRawData){\r
- digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow());\r
- }else{\r
- digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow());\r
- }\r
- digrowTmp1->ExpandBuffer(); //decrunch\r
- \r
- for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){\r
- Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad());\r
- AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT());\r
- clusterKr.AddDigitToCluster(vtpr);\r
- }\r
- \r
- clusterKr.SetCenter();//set center of the cluster\r
- \r
- //which one is bigger\r
- if( (mp2)->GetAdc() > maxDig ){\r
- maxDig =(mp2)->GetAdc() ;\r
- maxSumAdc =(mp2)->GetSum() ;\r
- maxTimeBin =(mp2)->GetTime();\r
- maxPad =(mp2)->GetPad() ;\r
- maxRow =(mp2)->GetRow() ;\r
- maxX =(mp2)->GetX() ;\r
- maxY =(mp2)->GetY() ;\r
- maxT =(mp2)->GetT() ;\r
- } else if ( (mp2)->GetAdc() == maxDig ){\r
- if( (mp2)->GetSum() > maxSumAdc){\r
- maxDig =(mp2)->GetAdc() ;\r
- maxSumAdc =(mp2)->GetSum() ;\r
- maxTimeBin =(mp2)->GetTime();\r
- maxPad =(mp2)->GetPad() ;\r
- maxRow =(mp2)->GetRow() ;\r
- maxX =(mp2)->GetX() ;\r
- maxY =(mp2)->GetY() ;\r
- maxT =(mp2)->GetT() ;\r
- }\r
- }\r
- delete maximaInSector->RemoveAt(it2);\r
- }\r
- }//inside loop\r
- delete maximaInSector->RemoveAt(it1); \r
- clusterKr.SetSize();\r
- //through out clusters on the edge and noise\r
- //if(clusterValue/clusterKr.fCluster.size()<fValueToSize)continue;\r
- if(clusterValue/(clusterKr.GetSize())<fValueToSize)continue;\r
- \r
- clusterKr.SetADCcluster(clusterValue);\r
- clusterKr.SetNPads(nUsedPads);\r
- clusterKr.SetMax(AliTPCvtpr(maxDig,maxTimeBin,maxPad,maxRow,maxX,maxY,maxT));\r
- clusterKr.SetSec(iSec);\r
- clusterKr.SetSize();\r
- \r
- nUsedRows.sort();\r
- nUsedRows.unique();\r
- clusterKr.SetNRows(nUsedRows.size());\r
- clusterKr.SetCenter();\r
- \r
- clusterKr.SetRMS();//Set pad,row,timebin RMS\r
- clusterKr.Set1D();//Set size in pads and timebins\r
-\r
- clusterKr.SetTimeStamp(fTimeStamp);\r
- clusterKr.SetRun(fRun);\r
-\r
- clusterCounter++;\r
- \r
- \r
- //save each cluster into file\r
- if (fOutput){\r
- (*fOutput)<<"Kr"<<\r
- "Cl.="<<&clusterKr<<\r
- "\n";\r
- }\r
- //end of save each cluster into file adc.root\r
- }//outer loop\r
-}\r
-\r
-\r
-\r
-////____________________________________________________________________________\r
-\r
-\r
-void AliTPCclustererKr::GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob){\r
- //\r
- //gives global XY coordinate of the pad\r
- //\r
-\r
- Double_t yLocal = fParam->GetPadRowRadii(sec,row);//radius of row in sector in cm\r
-\r
- Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row\r
- Float_t padXSize;\r
- if(sec<fParam->GetNInnerSector())padXSize=0.4;\r
- else padXSize=0.6;\r
- Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad\r
-\r
- Float_t sin,cos;\r
- fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector\r
-\r
- Double_t xGlob1 = xLocal * cos + yLocal * sin;\r
- Double_t yGlob1 = -xLocal * sin + yLocal * cos;\r
-\r
-\r
- Double_t rot=0;\r
- rot=TMath::Pi()/2.;\r
-\r
- xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot);\r
- yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot);\r
-\r
- yGlob=-1*yGlob;\r
- if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob;\r
-\r
-\r
- return;\r
-}\r
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */
+
+//-----------------------------------------------------------------
+// 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("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 <list>
+#include "TObject.h"
+#include "AliPadMax.h"
+#include "AliSimDigits.h"
+#include "AliTPCv4.h"
+#include "AliTPCParam.h"
+#include "AliTPCDigitsArray.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"
+#include "AliTPCCalPad.h"
+#include "AliTPCAltroMapping.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCRawStreamV3.h"
+#include "AliTPCRecoParam.h"
+#include "AliTPCReconstructor.h"
+#include "AliRawReader.h"
+#include "AliTPCCalROC.h"
+#include "AliRawEventHeaderBase.h"
+
+using std::cerr;
+using std::cout;
+using std::endl;
+using std::list;
+ClassImp(AliTPCclustererKr)
+
+
+AliTPCclustererKr::AliTPCclustererKr()
+ :TObject(),
+ fRawData(kFALSE),
+ fInput(0),
+ fOutput(0),
+ fParam(0),
+ fDigarr(0),
+ fRecoParam(0),
+ 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),
+ fIsolCut(3),
+ fDebugLevel(-1),
+ fHistoRow(0),
+ fHistoPad(0),
+ fHistoTime(0),
+ fHistoRowPad(0),
+ fTimeStamp(0),
+ fRun(0)
+{
+//
+// default constructor
+//
+}
+
+AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m)
+ :TObject(),
+ fRawData(kFALSE),
+ fInput(0),
+ fOutput(0),
+ fParam(0),
+ fDigarr(0),
+ fRecoParam(0),
+ 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),
+ fIsolCut(3),
+ fDebugLevel(-1),
+ fHistoRow(0),
+ fHistoPad(0),
+ fHistoTime(0),
+ fHistoRowPad(0),
+ fTimeStamp(0),
+ fRun(0)
+{
+//
+// copy constructor
+//
+ fParam = param.fParam;
+ 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;
+ fIsolCut = param.fIsolCut;
+ fDebugLevel = param.fDebugLevel;
+ fHistoRow = param.fHistoRow ;
+ fHistoPad = param.fHistoPad ;
+ fHistoTime = param.fHistoTime;
+ fHistoRowPad = param.fHistoRowPad;
+ fTimeStamp = param.fTimeStamp;
+ fRun = param.fRun;
+
+}
+
+AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param)
+{
+ //
+ // assignment operator
+ //
+ if (this == ¶m) return (*this);
+
+ fParam = param.fParam;
+ 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;
+ fIsolCut = param.fIsolCut;
+ fDebugLevel = param.fDebugLevel;
+ fHistoRow = param.fHistoRow ;
+ fHistoPad = param.fHistoPad ;
+ fHistoTime = param.fHistoTime;
+ fHistoRowPad = param.fHistoRowPad;
+ fTimeStamp = param.fTimeStamp;
+ fRun = param.fRun;
+ return (*this);
+}
+
+AliTPCclustererKr::~AliTPCclustererKr()
+{
+ //
+ // destructor
+ //
+ delete fOutput;
+}
+
+void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam)
+{
+ //
+ // set reconstruction parameters
+ //
+ if (recoParam) {
+ fRecoParam = recoParam;
+ }else{
+ //set default parameters if not specified
+ fRecoParam = AliTPCReconstructor::GetRecoParam();
+ if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam();
+ }
+ return;
+}
+
+
+////____________________________________________________________________________
+//// I/O
+void AliTPCclustererKr::SetInput(TTree * tree)
+{
+ //
+ // set input tree with digits
+ //
+ fInput = tree;
+ if (!fInput->GetBranch("Segment")){
+ cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n";
+ fInput=0;
+ return;
+ }
+}
+
+void AliTPCclustererKr::SetOutput(TTree * /*tree*/)
+{
+ //
+ //dummy
+ //
+ fOutput = new TTreeSRedirector("Krypton.root");
+}
+
+////____________________________________________________________________________
+//// with new I/O
+Int_t AliTPCclustererKr::FinderIO()
+{
+ // Krypton cluster finder for simulated events from MC
+
+ if (!fInput) {
+ Error("Digits2Clusters", "input tree not initialised");
+ return 10;
+ }
+
+ if (!fOutput) {
+ Error("Digits2Clusters", "output tree not initialised");
+ return 11;
+ }
+
+ FindClusterKrIO();
+ return 0;
+}
+
+
+
+Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader)
+{
+ // Krypton cluster finder for the TPC raw data
+ // this method is unsing AliAltroRawStreamV3
+ // fParam must be defined before
+ if (!rawReader) return 1;
+ //
+ fRawData=kTRUE; //set flag to data
+
+ if (!fOutput) {
+ Error("Digits2Clusters", "output tree not initialised");
+ return 11;
+ }
+
+ fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param
+ // used later for memory allocation
+
+ AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
+ if (eventHeader){
+ fTimeStamp = eventHeader->Get("Timestamp");
+ fRun = rawReader->GetRunNumber();
+ }
+
+
+ Bool_t isAltro=kFALSE;
+
+ AliTPCROC * roc = AliTPCROC::Instance();
+ AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
+ AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
+ //
+ AliTPCRawStreamV3 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
+
+
+ //crate TPC view
+ AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim
+ digarr->Setup(fParam);//as usually parameters
+
+ for(Int_t iSec = 0; iSec < kNS; iSec++) {
+ AliTPCCalROC * noiseROC;
+ AliTPCCalROC noiseDummy(iSec);
+ if(noiseTPC==0x0){
+ noiseROC = &noiseDummy;//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 (iSec < kNIS) {
+ nRows = fParam->GetNRowLow();
+ nDDLs = 2;
+ indexDDL = iSec * 2;
+ }else {
+ nRows = fParam->GetNRowUp();
+ nDDLs = 4;
+ indexDDL = (iSec-kNIS) * 4 + kNIS * 2;
+ }
+
+ //
+ // Load the raw data for corresponding DDLs
+ //
+ rawReader->Reset();
+ rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
+
+
+ while (input.NextDDL()){
+ // Allocate memory for rows in sector (pads(depends on row) x timebins)
+ if (!digarr->GetRow(iSec,0)){
+ for(Int_t iRow = 0; iRow < nRows; iRow++) {
+ digarr->CreateRow(iSec,iRow);
+ }//end loop over rows
+ }
+ //loop over pads
+ while ( input.NextChannel() ) {
+ Int_t iRow = input.GetRow();
+ Int_t iPad = input.GetPad();
+ //check row consistency
+ if (iRow < 0 ) continue;
+ if (iRow < 0 || iRow >= nRows){
+ AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
+ iRow, 0, nRows -1));
+ continue;
+ }
+
+ //check pad consistency
+ if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) {
+ AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
+ iPad, 0, roc->GetNPads(iSec,iRow) ));
+ continue;
+ }
+
+ //loop over bunches
+ while ( input.NextBunch() ){
+ Int_t startTbin = (Int_t)input.GetStartTimeBin();
+ Int_t bunchlength = (Int_t)input.GetBunchLength();
+ const UShort_t *sig = input.GetSignals();
+ isAltro=kTRUE;
+ for (Int_t iTime = 0; iTime<bunchlength; iTime++){
+ Int_t iTimeBin=startTbin-iTime;
+ //
+ 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 time consistency
+ if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
+ //cout<<iTimeBin<<endl;
+ continue;
+ AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
+ iTimeBin, 0, fRecoParam->GetLastBin() -1));
+ }
+ //signal
+ Float_t signal=(Float_t)sig[iTime];
+ 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
+ }
+ if(signal <= fMaxNoiseSigma * noiseOnPad){
+ digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad);
+ continue;
+ }
+ digarr->GetRow(iSec,iRow)->SetDigitFast(TMath::Nint(signal),iTimeBin,iPad);
+ }// end loop signals in bunch
+ }// end loop bunches
+ } // end loop pads
+ }// end ddl loop
+ }// end sector loop
+ SetDigArr(digarr);
+ if(isAltro) FindClusterKrIO();
+ delete digarr;
+
+ return 0;
+}
+
+void AliTPCclustererKr::CleanSector(Int_t sector){
+ //
+ // clean isolated digits
+ //
+ const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector
+ for(Int_t iRow=0; iRow<kNRows; ++iRow){
+ AliSimDigits *digrow;
+ if(fRawData){
+ digrow = (AliSimDigits*)fDigarr->GetRow(sector,iRow);//real data
+ }else{
+ digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC
+ }
+ if(!digrow) continue;
+ digrow->ExpandBuffer(); //decrunch
+ const Int_t kNPads = digrow->GetNCols(); // number of pads
+ const Int_t kNTime = digrow->GetNRows(); // number of timebins
+ for(Int_t iPad=1;iPad<kNPads-1;iPad++){
+ Short_t* val = digrow->GetDigitsColumn(iPad);
+
+ for(Int_t iTimeBin=1;iTimeBin<kNTime-1;iTimeBin++){
+ if (val[iTimeBin]<=0) continue;
+ if (val[iTimeBin-1]+val[iTimeBin+1]<fIsolCut) {val[iTimeBin]=0; continue;}
+ if (val[iTimeBin-kNTime]+val[iTimeBin+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
+ //
+ if (val[iTimeBin-1-kNTime]+val[iTimeBin+1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
+ if (val[iTimeBin+1-kNTime]+val[iTimeBin-1+kNTime]<fIsolCut) {val[iTimeBin]=0; continue;}
+
+ }
+ }
+ }
+}
+
+
+////____________________________________________________________________________
+Int_t AliTPCclustererKr::FindClusterKrIO()
+{
+
+ //
+ //fParam and fDigarr must be set to run this method
+ //
+
+ Int_t clusterCounter=0;
+ const Int_t nTotalSector=fParam->GetNSector();//number of sectors
+ for(Int_t iSec=0; iSec<nTotalSector; ++iSec){
+ CleanSector(iSec);
+
+ //vector of maxima for each sector
+ //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(iSec,iRow);//real data
+ }else{
+ digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC
+ }
+ if(digrow){//if pointer exist
+ digrow->ExpandBuffer(); //decrunch
+ 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++){
+
+ 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 for fMinAdc
+ 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
+ //Double_t xCord,yCord;
+ //GetXY(iSec,iRow,iPad,xCord,yCord);
+ Double_t x[]={iRow,iPad,iTimeBin};
+ Int_t i[]={iSec};
+ AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+
+ transform->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;
+ }
+
+
+
+
+
+
+ if(ifIncreaseBegin){
+ ifIncreaseBegin=false;
+ increaseBegin=iTimeBin;
+ }
+
+ if(adc>valueMaximum){
+ timeBinMax=iTimeBin;
+ valueMaximum=adc;
+ ifMaximum=true;
+ }
+ sumAdc+=adc;
+ if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge
+ //at least 3 timebins
+ //insert maximum, default values and set flags
+ //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);
+
+ AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+
+ transform->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;
+ }
+
+ }//end loop over timebins
+ }//end loop over pads
+// }else{
+// cout<<"Pointer does not exist!!"<<endl;
+ }//end if poiner exists
+ }//end loop over rows
+
+ MakeClusters(maximaInSector,iSec,clusterCounter);
+ //
+ maximaInSector->SetOwner(kTRUE);
+ maximaInSector->Delete();
+ delete maximaInSector;
+ }//end sector for
+ cout<<"Number of clusters in event: "<<clusterCounter<<endl;
+ return 0;
+}
+
+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();
+
+ nUsedPads++;
+ nUsedRows.push_back((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() ;
+ }
+ }
+ 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();
+
+ clusterKr.SetRMS();//Set pad,row,timebin RMS
+ clusterKr.Set1D();//Set size in pads and timebins
+
+ clusterKr.SetTimeStamp(fTimeStamp);
+ clusterKr.SetRun(fRun);
+
+ clusterCounter++;
+
+
+ //save each cluster into file
+ if (fOutput){
+ (*fOutput)<<"Kr"<<
+ "Cl.="<<&clusterKr<<
+ "\n";
+ }
+ //end of save each cluster into file adc.root
+ }//outer loop
+}
+
+
+
+////____________________________________________________________________________
+
+
+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