X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TRD%2FAliTRDclusterizer.cxx;h=6d57608efaf734bce1852d807809700ea2f64c2a;hp=2e7c152593b7e8671f2634236989953a97556a05;hb=7377baa39e2029aef48cbaedfbe1b6003c95da7c;hpb=6beac890e3f746c6cb0ad1843a5ca395ca109700 diff --git a/TRD/AliTRDclusterizer.cxx b/TRD/AliTRDclusterizer.cxx index 2e7c152593b..6d57608efaf 100644 --- a/TRD/AliTRDclusterizer.cxx +++ b/TRD/AliTRDclusterizer.cxx @@ -26,6 +26,7 @@ #include "AliRunLoader.h" #include "AliLoader.h" +#include "AliTreeLoader.h" #include "AliAlignObj.h" #include "AliTRDclusterizer.h" @@ -40,15 +41,19 @@ #include "AliTRDcalibDB.h" #include "AliTRDtransform.h" #include "AliTRDSignalIndex.h" -#include "AliTRDrawStreamBase.h" +#include "AliTRDrawStream.h" #include "AliTRDfeeParam.h" #include "AliTRDtrackletWord.h" +#include "AliTRDtrackletMCM.h" +#include "AliTRDtrackGTU.h" +#include "AliESDTrdTrack.h" #include "TTreeStream.h" #include "Cal/AliTRDCalROC.h" #include "Cal/AliTRDCalDet.h" #include "Cal/AliTRDCalSingleChamberStatus.h" +#include "Cal/AliTRDCalOnlineGainTableROC.h" ClassImp(AliTRDclusterizer) @@ -60,14 +65,15 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec) ,fClusterTree(NULL) ,fRecPoints(NULL) ,fTracklets(NULL) - ,fTrackletTree(NULL) + ,fTracks(NULL) ,fDigitsManager(new AliTRDdigitsManager()) - ,fTrackletContainer(NULL) ,fRawVersion(2) ,fTransform(new AliTRDtransform(0)) ,fDigits(NULL) + ,fDigitsRaw(NULL) ,fIndexes(NULL) ,fMaxThresh(0) + ,fMaxThreshTest(0) ,fSigThresh(0) ,fMinMaxCutSigma(0) ,fMinLeftRightCutSigma(0) @@ -81,20 +87,20 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec) ,fCalNoiseROC(NULL) ,fCalNoiseDetValue(0) ,fCalPadStatusROC(NULL) + ,fCalOnlGainROC(NULL) ,fClusterROC(0) ,firstClusterROC(0) ,fNoOfClusters(0) + ,fBaseline(0) + ,fRawStream(NULL) + ,fTrgFlags() { // // AliTRDclusterizer default constructor // SetBit(kLabels, kTRUE); - - AliTRDcalibDB *trd = 0x0; - if (!(trd = AliTRDcalibDB::Instance())) { - AliFatal("Could not get calibration object"); - } + SetBit(knewDM, kFALSE); fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion(); @@ -110,21 +116,24 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec) } //_____________________________________________________________________________ -AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, const AliTRDReconstructor *const rec) +AliTRDclusterizer::AliTRDclusterizer(const Text_t *name + , const Text_t *title + , const AliTRDReconstructor *const rec) :TNamed(name,title) ,fReconstructor(rec) ,fRunLoader(NULL) ,fClusterTree(NULL) ,fRecPoints(NULL) ,fTracklets(NULL) - ,fTrackletTree(NULL) + ,fTracks(NULL) ,fDigitsManager(new AliTRDdigitsManager()) - ,fTrackletContainer(NULL) ,fRawVersion(2) ,fTransform(new AliTRDtransform(0)) ,fDigits(NULL) + ,fDigitsRaw(NULL) ,fIndexes(NULL) ,fMaxThresh(0) + ,fMaxThreshTest(0) ,fSigThresh(0) ,fMinMaxCutSigma(0) ,fMinLeftRightCutSigma(0) @@ -138,20 +147,20 @@ AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, co ,fCalNoiseROC(NULL) ,fCalNoiseDetValue(0) ,fCalPadStatusROC(NULL) + ,fCalOnlGainROC(NULL) ,fClusterROC(0) ,firstClusterROC(0) ,fNoOfClusters(0) + ,fBaseline(0) + ,fRawStream(NULL) + ,fTrgFlags() { // // AliTRDclusterizer constructor // SetBit(kLabels, kTRUE); - - AliTRDcalibDB *trd = 0x0; - if (!(trd = AliTRDcalibDB::Instance())) { - AliFatal("Could not get calibration object"); - } + SetBit(knewDM, kFALSE); fDigitsManager->CreateArrays(); @@ -169,14 +178,15 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c) ,fClusterTree(NULL) ,fRecPoints(NULL) ,fTracklets(NULL) - ,fTrackletTree(NULL) + ,fTracks(NULL) ,fDigitsManager(NULL) - ,fTrackletContainer(NULL) ,fRawVersion(2) ,fTransform(NULL) ,fDigits(NULL) + ,fDigitsRaw(NULL) ,fIndexes(NULL) ,fMaxThresh(0) + ,fMaxThreshTest(0) ,fSigThresh(0) ,fMinMaxCutSigma(0) ,fMinLeftRightCutSigma(0) @@ -190,15 +200,20 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c) ,fCalNoiseROC(NULL) ,fCalNoiseDetValue(0) ,fCalPadStatusROC(NULL) + ,fCalOnlGainROC(NULL) ,fClusterROC(0) ,firstClusterROC(0) ,fNoOfClusters(0) + ,fBaseline(0) + ,fRawStream(NULL) + ,fTrgFlags() { // // AliTRDclusterizer copy constructor // SetBit(kLabels, kTRUE); + SetBit(knewDM, kFALSE); //FillLUT(); @@ -211,31 +226,25 @@ AliTRDclusterizer::~AliTRDclusterizer() // AliTRDclusterizer destructor // - if (fRecPoints/* && IsClustersOwner()*/){ - fRecPoints->Delete(); - delete fRecPoints; - } - - if (fTracklets){ - fTracklets->Delete(); - delete fTracklets; - } - if (fDigitsManager) { delete fDigitsManager; fDigitsManager = NULL; } - if (fTrackletContainer){ - delete fTrackletContainer; - fTrackletContainer = NULL; + if (fDigitsRaw) { + delete fDigitsRaw; + fDigitsRaw = NULL; } if (fTransform){ delete fTransform; - fTransform = NULL; + fTransform = NULL; } + if (fRawStream){ + delete fRawStream; + fRawStream = NULL; + } } //_____________________________________________________________________________ @@ -263,14 +272,14 @@ void AliTRDclusterizer::Copy(TObject &c) const ((AliTRDclusterizer &) c).fClusterTree = NULL; ((AliTRDclusterizer &) c).fRecPoints = NULL; - ((AliTRDclusterizer &) c).fTrackletTree = NULL; ((AliTRDclusterizer &) c).fDigitsManager = NULL; - ((AliTRDclusterizer &) c).fTrackletContainer = NULL; ((AliTRDclusterizer &) c).fRawVersion = fRawVersion; ((AliTRDclusterizer &) c).fTransform = NULL; - ((AliTRDclusterizer &) c).fDigits = NULL; + ((AliTRDclusterizer &) c).fDigits = NULL; + ((AliTRDclusterizer &) c).fDigitsRaw = NULL; ((AliTRDclusterizer &) c).fIndexes = NULL; ((AliTRDclusterizer &) c).fMaxThresh = 0; + ((AliTRDclusterizer &) c).fMaxThreshTest = 0; ((AliTRDclusterizer &) c).fSigThresh = 0; ((AliTRDclusterizer &) c).fMinMaxCutSigma= 0; ((AliTRDclusterizer &) c).fMinLeftRightCutSigma = 0; @@ -287,6 +296,9 @@ void AliTRDclusterizer::Copy(TObject &c) const ((AliTRDclusterizer &) c).fClusterROC = 0; ((AliTRDclusterizer &) c).firstClusterROC= 0; ((AliTRDclusterizer &) c).fNoOfClusters = 0; + ((AliTRDclusterizer &) c).fBaseline = 0; + ((AliTRDclusterizer &) c).fRawStream = NULL; + } //_____________________________________________________________________________ @@ -352,46 +364,6 @@ Bool_t AliTRDclusterizer::OpenOutput(TTree *const clusterTree) return kTRUE; } -//_____________________________________________________________________________ -Bool_t AliTRDclusterizer::OpenTrackletOutput() -{ - // - // Tracklet writing - // - - if (fReconstructor->IsWritingTracklets()){ - TString evfoldname = AliConfig::GetDefaultEventFolderName(); - fRunLoader = AliRunLoader::GetRunLoader(evfoldname); - - if (!fRunLoader) { - fRunLoader = AliRunLoader::Open("galice.root"); - } - if (!fRunLoader) { - AliError(Form("Can not open session for file galice.root.")); - return kFALSE; - } - - UInt_t **leaves = new UInt_t *[2]; - AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets"); - if (!dl) { - AliError("Could not get the tracklets data loader!"); - dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets"); - fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl); - } - fTrackletTree = dl->Tree(); - if (!fTrackletTree) - { - dl->MakeTree(); - fTrackletTree = dl->Tree(); - } - TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch"); - if (!trkbranch) - fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i"); - } - - return kTRUE; -} - //_____________________________________________________________________________ Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent) { @@ -419,80 +391,55 @@ Bool_t AliTRDclusterizer::WriteClusters(Int_t det) AliError(Form("Unexpected detector index %d.\n",det)); return kFALSE; } + Int_t nRecPoints = RecPoints()->GetEntriesFast(); + if(!nRecPoints) return kTRUE; TObjArray *ioArray = new TObjArray(400); TBranch *branch = fClusterTree->GetBranch("TRDcluster"); if (!branch) { - branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0); + fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0); } else branch->SetAddress(&ioArray); - Int_t nRecPoints = RecPoints()->GetEntriesFast(); + AliTRDcluster *c(NULL); if(det >= 0){ for (Int_t i = 0; i < nRecPoints; i++) { - AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i); + if(!(c = (AliTRDcluster *) RecPoints()->UncheckedAt(i))) continue; if(det != c->GetDetector()) continue; ioArray->AddLast(c); } fClusterTree->Fill(); + ioArray->Clear(); } else { - - Int_t detOld = -1; - for (Int_t i = 0; i < nRecPoints; i++) { - AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i); + if(!(c = (AliTRDcluster*)RecPoints()->UncheckedAt(0))){ + AliError("Missing first cluster."); + delete ioArray; + return kFALSE; + } + Int_t detOld(c->GetDetector()), nw(0); + ioArray->AddLast(c); + for (Int_t i(1); iUncheckedAt(i))) continue; if(c->GetDetector() != detOld){ + nw += ioArray->GetEntriesFast(); + // fill & clear previously detector set of clusters fClusterTree->Fill(); ioArray->Clear(); detOld = c->GetDetector(); } ioArray->AddLast(c); } + if(ioArray->GetEntriesFast()){ + nw += ioArray->GetEntriesFast(); + // fill & clear last detector set of clusters (if any) + fClusterTree->Fill(); + ioArray->Clear(); + } + AliDebug(2, Form("Clusters FOUND[%d] WRITTEN[%d] STATUS[%s]", nRecPoints, nw, nw==nRecPoints?"OK":"FAILED")); + if(nw!=nRecPoints) AliWarning(Form("Clusters FOUND[%d] WRITTEN[%d]", nRecPoints, nw)); } delete ioArray; return kTRUE; - -} - -//_____________________________________________________________________________ -Bool_t AliTRDclusterizer::WriteTracklets(Int_t det) -{ - // - // Write the raw data tracklets into seperate file - // - - UInt_t **leaves = new UInt_t *[2]; - for (Int_t i=0; i<2 ;i++){ - leaves[i] = new UInt_t[258]; - leaves[i][0] = det; // det - leaves[i][1] = i; // side - memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256); - } - - if (!fTrackletTree){ - AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets"); - dl->MakeTree(); - fTrackletTree = dl->Tree(); - } - - TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch"); - if (!trkbranch) { - trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i"); - } - - for (Int_t i=0; i<2; i++){ - if (leaves[i][2]>0) { - trkbranch->SetAddress(leaves[i]); - fTrackletTree->Fill(); - } - } - - AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets"); - dl->WriteData("OVERWRITE"); - //dl->Unload(); - delete [] leaves; - - return kTRUE; - } //_____________________________________________________________________________ @@ -543,6 +490,122 @@ Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader) } +Bool_t AliTRDclusterizer::ReadTracklets() +{ + // + // Reads simulated tracklets from the input aliroot file + // + + AliRunLoader *runLoader = AliRunLoader::Instance(); + if (!runLoader) { + AliError("No run loader available"); + return kFALSE; + } + + AliLoader* loader = runLoader->GetLoader("TRDLoader"); + + AliDataLoader *trackletLoader = loader->GetDataLoader("tracklets"); + if (!trackletLoader) { + return kFALSE; + } + trackletLoader->Load(); + + Bool_t loaded = kFALSE; + // look for simulated tracklets + TTree *trackletTree = trackletLoader->Tree(); + + if (trackletTree) { + TBranch *trklbranch = trackletTree->GetBranch("mcmtrklbranch"); + TClonesArray *trklArray = TrackletsArray("AliTRDtrackletMCM"); + if (trklbranch && trklArray) { + AliTRDtrackletMCM *trkl = 0x0; + trklbranch->SetAddress(&trkl); + Int_t nTracklets = trklbranch->GetEntries(); + for (Int_t iTracklet = 0; iTracklet < nTracklets; iTracklet++) { + trklbranch->GetEntry(iTracklet); + new ((*trklArray)[trklArray->GetEntries()]) AliTRDtrackletMCM(*trkl); + } + loaded = kTRUE; + } + } + else { + // if no simulated tracklets found, look for raw tracklets + AliTreeLoader *treeLoader = (AliTreeLoader*) trackletLoader->GetBaseLoader("tracklets-raw"); + trackletTree = treeLoader ? treeLoader->Load(), treeLoader->Tree() : 0x0; + + if (trackletTree) { + TClonesArray *trklArray = TrackletsArray("AliTRDtrackletWord"); + + Int_t hc; + TClonesArray *ar = 0x0; + trackletTree->SetBranchAddress("hc", &hc); + trackletTree->SetBranchAddress("trkl", &ar); + + Int_t nEntries = trackletTree->GetEntries(); + for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { + trackletTree->GetEntry(iEntry); + Int_t nTracklets = ar->GetEntriesFast(); + AliDebug(2, Form("%i tracklets in HC %i", nTracklets, hc)); + for (Int_t iTracklet = 0; iTracklet < nTracklets; iTracklet++) { + AliTRDtrackletWord *trklWord = (AliTRDtrackletWord*) (*ar)[iTracklet]; + new ((*trklArray)[trklArray->GetEntries()]) AliTRDtrackletWord(trklWord->GetTrackletWord(), hc); + } + } + loaded = kTRUE; + } + } + + trackletLoader->UnloadAll(); + trackletLoader->CloseFile(); + + return loaded; +} + +Bool_t AliTRDclusterizer::ReadTracks() +{ + // + // Reads simulated GTU tracks from the input aliroot file + // + + AliRunLoader *runLoader = AliRunLoader::Instance(); + + if (!runLoader) { + AliError("No run loader available"); + return kFALSE; + } + + AliLoader* loader = runLoader->GetLoader("TRDLoader"); + if (!loader) { + return kFALSE; + } + + AliDataLoader *trackLoader = loader->GetDataLoader("gtutracks"); + if (!trackLoader) { + return kFALSE; + } + + Bool_t loaded = kFALSE; + + trackLoader->Load(); + + TTree *trackTree = trackLoader->Tree(); + if (trackTree) { + TClonesArray *trackArray = TracksArray(); + AliTRDtrackGTU *trk = 0x0; + trackTree->SetBranchAddress("TRDtrackGTU", &trk); + for (Int_t iTrack = 0; iTrack < trackTree->GetEntries(); iTrack++) { + trackTree->GetEntry(iTrack); + new ((*trackArray)[trackArray->GetEntries()]) AliESDTrdTrack(*(trk->CreateTrdTrack())); + } + loaded = kTRUE; + } + + trackLoader->UnloadAll(); + trackLoader->CloseFile(); + + return loaded; +} + //_____________________________________________________________________________ Bool_t AliTRDclusterizer::MakeClusters() { @@ -562,21 +625,26 @@ Bool_t AliTRDclusterizer::MakeClusters() // This is to take care of switched off super modules if (!digitsIn->HasData()) continue; digitsIn->Expand(); - // Subtract the ADC baseline - digitsIn->SubtractBaseline(fDigitsManager->GetDigitsParam()->GetADCbaseline()); digitsIn->DeleteNegatives(); // Restore digits array to >=0 values AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i); if (indexes->IsAllocated() == kFALSE){ fDigitsManager->BuildIndexes(i); } - Bool_t fR = kFALSE; + Bool_t fR(kFALSE); if (indexes->HasEntry()){ if (TestBit(kLabels)){ + Int_t nDict(0); for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){ - AliTRDarrayDictionary *tracksIn = 0; //mod + AliTRDarrayDictionary *tracksIn(NULL); //mod tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict); //mod - tracksIn->Expand(); + // This is to take care of data reconstruction + if (!tracksIn->GetDim()) continue; + tracksIn->Expand(); nDict++; + } + if(!nDict){ + AliDebug(1, "MC labels not available. Switch them off."); + SetUseLabels(kFALSE); } } fR = MakeClusters(i); @@ -588,15 +656,16 @@ Bool_t AliTRDclusterizer::MakeClusters() // ResetRecPoints(); //} - // No compress just remove - fDigitsManager->RemoveDigits(i); - fDigitsManager->RemoveDictionaries(i); - fDigitsManager->ClearIndexes(i); + // Clear arrays of this chamber, to prepare for next event + fDigitsManager->ClearArrays(i); } if(fReconstructor->IsWritingClusters()) WriteClusters(-1); - AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast())); + AliInfo(Form("Found :: clusters[%d] tracklets[%d] tracks[%d]", + RecPoints()?RecPoints()->GetEntriesFast():0, + TrackletsArray()?TrackletsArray()->GetEntriesFast():0, + TracksArray()?TracksArray()->GetEntriesFast():0)); return fReturn; @@ -622,57 +691,51 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader) // Create the digits manager if (!fDigitsManager){ + SetBit(knewDM, kTRUE); fDigitsManager = new AliTRDdigitsManager(kTRUE); fDigitsManager->CreateArrays(); } fDigitsManager->SetUseDictionaries(TestBit(kLabels)); - // tracklet container for raw tracklet writing - if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) { - // maximum tracklets for one HC - const Int_t kTrackletChmb=256; - fTrackletContainer = new UInt_t *[2]; - fTrackletContainer[0] = new UInt_t[kTrackletChmb]; - fTrackletContainer[1] = new UInt_t[kTrackletChmb]; - } + if(!fRawStream) + fRawStream = new AliTRDrawStream(rawReader); + else + fRawStream->SetReader(rawReader); - AliTRDrawStreamBase::SetSubtractBaseline(10); - AliTRDrawStreamBase *input = AliTRDrawStreamBase::GetRawStream(rawReader); - if(fReconstructor->IsHLT()) - input->SetSharedPadReadout(kFALSE); + //if(fReconstructor->IsHLT()){ + fRawStream->DisableErrorStorage(); + //} - AliInfo(Form("Stream version: %s", input->IsA()->GetName())); - - Int_t det = 0; - while ((det = input->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){ - Bool_t iclusterBranch = kFALSE; + // register tracklet array for output + fRawStream->SetTrackletArray(TrackletsArray("AliTRDtrackletWord")); + fRawStream->SetTrackArray(TracksArray()); + + UInt_t det = 0; + while ((det = fRawStream->NextChamber(fDigitsManager)) < AliTRDgeometry::kNdet){ if (fDigitsManager->GetIndexes(det)->HasEntry()){ - iclusterBranch = MakeClusters(det); + MakeClusters(det); + fDigitsManager->ClearArrays(det); } - - fDigitsManager->ResetArrays(det); - - if (!fReconstructor->IsWritingTracklets()) continue; - if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det); } - - if (fTrackletContainer){ - delete [] fTrackletContainer[0]; - delete [] fTrackletContainer[1]; - delete [] fTrackletContainer; - fTrackletContainer = NULL; + + for (Int_t iSector = 0; iSector < AliTRDgeometry::kNsector; iSector++) { + fTrgFlags[iSector] = fRawStream->GetTriggerFlags(iSector); } if(fReconstructor->IsWritingClusters()) WriteClusters(-1); - delete fDigitsManager; - fDigitsManager = NULL; - - delete input; - input = NULL; + if(!TestBit(knewDM)){ + delete fDigitsManager; + fDigitsManager = NULL; + delete fRawStream; + fRawStream = NULL; + } - AliInfo(Form("Number of found clusters : %d", fNoOfClusters)); + AliInfo(Form("Found :: clusters[%d] tracklets[%d] tracks[%d]", + RecPoints()?RecPoints()->GetEntriesFast():0, + TrackletsArray()?TrackletsArray()->GetEntriesFast():0, + TracksArray()?TracksArray()->GetEntriesFast():0)); return kTRUE; } @@ -731,23 +794,19 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det) // // Get the digits - // digits should be expanded beforehand! - // digitsIn->Expand(); - fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod + fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod + fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline(det); // This is to take care of switched off super modules if (!fDigits->HasData()) return kFALSE; - // Subtract the ADC baseline - fDigits->SubtractBaseline(fDigitsManager->GetDigitsParam()->GetADCbaseline()); - fIndexes = fDigitsManager->GetIndexes(det); if (fIndexes->IsAllocated() == kFALSE) { AliError("Indexes do not exist!"); return kFALSE; } - AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); + AliTRDcalibDB* const calibration = AliTRDcalibDB::Instance(); if (!calibration) { AliFatal("No AliTRDcalibDB instance available\n"); return kFALSE; @@ -758,10 +817,14 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det) return kFALSE; } - fMaxThresh = fReconstructor->GetRecoParam()->GetClusMaxThresh(); - fSigThresh = fReconstructor->GetRecoParam()->GetClusSigThresh(); - fMinMaxCutSigma = fReconstructor->GetRecoParam()->GetMinMaxCutSigma(); - fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma(); + const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam(); + + fMaxThresh = recoParam->GetClusMaxThresh(); + fMaxThreshTest = (recoParam->GetClusMaxThresh()/2+fBaseline); + fSigThresh = recoParam->GetClusSigThresh(); + fMinMaxCutSigma = recoParam->GetMinMaxCutSigma(); + fMinLeftRightCutSigma = recoParam->GetMinLeftRightCutSigma(); + const Int_t iEveryNTB = recoParam->GetRecEveryNTB(); Int_t istack = fIndexes->GetStack(); fLayer = fIndexes->GetLayer(); @@ -771,10 +834,12 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det) fDet = AliTRDgeometry::GetDetector(fLayer,istack,isector); if (fDet != det) { - AliError("Strange Detector number Missmatch!"); + AliError(Form("Detector number missmatch! Request[%03d] RAW[%03d]", det, fDet)); return kFALSE; } + AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer)); + // TRD space point transformation fTransform->SetDetector(det); @@ -782,12 +847,41 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det) Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector; fVolid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule); - if(fReconstructor->IsProcessingTracklets() && fTrackletContainer) - AddTrackletsToArray(); - fColMax = fDigits->GetNcol(); - //Int_t nRowMax = fDigits->GetNrow(); - fTimeTotal = fDigits->GetNtime(); + fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det); + + // Check consistency between Geometry and raw data + AliTRDpadPlane *pp(fTransform->GetPadPlane()); + Int_t ncols(pp->GetNcols()), nrows(pp->GetNrows()); + if(ncols != fColMax) AliDebug(1, Form("N cols missmatch in Digits for Det[%3d] :: Geom[%3d] RAW[%3d]", fDet, ncols, fColMax)); + if(nrows != fDigits->GetNrow()) AliDebug(1, Form("N rows missmatch in Digits for Det[%3d] :: Geom[%3d] RAW[%3d]", fDet, nrows, fDigits->GetNrow())); + if(ncols != fIndexes->GetNcol()) AliDebug(1, Form("N cols missmatch in Digits for Det[%3d] :: Geom[%3d] RAW[%3d]", fDet, ncols, fIndexes->GetNcol())); + if(nrows != fIndexes->GetNrow()) AliDebug(1, Form("N rows missmatch in Digits for Det[%3d] :: Geom[%3d] RAW[%3d]", fDet, nrows, fIndexes->GetNrow())); + + // Check consistency between OCDB and raw data + Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS(); + if(fReconstructor->IsHLT()){ + if((nTimeOCDB > -1) && (fTimeTotal != nTimeOCDB)){ + AliWarning(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d]), using raw value" + ,fTimeTotal,nTimeOCDB)); + } + }else{ + if(nTimeOCDB == -1){ + AliDebug(1, "Undefined number of timebins in OCDB, using value from raw data."); + if(!(fTimeTotal>0)){ + AliError(Form("Number of timebins in raw data is negative, skipping chamber[%3d]!", fDet)); + return kFALSE; + } + }else if(nTimeOCDB == -2){ + AliError("Mixed number of timebins in OCDB, no reconstruction of TRD data!"); + return kFALSE; + }else if(fTimeTotal != nTimeOCDB){ + AliError(Form("Number of timebins in raw data does not match OCDB value (RAW[%d] OCDB[%d]), skipping chamber[%3d]!" + ,fTimeTotal,nTimeOCDB, fDet)); + return kFALSE; + } + } + AliDebug(1, Form("Using %2d number of timebins for Det[%03d].", fTimeTotal, fDet)); // Detector wise calibration object for the gain factors const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet(); @@ -805,39 +899,51 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det) // Calibration object with the pad status fCalPadStatusROC = calibration->GetPadStatusROC(fDet); + // Calibration object of the online gain + fCalOnlGainROC = 0x0; + if (calibration->HasOnlineFilterGain()) { + fCalOnlGainROC = calibration->GetOnlineGainTableROC(fDet); + } - SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT()); - SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS()); - SetBit(kHLT, fReconstructor->IsHLT()); - firstClusterROC = -1; fClusterROC = 0; + SetBit(kLUT, recoParam->UseLUT()); + SetBit(kGAUS, recoParam->UseGAUS()); + // Apply the gain and the tail cancelation via digital filter - if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation(); + // Use the configuration from the DCS to find out whether online + // tail cancellation was applied + if(!calibration->HasOnlineTailCancellation()){ + // save a copy of raw data + if(TestBit(kRawSignal)){ + if(fDigitsRaw){ + fDigitsRaw->~AliTRDarrayADC(); + new(fDigitsRaw) AliTRDarrayADC(*fDigits); + } else fDigitsRaw = new AliTRDarrayADC(*fDigits); + } + TailCancelation(recoParam); + } MaxStruct curr, last; Int_t nMaximas = 0, nCorrupted = 0; // Here the clusterfining is happening - for(curr.Time = 0; curr.Time < fTimeTotal; curr.Time++){ - while(fIndexes->NextRCIndex(curr.Row, curr.Col)){ - //printf("\nCHECK r[%2d] c[%3d] t[%d]\n", curr.Row, curr.Col, curr.Time); - if(IsMaximum(curr, curr.padStatus, &curr.Signals[0])){ - //printf("\tMAX s[%d %d %d]\n", curr.Signals[0], curr.Signals[1], curr.Signals[2]); - if(last.Row>-1){ - if(curr.Time==last.Time && curr.Row==last.Row && curr.Col==last.Col+2) FivePadCluster(last, curr); + for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){ + while(fIndexes->NextRCIndex(curr.row, curr.col)){ + if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){ + if(last.row>-1){ + if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr); CreateCluster(last); } - last=curr; curr.FivePad=kFALSE; + last=curr; curr.fivePad=kFALSE; } - //printf("\t--- s[%d %d %d]\n", curr.Signals[0], curr.Signals[1], curr.Signals[2]); } } - if(last.Row>-1) CreateCluster(last); + if(last.row>-1) CreateCluster(last); - if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){ + if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){ TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer); (*fDebugStream) << "MakeClusters" << "Detector=" << det @@ -853,54 +959,68 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det) } //_____________________________________________________________________________ -Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals) +Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Float_t *const Signals) { // // Returns true if this row,col,time combination is a maximum. // Gives back the padStatus and the signals of the center pad and the two neighbouring pads. // - Signals[1] = fDigits->GetData(Max.Row, Max.Col, Max.Time); - if(Signals[1] < fMaxThresh) return kFALSE; + Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row); + Float_t onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col) : 1; + Signals[1] = (fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) /(onlcf * gain) + 0.5f; + if(Signals[1] <= fMaxThresh) return kFALSE; - Float_t noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.Col, Max.Row); - if (Signals[1] < noiseMiddleThresh) return kFALSE; + if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE; - if (Max.Col + 1 >= fColMax || Max.Col < 1) return kFALSE; + Float_t noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row); + if (Signals[1] <= noiseMiddleThresh) return kFALSE; - UChar_t status[3]={ - fCalPadStatusROC->GetStatus(Max.Col-1, Max.Row) - ,fCalPadStatusROC->GetStatus(Max.Col, Max.Row) - ,fCalPadStatusROC->GetStatus(Max.Col+1, Max.Row) + Char_t status[3]={ + fCalPadStatusROC->GetStatus(Max.col-1, Max.row) + ,fCalPadStatusROC->GetStatus(Max.col, Max.row) + ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row) }; - Signals[0] = fDigits->GetData(Max.Row, Max.Col-1, Max.Time); - Signals[2] = fDigits->GetData(Max.Row, Max.Col+1, Max.Time); + Short_t signal(0); + if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){ + gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row); + onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col-1) : 1; + Signals[0] = (signal - fBaseline) /( onlcf * gain) + 0.5f; + } else Signals[0] = 0.; + if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){ + gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row); + onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col+1) : 1; + Signals[2] = (signal - fBaseline) /( onlcf * gain) + 0.5f; + } else Signals[2] = 0.; if(!(status[0] | status[1] | status[2])) {//all pads are good if ((Signals[2] <= Signals[1]) && (Signals[0] < Signals[1])) { - if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) { - Float_t noiseSumThresh = fMinLeftRightCutSigma - * fCalNoiseDetValue - * fCalNoiseROC->GetValue(Max.Col, Max.Row); - if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE; + if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) { + if(Signals[0]<0) Signals[0]=0.; + if(Signals[2]<0) Signals[2]=0.; + Float_t noiseSumThresh = fMinLeftRightCutSigma * fCalNoiseDetValue + * fCalNoiseROC->GetValue(Max.col, Max.row); + if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE; padStatus = 0; return kTRUE; } } } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad - if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) { + if(Signals[0]<0)Signals[0]=0; + if(Signals[2]<0)Signals[2]=0; + if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) { Signals[2]=0; SetPadStatus(status[2], padStatus); return kTRUE; } - else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) { + else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) { Signals[0]=0; SetPadStatus(status[0], padStatus); return kTRUE; } - else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) { - Signals[1]=TMath::Nint(fMaxThresh); + else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) { + Signals[1] = fMaxThresh; SetPadStatus(status[1], padStatus); return kTRUE; } @@ -915,27 +1035,34 @@ Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &Neighbou // Look for 5 pad cluster with minimum in the middle // Gives back the ratio // - if (ThisMax.Col >= fColMax - 3) return kFALSE; - if (ThisMax.Col < fColMax - 5){ - if (fDigits->GetData(ThisMax.Row, ThisMax.Col+4, ThisMax.Time) >= fSigThresh) + + if (ThisMax.col >= fColMax - 3) return kFALSE; + Float_t gain; + Float_t onlcf; + if (ThisMax.col < fColMax - 5){ + gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row); + onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(ThisMax.row,ThisMax.col+4) : 1; + if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain * onlcf) return kFALSE; } - if (ThisMax.Col > 1) { - if (fDigits->GetData(ThisMax.Row, ThisMax.Col-2, ThisMax.Time) >= fSigThresh) + if (ThisMax.col > 1) { + gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row); + onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(ThisMax.row,ThisMax.col-2) : 1; + if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain * onlcf) return kFALSE; } const Float_t kEpsilon = 0.01; - Double_t padSignal[5] = {ThisMax.Signals[0], ThisMax.Signals[1], ThisMax.Signals[2], - NeighbourMax.Signals[1], NeighbourMax.Signals[2]}; + Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2], + NeighbourMax.signals[1], NeighbourMax.signals[2]}; // Unfold the two maxima and set the signal on // the overlapping pad to the ratio Float_t ratio = Unfold(kEpsilon,fLayer,padSignal); - ThisMax.Signals[2] = TMath::Nint(ThisMax.Signals[2]*ratio); - NeighbourMax.Signals[0] = TMath::Nint(NeighbourMax.Signals[0]*(1-ratio)); - ThisMax.FivePad=kTRUE; - NeighbourMax.FivePad=kTRUE; + ThisMax.signals[2] = ThisMax.signals[2]*ratio + 0.5f; + NeighbourMax.signals[0] = NeighbourMax.signals[0]*(1-ratio) + 0.5f; + ThisMax.fivePad=kTRUE; + NeighbourMax.fivePad=kTRUE; return kTRUE; } @@ -944,39 +1071,59 @@ Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &Neighbou void AliTRDclusterizer::CreateCluster(const MaxStruct &Max) { // - // Creates a cluster at the given position and saves it in fRecPoints + // Creates a cluster at the given position and saves it in RecPoints // Int_t nPadCount = 1; - Short_t signals[7] = { 0, 0, Max.Signals[0], Max.Signals[1], Max.Signals[2], 0, 0 }; - if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount); + Short_t signals[7] = { 0, 0, (Short_t)Max.signals[0], (Short_t)Max.signals[1], (Short_t)Max.signals[2], 0, 0 }; + if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount); - AliTRDcluster cluster(fDet, ((UChar_t) Max.Col), ((UChar_t) Max.Row), ((UChar_t) Max.Time), signals, fVolid); + AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid); cluster.SetNPads(nPadCount); + cluster.SetQ(Max.signals[0]+Max.signals[1]+Max.signals[2]); if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT); else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS); else cluster.SetRPhiMethod(AliTRDcluster::kCOG); - cluster.SetFivePad(Max.FivePad); + cluster.SetFivePad(Max.fivePad); // set pads status for the cluster UChar_t maskPosition = GetCorruption(Max.padStatus); if (maskPosition) { cluster.SetPadMaskedPosition(maskPosition); cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus)); } + cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection()); // Transform the local cluster coordinates into calibrated // space point positions defined in the local tracking system. // Here the calibration for T0, Vdrift and ExB is applied as well. - if(!fTransform->Transform(&cluster)) return; + if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return; + // Store raw signals in cluster. This MUST be called after position reconstruction ! + // Xianguo Lu and Alex Bercuci 19.03.2012 + if(TestBit(kRawSignal) && fDigitsRaw){ + Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration + Short_t rawSignal[7] = {0}; + for(Int_t ipad(Max.col-3), iRawId(0); ipad<=Max.col+3; ipad++, iRawId++){ + if(ipad<0 || ipad>=fColMax) continue; + if(!fCalOnlGainROC){ + rawSignal[iRawId] = fDigitsRaw->GetData(Max.row, ipad, Max.time); + continue; + } + // Deconvolute online gain calibration when available + // Alex Bercuci 27.04.2012 + tmp = (fDigitsRaw->GetData(Max.row, ipad, Max.time) - fBaseline)/fCalOnlGainROC->GetGainCorrectionFactor(Max.row, ipad) + 0.5f; + rawSignal[iRawId] = (Short_t)TMath::Min(tmp, kMaxShortVal); + } + cluster.SetSignals(rawSignal, kTRUE); + } // Temporarily store the Max.Row, column and time bin of the center pad // Used to later on assign the track indices - cluster.SetLabel(Max.Row, 0); - cluster.SetLabel(Max.Col, 1); - cluster.SetLabel(Max.Time,2); + cluster.SetLabel(Max.row, 0); + cluster.SetLabel(Max.col, 1); + cluster.SetLabel(Max.time,2); //needed for HLT reconstruction - AddClusterToArray(&cluster); + AddClusterToArray(&cluster); // Store the index of the first cluster in the current ROC if (firstClusterROC < 0) firstClusterROC = fNoOfClusters; @@ -988,37 +1135,43 @@ void AliTRDclusterizer::CreateCluster(const MaxStruct &Max) //_____________________________________________________________________________ void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount) { +// Calculate number of pads/cluster and +// ADC signals at position 0, 1, 5 and 6 + + Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration + Float_t gain(1.); Float_t onlcf(1.); Short_t signal(0); + // Store the amplitudes of the pads in the cluster for later analysis + // and check whether one of these pads is masked in the database + signals[3]=Max.signals[1]; + Int_t ipad(1), jpad(0); // Look to the right - Int_t ii = 1; - while (fDigits->GetData(Max.Row, Max.Col-ii, Max.Time) >= fSigThresh) { + while((jpad = Max.col-ipad)){ + if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit ! + gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row); + onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,jpad) : 1; + tmp = (signal - fBaseline) / (onlcf * gain) + 0.5f; + signal = (Short_t)TMath::Min(tmp, kMaxShortVal); + if(signalGetData(Max.Row, Max.Col+ii, Max.Time) >= fSigThresh) { + while((jpad = Max.col+ipad)GetData(Max.row, jpad, Max.time))) break; // empty digit ! + gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row); + onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,jpad) : 1; + tmp = (signal - fBaseline) / (onlcf * gain) + 0.5f; + signal = (Short_t)TMath::Min(tmp, kMaxShortVal); + if(signal= fColMax) break; + if(ipad<=3) signals[3 + ipad] = signal; + ipad++; } - // Store the amplitudes of the pads in the cluster for later analysis - // and check whether one of these pads is masked in the database - signals[2]=Max.Signals[0]; - signals[3]=Max.Signals[1]; - signals[4]=Max.Signals[2]; - for(Int_t i = 0; i<2; i++) - { - if(Max.Col+i >= 3) - signals[i] = fDigits->GetData(Max.Row, Max.Col-3+i, Max.Time); - if(Max.Col+3-i < fColMax) - signals[6-i] = fDigits->GetData(Max.Row, Max.Col+3-i, Max.Time); - } - /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) { - if ((jPad >= 0) && (jPad < fColMax)) - signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time)); - }*/ + AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]." + , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount)); } //_____________________________________________________________________________ @@ -1033,27 +1186,6 @@ void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster) new((*RecPoints())[n]) AliTRDcluster(*cluster); } -//_____________________________________________________________________________ -void AliTRDclusterizer::AddTrackletsToArray() -{ - // - // Add the online tracklets of this chamber to the array - // - - UInt_t* trackletword; - for(Int_t side=0; side<2; side++) - { - Int_t trkl=0; - trackletword=fTrackletContainer[side]; - while(trackletword[trkl]>0){ - Int_t n = TrackletsArray()->GetEntriesFast(); - AliTRDtrackletWord tmp(trackletword[trkl]); - new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid); - trkl++; - } - } -} - //_____________________________________________________________________________ Bool_t AliTRDclusterizer::AddLabels() { @@ -1077,7 +1209,7 @@ Bool_t AliTRDclusterizer::AddLabels() // Loop through the dictionary arrays one-by-one // to keep memory consumption low - AliTRDarrayDictionary *tracksIn = 0; //mod + AliTRDarrayDictionary *tracksIn(NULL); //mod for (Int_t iDict = 0; iDict < kNdict; iDict++) { // tracksIn should be expanded beforehand! @@ -1171,8 +1303,11 @@ Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *con irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal); // Calculate new overlapping ratio - ratio = TMath::Min((Double_t) 1.0 - ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0])); + // Coverity + if (irc != 0) { + ratio = TMath::Min((Double_t) 1.0 + ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0])); + } } @@ -1181,89 +1316,122 @@ Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *con } //_____________________________________________________________________________ -void AliTRDclusterizer::TailCancelation() +void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam) { // - // Applies the tail cancelation and gain factors: - // Transform fDigits to fDigits + // Applies the tail cancelation // + Int_t nexp = recoParam->GetTCnexp(); + if(!nexp) return; + Int_t iRow = 0; Int_t iCol = 0; Int_t iTime = 0; - Double_t *inADC = new Double_t[fTimeTotal]; // ADC data before tail cancellation - Double_t *outADC = new Double_t[fTimeTotal]; // ADC data after tail cancellation - - fIndexes->ResetCounters(); TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer); + Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming(); while(fIndexes->NextRCIndex(iRow, iCol)) { - Float_t fCalGainFactorROCValue = fCalGainFactorROC->GetValue(iCol,iRow); - Double_t gain = fCalGainFactorDetValue - * fCalGainFactorROCValue; - - Bool_t corrupted = kFALSE; - for (iTime = 0; iTime < fTimeTotal; iTime++) - { - // Apply gain gain factor - inADC[iTime] = fDigits->GetData(iRow,iCol,iTime); - if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE; - inADC[iTime] /= gain; - outADC[iTime] = inADC[iTime]; - if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming()){ - (*fDebugStream) << "TailCancellation" - << "col=" << iCol - << "row=" << iRow - << "time=" << iTime - << "inADC=" << inADC[iTime] - << "gain=" << gain - << "outADC=" << outADC[iTime] - << "corrupted=" << corrupted - << "\n"; - } - } - if (!corrupted) - { - // Apply the tail cancelation via the digital filter - // (only for non-coorupted pads) - DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp()); - } + // if corrupted then don't make the tail cancallation + if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue; + + if(debugStreaming){ + for (iTime = 0; iTime < fTimeTotal; iTime++) + (*fDebugStream) << "TailCancellation" + << "col=" << iCol + << "row=" << iRow + << "\n"; + } + + // Apply the tail cancelation via the digital filter + //DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp); + ApplyTCTM(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp); + } // while irow icol - for(iTime = 0; iTime < fTimeTotal; iTime++)//while (fIndexes->NextTbinIndex(iTime)) - { - // Store the amplitude of the digit if above threshold - if (outADC[iTime] > 0) - fDigits->SetData(iRow,iCol,iTime,TMath::Nint(outADC[iTime])); - else - fDigits->SetData(iRow,iCol,iTime,0); - } // while itime + return; - } // while irow icol +} - delete [] inADC; - delete [] outADC; - return; +//_____________________________________________________________________________ +void AliTRDclusterizer::ApplyTCTM(Short_t *const arr, const Int_t nTime, const Int_t nexp) +{ + // + // Steer tail cancellation + // + + + switch(nexp) { + case 1: + case 2: + DeConvExp(arr,nTime,nexp); + break; + case -1: + ConvExp(arr,nTime); + break; + case -2: + DeConvExp(arr,nTime,1); + ConvExp(arr,nTime); + break; + default: + break; + } +} + + +//_____________________________________________________________________________ +void AliTRDclusterizer::ConvExp(Short_t *const arr, const Int_t nTime) +{ + // + // Tail maker + // + // Initialization (coefficient = alpha, rates = lambda) + Float_t slope = 1.0; + Float_t coeff = 0.5; + Float_t rate; + + Double_t par[4]; + fReconstructor->GetRecoParam()->GetTCParams(par); + slope = par[1]; + coeff = par[3]; + + Double_t dt = 0.1; + + rate = TMath::Exp(-dt/(slope)); + + Float_t reminder = .0; + Float_t correction = 0.0; + Float_t result = 0.0; + + for (int i = nTime-1; i >= 0; i--) { + + result = arr[i] + correction - fBaseline; // No rescaling + arr[i] = (Short_t)(result + fBaseline + 0.5f); + + correction = 0.0; + + correction += reminder = rate * (reminder + coeff * result); + } } + //_____________________________________________________________________________ -void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target - ,const Int_t n, const Int_t nexp) +void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp) { // // Tail cancellation by deconvolution for PASA v4 TRF // - Double_t rates[2]; - Double_t coefficients[2]; + Float_t rates[2]; + Float_t coefficients[2]; // Initialization (coefficient = alpha, rates = lambda) - Double_t r1 = 1.0; - Double_t r2 = 1.0; - Double_t c1 = 0.5; - Double_t c2 = 0.5; + Float_t r1 = 1.0; + Float_t r2 = 1.0; + Float_t c1 = 0.5; + Float_t c2 = 0.5; if (nexp == 1) { // 1 Exponentials r1 = 1.156; @@ -1286,48 +1454,21 @@ void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const Double_t dt = 0.1; rates[0] = TMath::Exp(-dt/(r1)); - rates[1] = TMath::Exp(-dt/(r2)); + rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2)); - Int_t i = 0; - Int_t k = 0; - - Double_t reminder[2]; - Double_t correction = 0.0; - Double_t result = 0.0; + Float_t reminder[2] = { .0, .0 }; + Float_t correction = 0.0; + Float_t result = 0.0; - // Attention: computation order is important - for (k = 0; k < nexp; k++) { - reminder[k] = 0.0; - } - - for (i = 0; i < n; i++) { + for (int i = 0; i < nTime; i++) { - result = (source[i] - correction); // No rescaling - target[i] = result; - - for (k = 0; k < nexp; k++) { - reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result); - } + result = arr[i] - correction - fBaseline; // No rescaling + arr[i] = (Short_t)(result + fBaseline + 0.5f); correction = 0.0; - for (k = 0; k < nexp; k++) { - correction += reminder[k]; + for (int k = 0; k < 2; k++) { + correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result); } - - } - -} - -//_____________________________________________________________________________ -void AliTRDclusterizer::ResetRecPoints() -{ - // - // Resets the list of rec points - // - - if (fRecPoints) { - fRecPoints->Delete(); - delete fRecPoints; } } @@ -1338,33 +1479,28 @@ TClonesArray *AliTRDclusterizer::RecPoints() // Returns the list of rec points // - if (!fRecPoints) { - if(!(fRecPoints = AliTRDReconstructor::GetClusters())){ - // determine number of clusters which has to be allocated - Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters(); - - fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters)); - } - //SetClustersOwner(kTRUE); - AliTRDReconstructor::SetClusters(0x0); - } + fRecPoints = AliTRDReconstructor::GetClusters(); + if (!fRecPoints) AliError("Missing cluster array"); return fRecPoints; - } //_____________________________________________________________________________ -TClonesArray *AliTRDclusterizer::TrackletsArray() +TClonesArray *AliTRDclusterizer::TrackletsArray(const TString &trkltype) { // - // Returns the list of rec points + // Returns the array of on-line tracklets // - - if (!fTracklets && fReconstructor->IsProcessingTracklets()) { - fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC); - //SetClustersOwner(kTRUE); - //AliTRDReconstructor::SetTracklets(0x0); - } + fTracklets = AliTRDReconstructor::GetTracklets(trkltype.Data()); + if (!fTracklets) AliError("Missing online tracklets array"); return fTracklets; - } +//_____________________________________________________________________________ +TClonesArray* AliTRDclusterizer::TracksArray() +{ + // return array of GTU tracks (create TClonesArray if necessary) + + fTracks = AliTRDReconstructor::GetTracks(); + if (!fTracks) AliError("Missing online tracks array"); + return fTracks; +}