From f8aae3779b58ff58faf59c4c1104d068068b06f6 Mon Sep 17 00:00:00 2001 From: hristov Date: Wed, 23 Jul 2003 07:15:03 +0000 Subject: [PATCH] New class AliTPCRawStream for reading raw data, AliTPCclustererMI adapted to NewIO and prepared for raw data input, AliTPCtrackerMI adapted to NewIO, AliTPCTracking.C adapted to NewIO (T.Kuhr) --- TPC/AliTPCCompression.cxx | 6 - TPC/AliTPCCompression.h | 6 + TPC/AliTPCRawStream.cxx | 166 +++++++++++++++++ TPC/AliTPCRawStream.h | 59 +++++++ TPC/AliTPCTracking.C | 362 ++++++++++++++------------------------ TPC/AliTPCclustererMI.cxx | 336 +++++++++++++++++++++++++---------- TPC/AliTPCclustererMI.h | 7 +- TPC/AliTPCtrackerMI.cxx | 63 ++----- TPC/AliTPCtrackerMI.h | 14 +- TPC/TPCLinkDef.h | 3 +- TPC/libTPC.pkg | 1 + 11 files changed, 644 insertions(+), 379 deletions(-) create mode 100644 TPC/AliTPCRawStream.cxx create mode 100644 TPC/AliTPCRawStream.h diff --git a/TPC/AliTPCCompression.cxx b/TPC/AliTPCCompression.cxx index 21b7fa7bf94..05ee15def64 100644 --- a/TPC/AliTPCCompression.cxx +++ b/TPC/AliTPCCompression.cxx @@ -1048,7 +1048,6 @@ Int_t AliTPCCompression::Decompress(AliTPCHNode *RootNode[],const Int_t NumTable cout<<"First one has been found "<1){ - // - bufferFile.FillBuffer(symbol); // //ftxt<Select(0); + fData = new UShort_t[kDataMax]; + fDataSize = fPosition = 0; + fCount = fBunchLength = 0; + + if (!fgRootNode) { + fgRootNode = new AliTPCHNode*[kNumTables]; + fCompression.CreateTreesFromFile(fgRootNode, kNumTables); + } + + fSector = fPrevSector = fRow = fPrevRow = fPad = fPrevPad = fTime = fSignal = -1; +} + +AliTPCRawStream::~AliTPCRawStream() +{ +// clean up + + delete[] fData; +} + + +Bool_t AliTPCRawStream::Next() +{ +// read the next raw digit +// returns kFALSE if there is no digit left + + fPrevSector = fSector; + fPrevRow = fRow; + fPrevPad = fPad; + + while (fCount == 0) { // next trailer + if (fPosition >= fDataSize) { // next payload + UChar_t* data; + do { + if (!fRawReader->ReadNextData(data)) return kFALSE; + } while (fRawReader->GetDataSize() == 0); + + if (fRawReader->IsCompressed()) { // compressed data + ULong_t size = 0; + fCompression.Decompress(fgRootNode, kNumTables, + (char*) data, fRawReader->GetDataSize(), + fData, size); + fDataSize = size; + + } else { // uncompressed data + fDataSize = 0; + Int_t pos = (fRawReader->GetDataSize() * 8) / 10; + while (Get10BitWord(data, pos-1) == 0x2AA) pos--; + while (pos > 0) { + for (Int_t i = 0; i < 4; i++) { // copy trailer + fData[fDataSize++] = Get10BitWord(data, pos-4+i); + } + pos -= 4; + Int_t count = fData[fDataSize-4]; + pos -= (4 - (count % 4)) % 4; // skip fill words + + while (count > 0) { + UShort_t bunchLength = Get10BitWord(data, pos-1); + fData[fDataSize++] = bunchLength; + fData[fDataSize++] = Get10BitWord(data, pos-2); // time bin + + // copy signal amplitudes in increasing order on time + for (Int_t i = 0; i < bunchLength-2; i++) { + fData[fDataSize++] = Get10BitWord(data, pos - bunchLength + i); + } + pos -= bunchLength; + count -= bunchLength; + } + } + } + + fPosition = 0; + } + if (fPosition + 4 >= fDataSize) { + Error("Next", "could not read trailer"); + return kFALSE; + } + fCount = fData[fPosition++]; + fPad = fData[fPosition++]; + fRow = fData[fPosition++]; + fSector = fData[fPosition++]; + fBunchLength = 0; + } + + if (fBunchLength == 0) { + if (fPosition >= fDataSize) { + Error("Next", "could not read bunch length"); + return kFALSE; + } + fBunchLength = fData[fPosition++] - 2; + fCount--; + + if (fPosition >= fDataSize) { + Error("Next", "could not read time bin"); + return kFALSE; + } + fTime = fData[fPosition++] - fBunchLength; + fCount--; + } + + fTime++; + if (fPosition >= fDataSize) { + Error("Next", "could not read sample amplitude"); + return kFALSE; + } + fSignal = fData[fPosition++] + kOffset; + fCount--; + fBunchLength--; + + return kTRUE; +} + + +UShort_t AliTPCRawStream::Get10BitWord(UChar_t* buffer, Int_t position) +// return a word in a 10 bit array as an UShort_t +{ + Int_t iBit = position * 10; + Int_t iByte = iBit / 8; + Int_t shift = iBit % 8; +// return ((buffer[iByte+1] * 256 + buffer[iByte]) >> shift) & 0x03FF; + + // recalculate the byte numbers and the shift because + // the raw data is written as integers where the high bits are filled first + // -> little endian is assumed here ! + Int_t iByteHigh = 4 * (iByte / 4) + 3 - (iByte % 4); + iByte++; + Int_t iByteLow = 4 * (iByte / 4) + 3 - (iByte % 4); + shift = 6 - shift; + return ((buffer[iByteHigh] * 256 + buffer[iByteLow]) >> shift) & 0x03FF; +} diff --git a/TPC/AliTPCRawStream.h b/TPC/AliTPCRawStream.h new file mode 100644 index 00000000000..68781aaa2ce --- /dev/null +++ b/TPC/AliTPCRawStream.h @@ -0,0 +1,59 @@ +#ifndef ALITPCRAWSTREAM_H +#define ALITPCRAWSTREAM_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +#include +#include "AliRawReader.h" +#include "AliTPCCompression.h" + + +class AliTPCRawStream: public TObject { + public : + AliTPCRawStream(AliRawReader* rawReader); + virtual ~AliTPCRawStream(); + + virtual Bool_t Next(); + + inline Int_t GetSector() const {return fSector;}; + inline Int_t GetPrevSector() const {return fPrevSector;}; + inline Bool_t IsNewSector() const {return fSector != fPrevSector;}; + inline Int_t GetRow() const {return fRow;}; + inline Int_t GetPrevRow() const {return fPrevRow;}; + inline Bool_t IsNewRow() const {return (fRow != fPrevRow) || IsNewSector();}; + inline Int_t GetPad() const {return fPad;}; + inline Int_t GetPrevPad() const {return fPrevPad;}; + inline Bool_t IsNewPad() const {return (fPad != fPrevPad) || IsNewRow();}; + inline Int_t GetTime() const {return fTime;}; + inline Int_t GetSignal() const {return fSignal;}; + + protected : + UShort_t Get10BitWord(UChar_t* buffer, Int_t position); + + AliRawReader* fRawReader; // object for reading the raw data + + AliTPCCompression fCompression; // object for decompression + static const Int_t kNumTables = 5; // number of (de)compression tables + static AliTPCHNode** fgRootNode; // (de)compression tables + + static const Int_t kDataMax = 10000000; // size of array for uncompressed raw data + UShort_t* fData; // uncompressed raw data + Int_t fDataSize; // actual size of the uncompressed raw data + Int_t fPosition; // current position in fData + Int_t fCount; // counter of words to be read for current trailer + Int_t fBunchLength; // remaining number of signal bins in the current bunch + static const Int_t kOffset = 1; // offset of signal + + Int_t fSector; // index of current sector + Int_t fPrevSector; // index of previous sector + Int_t fRow; // index of current row + Int_t fPrevRow; // index of previous row + Int_t fPad; // index of current pad + Int_t fPrevPad; // index of previous pad + Int_t fTime; // index of current time bin + Int_t fSignal; // signal in ADC counts + + ClassDef(AliTPCRawStream, 0) // base class for reading TPC raw digits +}; + +#endif diff --git a/TPC/AliTPCTracking.C b/TPC/AliTPCTracking.C index 71656ac8401..c10c0efde3c 100644 --- a/TPC/AliTPCTracking.C +++ b/TPC/AliTPCTracking.C @@ -6,21 +6,22 @@ // author: Jiri Chudoba based on code of Jourij Belikov // version: 1.0 // description: -// reconstructs of tracks in TPC inthe following steps: +// reconstructs of tracks in TPC in the following steps: // TPC cluster finding // TPC track finding // input parameters: -// Int_t nEvents ... nr of events to process +// Int_t nEvents ... nr of events to process (<0 means all) // Int_t firstEventNr ... first event number (starts from 0) -// const char* fileNameHits ... name of file with hits -// const char* fileNameDigits .. name of file with TPC digits -// const char* fileNameClusters .. name of file with TPC clusters (output) -// const char* fileNameTracks .. name of file with TPC tracks (output) -// -// default file names correspond to pp production (2002-04) +// const char* fileName ... name of galice file +// Bool_t makeClusters ... run the cluster finder or not? +// Bool_t makeTracks ... run the track finder or not? +// const char* fileNameRaw ... if not NULL, the cluster finder uses +// the given file as raw data input // // History: // +// 21.07.2003 ... NewIO +// // 18.03.2003 ... Char_t* replaced by const char* // // 03.03.2003 ... SetFieldFactor moved to AliTracker class and @@ -35,263 +36,162 @@ //////////////////////////////////////////////////////////////////////// #if !defined(__CINT__) || defined(__MAKECINT__) -#include "Riostream.h" -#include "TFile.h" -#include "TTree.h" -#include "TBenchmark.h" -#include "AliTPCtracker.h" -#include "AliTPCtrackerMI.h" -#include "AliTPCclusterer.h" -#include "AliTPCclustererMI.h" +#include +#include +#include +#include +#include "AliRunLoader.h" #include "AliTPC.h" -#include "AliRun.h" -#include "AliHeader.h" -#include "AliGenEventHeader.h" -#include "AliTracker.h" - +#include "AliTPCParam.h" +#include "AliTPCclustererMI.h" +#include "AliTPCtrackerMI.h" +#include "AliRawReaderRoot.h" #endif -Int_t gDEBUG = 2; +Bool_t AliTPCTracking(Int_t nEvents = -1, Int_t firstEvent = 0, + const char* fileName = "galice.root", + Bool_t makeClusters = kTRUE, + Bool_t makeTracks = kTRUE, + const char* fileNameRaw = NULL); -Int_t AliTPCTracking(Int_t nEvents=1, Int_t firstEvent=0, - const char* fileNameHits="galice.root", - const char* fileNameDigits="tpc.digits.root", +/* +Int_t TPCRefitInward(Int_t nEvents=1, Int_t firstEvent=0, const char* fileNameClusters="tpc.clusters.root", const char* fileNameTracks="tpc.tracks.root", - Bool_t versionMI = kFALSE); - -Int_t TPCFindClusters(Int_t nEvents=1, Int_t firstEvent=0, - const char* fileNameDigits="tpc.digits.root", - const char* fileNameClusters="tpc.clusters.root", - Bool_t versionMI = kFALSE); -Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent, - TFile* fileDigits, TFile* fileClusters, - AliTPCParam* paramTPC=0); -Int_t TPCFindClustersMI(Int_t nEvents, Int_t firstEvent, - TFile* fileDigits, TFile* fileClusters, - AliTPCParam* paramTPC=0); -Int_t TPCFindTracks(Int_t nEvents=1, Int_t firstEvent=0, - const char* fileNameClusters="tpc.clusters.root", - const char* fileNameTracks="tpc.tracks.root", - Bool_t versionMI = kFALSE); -Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent, - TFile* fileClusters, TFile* fileTracks, - AliTPCParam* paramTPC=0); -Int_t TPCFindTracksMI(Int_t nEvents, Int_t firstEvent, - TFile* fileClusters, TFile* fileTracks, - AliTPCParam* paramTPC=0); - -void FindVertex(Int_t iEvent, Double_t *vertex); -void PrintVertex(TArrayF &primaryVertex); + const char* fileNameTRDTracks="trd.tracks.root", + const char* fileNameRefittedTracks="tpc.refitted.tracks.root"); +*/ //////////////////////////////////////////////////////////////////////// -Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent, - const char* fileNameHits, - const char* fileNameDigits, - const char* fileNameClusters, - const char* fileNameTracks, - Bool_t versionMI) { - - AliTracker::SetFieldFactor(fileNameHits,kFALSE); - -// ********** Find TPC clusters *********** // - if (fileNameDigits && fileNameClusters) { - if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters,versionMI)) { - cerr<<"Failed to get TPC clusters: !\n"; - return 1; - } - } - -// ********** Find TPC tracks *********** // - if (fileNameClusters && fileNameTracks) { - if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks,versionMI)) { - cerr<<"Failed to get TPC tracks !\n"; - return 2; - } +Bool_t AliTPCTracking(Int_t nEvents, Int_t firstEvent, + const char* fileName, + Bool_t makeClusters, + Bool_t makeTracks, + const char* fileNameRaw) +{ + // get the loaders + AliRunLoader* runLoader = AliRunLoader::Open(fileName); + if (!runLoader) { + cerr << "AliTPCTracking: no run loader found\n"; + return kFALSE; + } + AliLoader* tpcLoader = runLoader->GetLoader("TPCLoader"); + if (!tpcLoader) { + cerr << "AliTPCTracking: no TPC loader found\n"; + return kFALSE; } - return 0; -} + // get the TPC parameters + runLoader->CdGAFile(); + AliTPCParam* param = AliTPC::LoadTPCParam(gFile); + if (!param) { + cerr << "AliTPCTracking: no TPC parameters found\n"; + return kFALSE; + } -//////////////////////////////////////////////////////////////////////// -Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent, - const char* fileNameDigits, const char* fileNameClusters, - Bool_t versionMI) { - - Int_t rc; - const Char_t *name="TPCFindClusters"; - if (gDEBUG>1) cout<1) gBenchmark->Start(name); - TFile *fileClusters = TFile::Open(fileNameClusters,"recreate"); - TFile *fileDigits = TFile::Open(fileNameDigits); - if (!fileDigits->IsOpen()) { - cerr<<"Cannnot open "<LoadDigits(); + tpcLoader->LoadRecPoints("recreate"); } - if (!fileClusters->IsOpen()) { - cerr<<"Cannnot open "<LoadRecPoints(); + tpcLoader->LoadTracks("recreate"); } - if (versionMI) { - rc = TPCFindClustersMI(nEvents,firstEvent,fileDigits,fileClusters); + // get the event number range + Int_t maxEvent = 0; + if (fileNameRaw) { + TFile* file = TFile::Open(fileNameRaw); + if (file && file->IsOpen()) { + TTree* tree = (TTree*) file->Get("T"); + if (tree) maxEvent = (Int_t) tree->GetEntries(); + } } else { - rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters); + maxEvent = runLoader->GetNumberOfEvents(); } + if (nEvents < 0) nEvents = maxEvent - firstEvent; + Int_t lastEvent = firstEvent + nEvents; + if (lastEvent > maxEvent) lastEvent = maxEvent; + + // loop over the events + for (Int_t iEvent = firstEvent; iEvent < lastEvent; iEvent++) { + + runLoader->GetEvent(iEvent); + + // run the cluster finder + if (makeClusters) { + if (!tpcLoader->TreeR()) tpcLoader->MakeRecPointsContainer(); + clusterer->SetOutput(tpcLoader->TreeR()); + if (fileNameRaw) { + AliRawReaderRoot rawReader(fileNameRaw, iEvent); + clusterer->Digits2Clusters(&rawReader); + } else { + clusterer->SetInput(tpcLoader->TreeD()); + clusterer->Digits2Clusters(); + } + tpcLoader->WriteRecPoints("OVERWRITE"); + } - fileDigits->Close(); - fileClusters->Close(); - delete fileDigits; - delete fileClusters; - if (gDEBUG>1) gBenchmark->Stop(name); - if (gDEBUG>1) gBenchmark->Show(name); - - return rc; -} -//////////////////////////////////////////////////////////////////////// -Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent, - TFile* fileDigits, TFile* fileClusters, - AliTPCParam* paramTPC) { - - fileDigits->cd(); - if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileDigits); - if (!paramTPC) return 1; - - for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){ - if (gDEBUG > 2) cout<<"TPCFindClusters: event "<Clusters2Tracks(); + delete tracker; + } } - return 0; -} -//////////////////////////////////////////////////////////////////////// -Int_t TPCFindClustersMI(Int_t nEvents, Int_t firstEvent, - TFile* fileDigits, TFile* fileClusters, - AliTPCParam* paramTPC) { - fileDigits->cd(); - if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileDigits); - if (!paramTPC) return 1; + if (tracker) delete tracker; + if (clusterer) delete clusterer; - AliTPCclustererMI clusterer; - for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){ - if (gDEBUG > 2) cout<<"TPCFindClustersMI: event "<Get(treeName); - fileClusters->cd(); - sprintf(treeName, "TreeC_TPC_%d", iEvent); - TTree* output = new TTree(treeName, treeName); - clusterer.SetInput(input); - clusterer.SetOutput(output); - clusterer.Digits2Clusters(paramTPC, iEvent); - } - - return 0; + return kTRUE; } -//////////////////////////////////////////////////////////////////////// -Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent, - const char* fileNameClusters, const char* fileNameTracks, - Bool_t versionMI) { +/* +//////////////////////////////////////////////////////////////////////// +Int_t TPCRefitInward(Int_t nEvents, Int_t firstEvent, + const char* fileNameClusters, + const char* fileNameTracks, + const char* fileNameTRDTracks, + const char* fileNameRefittedTracks) +{ Int_t rc = 0; - const Char_t *name="TPCFindTracks"; + const Char_t *name="TPCRefitInward"; if (gDEBUG>1) cout<1) gBenchmark->Start(name); - TFile *fileTracks = TFile::Open(fileNameTracks,"recreate"); - TFile *fileClusters =TFile::Open(fileNameClusters); - - if (versionMI) { - rc = TPCFindTracksMI(nEvents, firstEvent, fileClusters, fileTracks); - } else { - rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks); - } - - fileClusters->Close(); - fileTracks->Close(); - delete fileClusters; - delete fileTracks; - if (gDEBUG>1) gBenchmark->Stop(name); - if (gDEBUG>1) gBenchmark->Show(name); - return rc; - -} -//////////////////////////////////////////////////////////////////////// -Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent, - TFile *fileClusters, TFile * fileTracks, - AliTPCParam* paramTPC) { + TFile *fileClusters = TFile::Open(fileNameClusters); + TFile *fileTracks = TFile::Open(fileNameTracks); + TFile *fileTRDTracks = TFile::Open(fileNameTRDTracks); + TFile *fileRefittedTracks = TFile::Open(fileNameRefittedTracks, "recreate"); - Int_t rc = 0; - if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileClusters); + AliTPCParam* paramTPC = AliTPC::LoadTPCParam(fileClusters); if (!paramTPC) return 1; for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){ - if (gDEBUG > 2) cout<<"TPCFindTracks: event "< 2) cout<<"TPCRefitInward: event "<SetEventNumber(iEvent); - Double_t vertex[3]; - FindVertex(iEvent,vertex); - tracker->SetVertex(vertex); fileClusters->cd(); - rc = tracker->Clusters2Tracks(0,fileTracks); + rc = tracker->RefitInward(fileTRDTracks, fileTracks, fileRefittedTracks); delete tracker; + if (rc) return rc; } - return rc; -} -//////////////////////////////////////////////////////////////////////// -Int_t TPCFindTracksMI(Int_t nEvents, Int_t firstEvent, - TFile *fileClusters, TFile * fileTracks, - AliTPCParam* paramTPC) { - Int_t rc = 0; - if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileClusters); - if (!paramTPC) return 1; - - for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){ - if (gDEBUG > 2) cout<<"TPCFindTracksMI: event "<Close(); + fileTracks->Close(); + fileTRDTracks->Close(); + fileRefittedTracks->Close(); + delete fileClusters; + delete fileTracks; + delete fileTRDTracks; + delete fileRefittedTracks; + if (gDEBUG>1) gBenchmark->Show(name); return rc; } - -//////////////////////////////////////////////////////////////////////// -void FindVertex(Int_t eventNr, Double_t *vertex) { - - vertex[0] = vertex[1] = vertex[2] = 0.; - if (!gAlice) { - cerr<<"gAlice was not found! Using vertex position (0,0,0).\n"; - return; - } - - gAlice->GetEvent(eventNr); - AliHeader *header = gAlice->GetHeader(); - if (!header) { - cerr<<"header was not found!\n"; - return; - } - AliGenEventHeader* genEventHeader = header->GenEventHeader(); - if (!genEventHeader) { - cerr<<"AliGenEventHeader was not found!\n"; - return; - } - - TArrayF primaryVertex(3); - genEventHeader->PrimaryVertex(primaryVertex); - PrintVertex(primaryVertex); - vertex[0] = static_cast(primaryVertex[0]); - vertex[1] = static_cast(primaryVertex[1]); - vertex[2] = static_cast(primaryVertex[2]); -// delete header; - delete genEventHeader; - return; - -} -//////////////////////////////////////////////////////////////////////// -void PrintVertex(TArrayF &primaryVertex) -{ - cout <<"Vertex: " - < #include "AliTPCClustersArray.h" #include "AliTPCClustersRow.h" +#include "AliTPCRawStream.h" #include "AliDigits.h" #include "AliSimDigits.h" #include "AliTPCParam.h" +#include "AliRawReader.h" +#include "AliTPCRawStream.h" +#include "AliRunLoader.h" +#include "AliLoader.h" #include "Riostream.h" #include @@ -37,10 +42,11 @@ ClassImp(AliTPCclustererMI) -AliTPCclustererMI::AliTPCclustererMI() +AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par) { fInput =0; fOutput=0; + fParam = par; } void AliTPCclustererMI::SetInput(TTree * tree) { @@ -403,9 +409,11 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){ if (kj<0) kj=0; if (kj>=fMaxTime-3) kj=fMaxTime-4; // ki and kj shifted to "real" coordinata - c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0); - c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1); - c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2); + if (fRowDig) { + c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0); + c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1); + c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2); + } Float_t s2 = c.GetSigmaY2(); @@ -436,21 +444,20 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){ //_____________________________________________________________________________ -void AliTPCclustererMI::Digits2Clusters(const AliTPCParam *par, Int_t eventn) +void AliTPCclustererMI::Digits2Clusters() { //----------------------------------------------------------------- // This is a simple cluster finder. //----------------------------------------------------------------- - TDirectory *savedir=gDirectory; - if (fInput==0){ - cerr<<"AliTPC::Digits2Clusters(): input tree not initialised !\n"; + if (!fInput) { + Error("Digits2Clusters", "input tree not initialised"); return; } - if (fOutput==0) { - cerr<<"AliTPC::Digits2Clusters(): output tree not initialised !\n"; - return; + if (!fOutput) { + Error("Digits2Clusters", "output tree not initialised"); + return; } AliSimDigits digarr, *dummy=&digarr; @@ -458,17 +465,14 @@ void AliTPCclustererMI::Digits2Clusters(const AliTPCParam *par, Int_t eventn) fInput->GetBranch("Segment")->SetAddress(&dummy); Stat_t nentries = fInput->GetEntries(); - fMaxTime=par->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after + fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after - fParam = par; - ((AliTPCParam*)par)->Write(par->GetTitle()); - Int_t nclusters = 0; for (Int_t n=0; nGetEvent(n); Int_t row; - if (!par->AdjustSectorRow(digarr.GetID(),fSector,row)) { + if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,row)) { cerr<<"AliTPC warning: invalid segment ID ! "<SetID(digarr.GetID()); fOutput->GetBranch("Segment")->SetAddress(&clrow); - fRx=par->GetPadRowRadii(fSector,row); + fRx=fParam->GetPadRowRadii(fSector,row); - const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector(); + const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector(); fZWidth = fParam->GetZWidth(); if (fSector < kNIS) { - fMaxPad = par->GetNPadsLow(row); + fMaxPad = fParam->GetNPadsLow(row); fSign = (fSector < kNIS/2) ? 1 : -1; - fPadLength = par->GetPadPitchLength(fSector,row); - fPadWidth = par->GetPadPitchWidth(); + fPadLength = fParam->GetPadPitchLength(fSector,row); + fPadWidth = fParam->GetPadPitchWidth(); } else { - fMaxPad = par->GetNPadsUp(row); + fMaxPad = fParam->GetNPadsUp(row); fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; - fPadLength = par->GetPadPitchLength(fSector,row); - fPadWidth = par->GetPadPitchWidth(); + fPadLength = fParam->GetPadPitchLength(fSector,row); + fPadWidth = fParam->GetPadPitchWidth(); } @@ -506,75 +510,13 @@ void AliTPCclustererMI::Digits2Clusters(const AliTPCParam *par, Int_t eventn) if (digarr.First()) //MI change do { Short_t dig=digarr.CurrentDigit(); - if (dig<=par->GetZeroSup()) continue; + if (dig<=fParam->GetZeroSup()) continue; Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3; fBins[i*fMaxTime+j]=dig; } while (digarr.Next()); digarr.ExpandTrackBuffer(); - //add virtual charge at the edge - for (Int_t i=0; i0){ - Float_t amp2 = fBins[i+4*fMaxTime]; - if (amp2==0) amp2=0.5; - Float_t sigma2 = GetSigmaY2(i); - amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); - if (gDebug>4) printf("\n%f\n",amp0); - } - fBins[i+2*fMaxTime] = Int_t(amp0); - amp0 = 0; - amp1 = fBins[(fMaxPad+2)*fMaxTime+i]; - if (amp1>0){ - Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime]; - if (amp2==0) amp2=0.5; - Float_t sigma2 = GetSigmaY2(i); - amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); - if (gDebug>4) printf("\n%f\n",amp0); - } - fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0); - } - - memcpy(fResBins,fBins, fMaxBin*2); - // - fNcluster=0; - //first loop - for "gold cluster" - fLoop=1; - Int_t *b=&fBins[-1]+2*fMaxTime; - Int_t crtime = Int_t((fParam->GetZLength()-1.05*fRx)/fZWidth-5); - - for (Int_t i=2*fMaxTime; iFill(); delete clrow; @@ -582,8 +524,220 @@ void AliTPCclustererMI::Digits2Clusters(const AliTPCParam *par, Int_t eventn) delete[] fBins; delete[] fResBins; } - cerr<<"Number of found clusters : "<Reset(); + AliTPCRawStream input(rawReader); + + fRowDig = NULL; + + Int_t nclusters = 0; - fOutput->Write(); - savedir->cd(); + fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after + const Int_t kNIS = fParam->GetNInnerSector(); + const Int_t kNOS = fParam->GetNOuterSector(); + const Int_t kNS = kNIS + kNOS; + fZWidth = fParam->GetZWidth(); + Int_t zeroSup = fParam->GetZeroSup(); + + fBins = NULL; + Int_t** splitRows = new Int_t* [kNS*2]; + Int_t** splitRowsRes = new Int_t* [kNS*2]; + for (Int_t iSector = 0; iSector < kNS*2; iSector++) + splitRows[iSector] = NULL; + Int_t iSplitRow = -1; + + Bool_t next = kTRUE; + while (next) { + next = input.Next(); + + // when the sector or row number has changed ... + if (input.IsNewRow() || !next) { + + // ... find clusters in the previous pad row, and ... + if (fBins) { + if ((iSplitRow < 0) || splitRows[fSector + kNS*iSplitRow]) { + fRowCl = new AliTPCClustersRow; + fRowCl->SetClass("AliTPCclusterMI"); + fRowCl->SetArray(1); + fRowCl->SetID(fParam->GetIndex(fSector, input.GetPrevRow())); + fOutput->GetBranch("Segment")->SetAddress(&fRowCl); + + FindClusters(); + + fOutput->Fill(); + delete fRowCl; + nclusters += fNcluster; + delete[] fBins; + delete[] fResBins; + if (iSplitRow >= 0) splitRows[fSector + kNS*iSplitRow] = NULL; + + } else if (iSplitRow >= 0) { + splitRows[fSector + kNS*iSplitRow] = fBins; + splitRowsRes[fSector + kNS*iSplitRow] = fResBins; + } + } + + if (!next) break; + + // ... prepare for the next pad row + fSector = input.GetSector(); + Int_t iRow = input.GetRow(); + fRx = fParam->GetPadRowRadii(fSector, iRow); + + iSplitRow = -1; + if (fSector < kNIS) { + fMaxPad = fParam->GetNPadsLow(iRow); + fSign = (fSector < kNIS/2) ? 1 : -1; + if (iRow == 30) iSplitRow = 0; + } else { + fMaxPad = fParam->GetNPadsUp(iRow); + fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; + if (iRow == 27) iSplitRow = 0; + else if (iRow == 76) iSplitRow = 1; + } + fPadLength = fParam->GetPadPitchLength(fSector, iRow); + fPadWidth = fParam->GetPadPitchWidth(); + + fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after + if ((iSplitRow < 0) || !splitRows[fSector + kNS*iSplitRow]) { + fBins = new Int_t[fMaxBin]; + fResBins = new Int_t[fMaxBin]; //fBins with residuals after 1 finder loop + memset(fBins, 0, sizeof(Int_t)*fMaxBin); + } else { + fBins = splitRows[fSector + kNS*iSplitRow]; + fResBins = splitRowsRes[fSector + kNS*iSplitRow]; + } + } + + // fill fBins with digits data + if (input.GetSignal() <= zeroSup) continue; + Int_t i = input.GetPad() + 3; + Int_t j = input.GetTime() + 3; + fBins[i*fMaxTime+j] = input.GetSignal(); + } + + // find clusters in split rows that were skipped until now. + // this can happen if the rows were not splitted + for (fSector = 0; fSector < kNS; fSector++) + for (Int_t iSplit = 0; iSplit < 2; iSplit++) + if (splitRows[fSector + kNS*iSplit]) { + + Int_t iRow = -1; + if (fSector < kNIS) { + iRow = 30; + fMaxPad = fParam->GetNPadsLow(iRow); + fSign = (fSector < kNIS/2) ? 1 : -1; + } else { + if (iSplit == 0) iRow = 27; else iRow = 76; + fMaxPad = fParam->GetNPadsUp(iRow); + fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; + } + fRx = fParam->GetPadRowRadii(fSector, iRow); + fPadLength = fParam->GetPadPitchLength(fSector, iRow); + fPadWidth = fParam->GetPadPitchWidth(); + + fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after + fBins = splitRows[fSector + kNS*iSplit]; + fResBins = splitRowsRes[fSector + kNS*iSplit]; + + fRowCl = new AliTPCClustersRow; + fRowCl->SetClass("AliTPCclusterMI"); + fRowCl->SetArray(1); + fRowCl->SetID(fParam->GetIndex(fSector, iRow)); + fOutput->GetBranch("Segment")->SetAddress(&fRowCl); + + FindClusters(); + + fOutput->Fill(); + delete fRowCl; + nclusters += fNcluster; + delete[] fBins; + delete[] fResBins; + } + + delete[] splitRows; + delete[] splitRowsRes; + Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters); +} + +void AliTPCclustererMI::FindClusters() +{ + //add virtual charge at the edge + for (Int_t i=0; i0){ + Float_t amp2 = fBins[i+4*fMaxTime]; + if (amp2==0) amp2=0.5; + Float_t sigma2 = GetSigmaY2(i); + amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); + if (gDebug>4) printf("\n%f\n",amp0); + } + fBins[i+2*fMaxTime] = Int_t(amp0); + amp0 = 0; + amp1 = fBins[(fMaxPad+2)*fMaxTime+i]; + if (amp1>0){ + Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime]; + if (amp2==0) amp2=0.5; + Float_t sigma2 = GetSigmaY2(i); + amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); + if (gDebug>4) printf("\n%f\n",amp0); + } + fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0); + } + +// memcpy(fResBins,fBins, fMaxBin*2); + memcpy(fResBins,fBins, fMaxBin); + // + fNcluster=0; + //first loop - for "gold cluster" + fLoop=1; + Int_t *b=&fBins[-1]+2*fMaxTime; + Int_t crtime = Int_t((fParam->GetZLength()-1.05*fRx)/fZWidth-5); + + for (Int_t i=2*fMaxTime; iGetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2) { //--------------------------------------------------------------------- @@ -216,17 +217,6 @@ AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2) fClustersArray.Setup(par); fClustersArray.SetClusterType("AliTPCclusterMI"); - char cname[100]; - if (eventn==-1) { - sprintf(cname,"TreeC_TPC"); - } - else { - sprintf(cname,"TreeC_TPC_%d",eventn); - } - - fClustersArray.ConnectTree(cname); - - fEventN = eventn; fSeeds=0; fNtracks = 0; fParam = par; @@ -608,6 +598,10 @@ Int_t AliTPCtrackerMI::LoadClusters() for (Int_t i=0; iIsOpen()) { - cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n"; - return 1; - } - } - - if (!out->IsOpen()) { - cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n"; - return 2; + TTree* clustersTree = AliRunLoader::GetTreeR("TPC", kFALSE); + if (!clustersTree) { + Error("Clusters2Tracks", "no clusters found"); + return 1; } + fClustersArray.ConnectTree(clustersTree); - out->cd(); - - char tname[100]; - sprintf(tname,"TreeT_TPC_%d",fEventN); - TTree tracktree(tname,"Tree with TPC tracks"); - TTree seedtree("Seeds","Seeds"); + TTree* tracksTree = AliRunLoader::GetTreeT("TPC", kTRUE); + TTree& tracktree = *tracksTree; +// TTree seedtree("Seeds","Seeds"); AliTPCtrack *iotrack=0; AliTPCseed *ioseed=0; tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0); @@ -1775,13 +1758,6 @@ Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) { LoadClusters(); printf("Time for loading clusters: \t");timer.Print();timer.Start(); - printf("Loading outer sectors\n"); - LoadOuterSectors(); - printf("Time for loading outer sectors: \t");timer.Print();timer.Start(); - - printf("Loading inner sectors\n"); - LoadInnerSectors(); - printf("Time for loading inner sectors: \t");timer.Print();timer.Start(); fSectors = fOuterSec; fN=fkNOS; @@ -1836,7 +1812,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) { vseed->fPoints->ExpandCreateFast(2); //TBranch * seedbranch = - seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99); +// seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99); //delete vseed; nseed=fSeeds->GetEntriesFast(); @@ -1873,12 +1849,11 @@ Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) { UnloadClusters(); printf("Time for unloading cluster: \t"); timer.Print();timer.Start(); - tracktree.Write(); - seedtree.Write(); +// seedtree.Write(); cerr<<"Number of found tracks : "<<"\t"<cd(); - + AliRunLoader::GetDetectorLoader("TPC")->WriteTracks("OVERWRITE"); + return 0; } diff --git a/TPC/AliTPCtrackerMI.h b/TPC/AliTPCtrackerMI.h index e4f46b150c8..23caddef30b 100644 --- a/TPC/AliTPCtrackerMI.h +++ b/TPC/AliTPCtrackerMI.h @@ -126,19 +126,25 @@ public: AliTPCtrackerMI():AliTracker(),fkNIS(0),fkNOS(0) { fInnerSec=fOuterSec=0; fSeeds=0; } - AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn=0); + AliTPCtrackerMI(const AliTPCParam *par); ~AliTPCtrackerMI(); + virtual Int_t Clusters2Tracks(AliESD *event) {return -1;}; + virtual Int_t PropagateBack(AliESD *event) {return -1;}; + virtual Int_t RefitInward(AliESD *event) {return -1;}; + virtual Int_t LoadClusters(TTree *) {return -1;}; + Int_t ReadSeeds(const TFile *in); Int_t LoadClusters(); void UnloadClusters(); void LoadInnerSectors(); void LoadOuterSectors(); - AliCluster * GetCluster (int) const {return 0;} + AliCluster * GetCluster (Int_t index) const + {return (AliCluster*) GetClusterMI(index);} AliTPCclusterMI *GetClusterMI(Int_t index) const; - Int_t Clusters2Tracks(const TFile *in, TFile *out); - // Int_t PropagateBack(const TFile *in, TFile *out); + Int_t Clusters2Tracks(); + Int_t PropagateBack() {return -1;}; virtual void CookLabel(AliKalmanTrack *t,Float_t wrong) const; void RotateToLocal(AliTPCseed *seed); diff --git a/TPC/TPCLinkDef.h b/TPC/TPCLinkDef.h index 7f5e277e4e8..a78246831e7 100644 --- a/TPC/TPCLinkDef.h +++ b/TPC/TPCLinkDef.h @@ -17,7 +17,6 @@ #pragma link C++ class AliTPChit+; #pragma link C++ class AliTPCdigit+; #pragma link C++ class AliTPCcluster+; -#pragma link C++ class AliTPCclusterer+; #pragma link C++ class AliTPCtrack+; #pragma link C++ class AliTPCtracker+; @@ -56,6 +55,7 @@ #pragma link C++ class AliHitInfo+; +#pragma link C++ class AliTPCclusterer+; #pragma link C++ class AliTPCDigitizer; #pragma link C++ class AliTPCtrackerParam; @@ -81,6 +81,7 @@ #pragma link C++ class AliTPCBuffer160+; #pragma link C++ class AliTPCCompression+; #pragma link C++ class AliTPCDDLRawData+; +#pragma link C++ class AliTPCRawStream+; #pragma link C++ class AliTPCHNode+; #pragma link C++ class AliTPCHTable+; #pragma link C++ class AliTPCKalmanSegment+; diff --git a/TPC/libTPC.pkg b/TPC/libTPC.pkg index 29acc3a9619..d6ad10a39a0 100644 --- a/TPC/libTPC.pkg +++ b/TPC/libTPC.pkg @@ -16,6 +16,7 @@ SRCS:= AliTPCLoader.cxx\ AliTPCpolyTrack.cxx \ AliTPCBuffer.cxx AliTPCBuffer160.cxx \ AliTPCCompression.cxx AliTPCDDLRawData.cxx \ + AliTPCRawStream.cxx \ AliTPCHuffman.cxx AliTPCPid.cxx AliTPCtrackPid.cxx AliTPCpidESD.cxx HDRS:= $(SRCS:.cxx=.h) -- 2.39.3