From 12f6dd57adc26c26213ca4b80b3821a06ae3fc67 Mon Sep 17 00:00:00 2001 From: policheh Date: Tue, 20 Oct 2009 00:30:49 +0000 Subject: [PATCH] DA for calibrating energy by pi0 and MIP and related classes (Hisayuki Torii). --- PHOS/AliPHOSDATreeCluster.cxx | 224 +++++++++++++++++++++ PHOS/AliPHOSDATreeCluster.h | 66 +++++++ PHOS/AliPHOSDATreeDigit.cxx | 38 ++++ PHOS/AliPHOSDATreeDigit.h | 49 +++++ PHOS/AliPHOSDATreeEvent.cxx | 203 +++++++++++++++++++ PHOS/AliPHOSDATreeEvent.h | 63 ++++++ PHOS/AliPHOSDApi0mip.cxx | 321 +++++++++++++++++++++++++++++++ PHOS/AliPHOSDApi0mip.h | 59 ++++++ PHOS/AliPHOSPreprocessorPHYS.cxx | 110 +++++++++++ PHOS/AliPHOSPreprocessorPHYS.h | 28 +++ PHOS/PHOSPHYSda.cxx | 193 +++++++++++++++++++ 11 files changed, 1354 insertions(+) create mode 100644 PHOS/AliPHOSDATreeCluster.cxx create mode 100644 PHOS/AliPHOSDATreeCluster.h create mode 100644 PHOS/AliPHOSDATreeDigit.cxx create mode 100644 PHOS/AliPHOSDATreeDigit.h create mode 100644 PHOS/AliPHOSDATreeEvent.cxx create mode 100644 PHOS/AliPHOSDATreeEvent.h create mode 100644 PHOS/AliPHOSDApi0mip.cxx create mode 100644 PHOS/AliPHOSDApi0mip.h create mode 100644 PHOS/AliPHOSPreprocessorPHYS.cxx create mode 100644 PHOS/AliPHOSPreprocessorPHYS.h create mode 100644 PHOS/PHOSPHYSda.cxx diff --git a/PHOS/AliPHOSDATreeCluster.cxx b/PHOS/AliPHOSDATreeCluster.cxx new file mode 100644 index 00000000000..ad72d9b5939 --- /dev/null +++ b/PHOS/AliPHOSDATreeCluster.cxx @@ -0,0 +1,224 @@ +// -- +// -- +// Implementation for TTree output in PHOS DA +// for calibrating energy by pi0 and MIP. +// -- +// -- Author: Hisayuki Torii (Hiroshima Univ.) +// -- + +#include +#include +#include +#include "AliPHOSDATreeDigit.h" +#include "AliPHOSDATreeCluster.h" +ClassImp(AliPHOSDATreeCluster) +//------------------------------------------------------------------------ +AliPHOSDATreeCluster::AliPHOSDATreeCluster(const AliPHOSDATreeCluster& cluster): +//fEnergy(cluster.fEnergy),fX(cluster.fX),fY(cluster.fY),fZ(cluster.fZ),fNDigits(cluster.fNDigits){ +fEnergy(cluster.fEnergy),fRow(cluster.fRow),fCol(cluster.fCol),fNDigits(cluster.fNDigits),fDigits(0){ + // Copy Constructor + + if( fNDigits > 0 ){ + fDigits = new AliPHOSDATreeDigit[fNDigits]; + int ndigits = fNDigits; + while( ndigits-- ){ + fDigits[ndigits] = cluster.fDigits[ndigits]; + } + } else { + fDigits = 0; + } +} +//------------------------------------------------------------------------ +AliPHOSDATreeCluster& AliPHOSDATreeCluster::operator=(const AliPHOSDATreeCluster& cluster){ + // Copy Operator + + if( fNDigits> 0 ) delete[] fDigits; + fEnergy = cluster.fEnergy; + fNDigits = cluster.fNDigits; + fRow = cluster.fRow; + fCol = cluster.fCol; + //fX = cluster.fX; + //fY = cluster.fY; + //fZ = cluster.fZ; + if( fNDigits > 0 ){ + fDigits = new AliPHOSDATreeDigit[fNDigits]; + int ndigits = fNDigits; + while( ndigits-- ){ + fDigits[ndigits] = cluster.fDigits[ndigits]; + } + } else { + fDigits = 0; + } + return *this; +} +//------------------------------------------------------------------------ +void AliPHOSDATreeCluster::Print(char* opt){ + // Print out + std::cout<<" AliPHOSDATreeCluster:: Energy="<["<["< 0 ) delete[] fDigits; + fNDigits = 0; +} +//------------------------------------------------------------------------ +bool AliPHOSDATreeCluster::Append(AliPHOSDATreeDigit& digit){ + // Add digit information and sum all energy + // + if(! digit.IsValid() ){ + std::cout<<" AliPHOSDATreeCluster::Append():: Error!! Digit is not valid.."<0 ) delete[] fDigits; + fNDigits++; + fDigits = newfDigits; + fEnergy += digit.GetEnergy(); + return true; +} +//------------------------------------------------------------------------ +bool AliPHOSDATreeCluster::Append(AliPHOSDATreeCluster& cluster){ + // Add another cluster information and sum all energy + // + AliPHOSDATreeDigit* newfDigits = new AliPHOSDATreeDigit[fNDigits+cluster.fNDigits]; + int ndigits1 = fNDigits; + int ndigits2 = cluster.fNDigits; + int ndigitsall = ndigits1 + ndigits2; + while( ndigitsall-- ){ + //std::cout<<" ------ ndigits1:"<0) delete[] fDigits; + fDigits = newfDigits; + fNDigits += cluster.fNDigits; + fEnergy += cluster.GetEnergy(); + return true; + +} +//------------------------------------------------------------------------ +bool AliPHOSDATreeCluster::IsNeighbor(const AliPHOSDATreeDigit& digit) const{ + // Check wether the given digit is neighboring to this cluster. + // Return true if yes. + + bool status = false; + int ndigits = fNDigits; + while( ndigits-- && !status ){ + status = digit.IsNeighbor(fDigits[ndigits]); + } + return status; +} +//------------------------------------------------------------------------ +bool AliPHOSDATreeCluster::IsNeighbor(const AliPHOSDATreeCluster& cluster) const{ + // Check wether the given cluster is neighboring to this cluster. + // Return true if yes. + + bool status = false; + int ndigits = fNDigits; + while( ndigits-- && !status ){ + status = cluster.IsNeighbor(fDigits[ndigits]); + } + return status; +} +//------------------------------------------------------------------------ +bool AliPHOSDATreeCluster::CalculateProperty(){ + // Calculate the hit position + // (calculation of dispersion is not valid) + + fCol = 0; + fRow = 0; + float totweight = 0; + float weight; + int ndigits = fNDigits; + while( ndigits-- ){ + weight = log(fDigits[ndigits].GetEnergy()/fEnergy) + 4.5; //4.5 is for PHOS + //std::cout<<" AliPHOSDATreeCluster::CalculateProperty() DEBUG: ndigits="<0) delete[] fDigits; }; + AliPHOSDATreeCluster(const AliPHOSDATreeCluster& cluster); + AliPHOSDATreeCluster& operator=(const AliPHOSDATreeCluster& cluster); + void Set(float energy,float row,float col){fEnergy=energy; fRow=row; fCol=col; }; + void SetEnergy(float energy){fEnergy=energy;}; + float GetEnergy() const{ return fEnergy; }; + float GetRow() const{ return fRow; }; + float GetCol() const{ return fCol; }; + bool CalculateProperty(); + int GetNDigits() const{ return fNDigits; }; + AliPHOSDATreeDigit& GetDigit(int ndigit){ + if( ndigit >= 0 && ndigit < fNDigits ) return fDigits[ndigit]; + else std::cout<<" AliPHOSDATreeCluster::GetDigit("< "<= 0 ) return fDigits[0]; + else std::cout<<" AliPHOSDATreeCluster::GetMaxDigit()::Warning No digit information."< +#include +#include +#include "AliPHOSDATreeDigit.h" + +ClassImp(AliPHOSDATreeDigit) +//------------------------------------------------------------------------ +bool AliPHOSDATreeDigit::IsNeighbor(const AliPHOSDATreeDigit& digit) const{ + // Check wether the given digit is neighboring to this digit. + // Return true if yes. + if( fabs(this->GetRow() - digit.GetRow()) < 2 && fabs(this->GetCol() - digit.GetCol()) < 2 ){ + return true; + } + return false; +} +//------------------------------------------------------------------------ +void AliPHOSDATreeDigit::Print(char* opt){ + // Print out + std::cout<<" AliPHOSDATreeDigit:: "<=0 && fAbsId < 17920 ) return true; return false; }; + + private: + + float fEnergy; // Energy in GeV + short fAbsId; // col + row * 56 [+mod*3584 in some cases] + + ClassDef(AliPHOSDATreeDigit,1) // Simple Digit Structure for PHOS DA +}; + +#endif diff --git a/PHOS/AliPHOSDATreeEvent.cxx b/PHOS/AliPHOSDATreeEvent.cxx new file mode 100644 index 00000000000..3ef52a8612f --- /dev/null +++ b/PHOS/AliPHOSDATreeEvent.cxx @@ -0,0 +1,203 @@ +// -- +// -- +// Implementation for TTree output in PHOS DA +// for calibrating energy by pi0 and MIP. +// -- +// -- Author: Hisayuki Torii (Hiroshima Univ.) +// -- + +#include +#include "AliPHOSDATreeEvent.h" +ClassImp(AliPHOSDATreeEvent) +//------------------------------------------------------------------------ +AliPHOSDATreeEvent::AliPHOSDATreeEvent(const AliPHOSDATreeEvent& evt){ + // Copy Constructor + + fNDigits = evt.fNDigits; + if( fNDigits > 0 ){ + fDigits = new AliPHOSDATreeDigit[fNDigits]; + int ndigits = fNDigits; + while( ndigits-- ){ + fDigits[ndigits] = evt.fDigits[ndigits]; + } + } else { + fDigits = 0; + } + // + fNClusters = evt.fNClusters; + if( fNClusters > 0 ){ + fClusters = new AliPHOSDATreeCluster[fNClusters]; + int nclusters = fNClusters; + while( nclusters-- ){ + fClusters[nclusters] = evt.fClusters[nclusters]; + } + } else { + fClusters = 0; + } +} +//------------------------------------------------------------------------ +AliPHOSDATreeEvent& AliPHOSDATreeEvent::operator=(const AliPHOSDATreeEvent& evt){ + // Copy Operator + + if( fNDigits > 0 ) delete[] fDigits; + fNDigits = evt.fNDigits; + if( fNDigits > 0 ){ + fDigits = new AliPHOSDATreeDigit[fNDigits]; + int ndigits = fNDigits; + while( ndigits-- ){ + fDigits[ndigits] = evt.fDigits[ndigits]; + } + } else { + fDigits = 0; + } + // + if( fNClusters > 0 ) delete[] fClusters; + fNClusters = evt.fNClusters; + if( fNClusters > 0 ){ + fClusters = new AliPHOSDATreeCluster[fNClusters]; + int nclusters = fNClusters; + while( nclusters-- ){ + fClusters[nclusters] = evt.fClusters[nclusters]; + } + } else { + fClusters = 0; + } + return *this; +} +//------------------------------------------------------------------------ +bool AliPHOSDATreeEvent::Fill(float energy,int row,int col){ + // Fill new digit information + + AliPHOSDATreeDigit digit(energy,row,col); + return Fill(digit); +} +//------------------------------------------------------------------------ +bool AliPHOSDATreeEvent::Fill(AliPHOSDATreeDigit& digit){ + // Fill new digit information + + AliPHOSDATreeDigit* newDigits = new AliPHOSDATreeDigit[fNDigits+1]; + int ndigits = fNDigits; + newDigits[ndigits] = digit; + while( ndigits-- ){ + newDigits[ndigits] = fDigits[ndigits]; + } + if( fNDigits>0 ) delete[] fDigits; + fDigits = newDigits; + fNDigits++; + return true; +} +//------------------------------------------------------------------------ +bool AliPHOSDATreeEvent::Clusterize(AliPHOSDATreeDigit& digit){ + //Input digit information into clustering + + bool status = false; + int nclusters = fNClusters; + while( nclusters-- && !status ){ + if( fClusters[nclusters].IsNeighbor(digit) ){ + status = fClusters[nclusters].Append(digit); + } + } + //std::cout<<" DEBUGDEBUG:: status = "< "<0 ) delete[] fClusters; + fClusters = newClusters; + } else { + if( fNClusters>0 ) delete[] fClusters; + fClusters = 0; + fNClusters = 0; + } + return true; +} +//------------------------------------------------------------------------ +void AliPHOSDATreeEvent::Print(char* opt){ + // Print Out + + char* when = ctime(&fTime); + std::cout<<" AliPHOSDATreeEvent:: fNDigits = "<["<["<["<["< +#include +#include +#include "AliPHOSDATreeCluster.h" + +class AliPHOSDATreeEvent{ + + friend std::ostream& operator<<(std::ostream& out,const AliPHOSDATreeEvent& event); + + public: + + AliPHOSDATreeEvent(): fTime(0), fNDigits(0), fDigits(0), fNClusters(0), fClusters(0){/* */}; + AliPHOSDATreeEvent(const AliPHOSDATreeEvent& evt); + AliPHOSDATreeEvent& operator=(const AliPHOSDATreeEvent& evt); + virtual ~AliPHOSDATreeEvent(){ + if( fNDigits > 0 ) delete[] fDigits; + if( fNClusters > 0 ) delete[] fClusters; + }; + time_t GetTime() const{return fTime;}; + void SetTime(time_t time){fTime=time;}; + int GetNDigits() const{ return fNDigits; }; + int GetNClusters() const{ return fNClusters; }; + AliPHOSDATreeCluster& GetCluster(int nclusters){ + if( nclusters < fNClusters ) return fClusters[nclusters]; + }; + AliPHOSDATreeDigit& GetDigit(int ndigits){ + if( ndigits < fNDigits ) return fDigits[ndigits]; + }; + bool Fill(float fenergy,int row,int col); + bool Fill(AliPHOSDATreeDigit& digit); + bool ExecuteClustering(); + void Reset(){ + if( fNDigits > 0 ) delete[] fDigits; + if( fNClusters > 0 ) delete[] fClusters; + fTime = 0; + fNDigits = 0; + fNClusters = 0; + }; + void Print(char* opt=""); + + private: + bool Clusterize(AliPHOSDATreeDigit& digit); + + time_t fTime; // Time information + int fNDigits; // Number of digits in event + AliPHOSDATreeDigit* fDigits; //[fNDigits] + int fNClusters; // Number of clusters in event + AliPHOSDATreeCluster* fClusters; //[fNClusters] + + ClassDef(AliPHOSDATreeEvent,1) // Simple Event Structure for PHOS DA +}; +#endif + diff --git a/PHOS/AliPHOSDApi0mip.cxx b/PHOS/AliPHOSDApi0mip.cxx new file mode 100644 index 00000000000..d208950935e --- /dev/null +++ b/PHOS/AliPHOSDApi0mip.cxx @@ -0,0 +1,321 @@ +// -- +// -- +// Implementation for TTree output in PHOS DA +// for calibrating energy by pi0 and MIP. +// -- +// -- Author: Hisayuki Torii (Hiroshima Univ.) +// -- + +#include +#include +#include +#include "AliPHOSDApi0mip.h" +ClassImp(AliPHOSDApi0mip) +//---------------------------------------------------------------- +AliPHOSDApi0mip::AliPHOSDApi0mip(int module,int iterid,char* fopt) : + TNamed(), fCreateTree(false), fCreateHist(false), fMod(0), fIterId(0), + fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0), + fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0), + fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) { + // Constructor + + char hname[1024], htitle[1024]; + sprintf(hname,"AliPHOSDApi0mip_mod%d_iter%d",module,iterid); + SetName(hname); + sprintf(htitle,"PHOS MIP/pi0 Calibration DA for Module:%d Iteration:%d",module,iterid); + SetTitle(htitle); + + fMod = module; + fIterId = iterid; + fEvent = 0; + fEventClustered = false; + fTime = 0; + fCreateTree = false; + fCreateHist = false; + + char fname[1024]; + sprintf(fname,"AliPHOSDApi0mip_mod%d.root",module); + fTFile = TFile::Open(fname,fopt); +} +//---------------------------------------------------------------- +AliPHOSDApi0mip::AliPHOSDApi0mip(const AliPHOSDApi0mip& da): + TNamed(da), fCreateTree(false), fCreateHist(false), fMod(da.fMod), fIterId(da.fIterId), + fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0), + fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0), + fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) { + // Copy Constructor + + char fname[1024], hname[1024], htitle[1024];; + sprintf(fname,"%s.root",GetName()); + fTFile = TFile::Open(fname,"RECREATE"); + fEvent = 0; + fEventClustered = false; + fTime = 0; + + if( da.fCreateHist ){ + fCreateHist = true; + fH1Time = (TH1I*) da.fH1Time->Clone(); + fH1DigitNum = (TH1F*) da.fH1DigitNum->Clone(); + fH1ClusterNum = (TH1F*) da.fH1ClusterNum->Clone(); + fH2EneDigitId = (TH2F*) da.fH2EneDigitId->Clone(); + fH2MipDigitId = (TH2F*) da.fH2MipDigitId->Clone(); + fH2Pi0DigitId = (TH2F*) da.fH2Pi0DigitId->Clone(); + fH3Pi0AsymPt = (TH3F*) da.fH3Pi0AsymPt->Clone(); + } else { + fCreateHist = false; + fH1Time = 0; + fH1DigitNum = 0; + fH1ClusterNum = 0; + fH2EneDigitId = 0; + fH2MipDigitId = 0; + fH2Pi0DigitId = 0; + fH3Pi0AsymPt = 0; + } + if( da.fCreateTree ){ + // Create new ttree instead of copy. + fCreateTree = true; + sprintf(hname,"tevt_mod%d_iter%d",fMod,fIterId); + sprintf(htitle,"Calibration for Module:%d Iteration:%d",fMod,fIterId); + fTTree = new TTree(hname,htitle); + fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent); + } else { + fCreateTree = false; + fTTree = 0; + } +} +//------------------------------------------------------------------- +AliPHOSDApi0mip::~AliPHOSDApi0mip(){ + // Destructor + + //std::cout<<" DEBUGDEBUG :: AliPHOSDApi0mip:: deconstructor "<Write(0,TObject::kOverwrite); + fH1DigitNum->Write(0,TObject::kOverwrite); + fH1ClusterNum->Write(0,TObject::kOverwrite); + fH2EneDigitId->Write(0,TObject::kOverwrite); + fH2MipDigitId->Write(0,TObject::kOverwrite); + fH2Pi0DigitId->Write(0,TObject::kOverwrite); + fH3Pi0AsymPt->Write(0,TObject::kOverwrite); + delete fH1Time; + delete fH1DigitNum; + delete fH1ClusterNum; + delete fH2EneDigitId; + delete fH2MipDigitId; + delete fH2Pi0DigitId; + delete fH3Pi0AsymPt; + } + if( fCreateTree ){ + fTTree->Write(); + } + if( fEvent ) delete fEvent; + fTFile->Close(); + delete fTFile; +} +//------------------------------------------------------------------- +void AliPHOSDApi0mip::FillDigit(float adc,int row,int col){ + // Put new digit information in event + // + + if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent; + float energy = adc * 0.005; // 5MeV/ch is fixed. + if( energy > 0.020 ){ // Threshold at 20MeV. + //std::cout<<" AliPHOSDApi0mip::AppendDigit() DEBUG:: (row,col)=("<Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent); + fCreateTree = true; + return true; +} +//---------------------------------------------------------------- +bool AliPHOSDApi0mip::CreateHist(){ + // Create histogram routine called by the constructor only. + + char hname[1024], htitle[1024]; + + sprintf(hname,"h1_time_mod%d_iter%d",fMod,fIterId); + sprintf(htitle,"Time : Mod:%d Iter:%d",fMod,fIterId); + fH1Time = (TH1I*) gDirectory->Get(hname); + if( fH1Time>0 ){ + std::cout<<" AliPHOSDApi0mip:Warning!! Output object already exist : "<GetName()<ExecuteClustering(); + //fEventClustered = true; + //} + fEvent->SetTime(fTime); + fTTree->Fill(); + } else { + if( fEvent ) delete fEvent; + fEvent = event; + fTTree->Fill(); + fEvent = 0; + } +} +//------------------------------------------------------------------- +void AliPHOSDApi0mip::FillHist(AliPHOSDATreeEvent* event){ + // Fill all information (AliPHOSDATreeEvent) into histograms + // if event==NULL, the obtained information (=fEvent) is used. + + if(! fCreateHist ) CreateHist(); + if( event == 0 ){ + if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent; + if(! fEventClustered ){ + fEvent->ExecuteClustering(); + fEventClustered = true; + } + fEvent->SetTime(fTime); + event = fEvent; + } + + // -- Constant + int nclt1, nclt2; + int ndigits; + int ndigitsall = 0; + AliPHOSDATreeDigit digit1, digit2; + const float kDistanceFromVertex = 460.; + float dist1, dist2, px, py; + float mass, cosine, pt, asym; + + // -- Filling into Time Histogram + double start = fH1Time->GetBinContent(1); + double stop = fH1Time->GetBinContent(2); + double currenttime = (double) (event->GetTime()); + if( start == 0 && stop == 0 ) start = stop = currenttime; + if( currenttime < start ) start = currenttime; + if( currenttime > stop ) stop = currenttime; + fH1Time->SetBinContent(1,start); + fH1Time->SetBinContent(2,stop); + + // -- Start looping clusters + nclt1 = event->GetNClusters(); + fH1ClusterNum->Fill(nclt1); + //std::cout<<" AliPHOSDApi0mip::FillHisto() DEBUG:: Number of clusters: "<GetCluster(nclt1); + ndigits = clt1.GetNDigits(); + while( ndigits-- ){ + AliPHOSDATreeDigit& digit = clt1.GetDigit(ndigits); + fH2EneDigitId->Fill(digit.GetDigitId(),digit.GetEnergy()); + ndigitsall++; + } + digit1 = clt1.GetMaxDigit(); + if( clt1.GetEnergy() > 0 && clt1.GetNDigits() >= 2 ){ + fH2MipDigitId->Fill(digit1.GetDigitId(),clt1.GetEnergy()); + } + if( clt1.GetEnergy() > 0.2 ){ + nclt2 = nclt1; + while( nclt2-- ){ + AliPHOSDATreeCluster& clt2 = event->GetCluster(nclt2); + digit2 = clt2.GetMaxDigit(); + if( clt2.GetEnergy() > 0.2 ){ + // + // Following approximation is applied in order to avoid many calling "sqrt" function. + // sqrt( 1 + x^2 ) =~ 1 + 0.5x^2 (for x<<1) + // + //cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex) + // / ( kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*2.42/kDistanceFromVertex ) + // / ( kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*2.42/kDistanceFromVertex ); + // + + + dist1 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*4.84 ); + dist2 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 ); + cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex) + / dist1 / dist2; + mass = sqrt( 2.0 * clt1.GetEnergy() * clt2.GetEnergy() * (1 - cosine) ); + //std::cout<<" DEBUG: dist1:"<Fill(asym,pt,mass); + if( digit1.IsValid() && digit1.IsValid() && + digit1.GetDigitId() != digit2.GetDigitId() && + //pt > 1.5 && clt1.GetEnergy() > 0.3 && clt2.GetEnergy() > 0.3 ){ + //pt > 1.3 && clt1.GetEnergy() > 0.25 && clt2.GetEnergy() > 0.25 ){ + pt > 1.0 && clt1.GetEnergy() > 0.2 && clt2.GetEnergy() > 0.2 ){ + fH2Pi0DigitId->Fill(digit1.GetDigitId(),mass); + fH2Pi0DigitId->Fill(digit2.GetDigitId(),mass); + } + } // End of energy cut + } // End of second cluster loop (clt2) + } // End of energy cut + } // End of first cluster loop (clt1) + fH1DigitNum->Fill(ndigitsall); + // +} +//------------------------------------------------------------------- +void AliPHOSDApi0mip::Print(char* opt){ + // Print Out + + //fTFile->ls(); + //fTTree->Print(); + if( fEvent ){ + //fEvent->ExecuteClustering(); + fEvent->Print(opt); + } +} +//------------------------------------------------------------------- diff --git a/PHOS/AliPHOSDApi0mip.h b/PHOS/AliPHOSDApi0mip.h new file mode 100644 index 00000000000..d67a82cb082 --- /dev/null +++ b/PHOS/AliPHOSDApi0mip.h @@ -0,0 +1,59 @@ +// -- +// -- +// Implementation for TTree output in PHOS DA +// for calibrating energy by pi0 and MIP. +// -- +// -- Author: Hisayuki Torii (Hiroshima Univ.) +// -- +#ifndef AliPHOSDApi0mip_H +#define AliPHOSDApi0mip_H + +#include + +#include "TNamed.h" +#include "TH1I.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TH3F.h" +#include "TFile.h" +#include "TTree.h" +#include "AliPHOSDATreeEvent.h" + +class AliPHOSDApi0mip : public TNamed { + public: + AliPHOSDApi0mip(int module,int iterid=0,char* fopt="RECREATE"); + AliPHOSDApi0mip(const AliPHOSDApi0mip& da); + ~AliPHOSDApi0mip(); + + void NewEvent(); + void FillDigit(float adc,int row,int col); + void SetTime(time_t& intime){fTime=intime;}; + time_t GetTime(){return fTime;}; + void FillTree(AliPHOSDATreeEvent* event=0); + void FillHist(AliPHOSDATreeEvent* event=0); + void Print(char* opt=""); + + private: + bool CreateTree(); + bool CreateHist(); + bool fCreateTree; //! Flag of tree initialization + bool fCreateHist; //! Flag of hist initialization + int fMod; // Module ID [0-4] ([2-4] for 2009) + int fIterId; // Iteration step [0-*] + // + TFile* fTFile; //! output file + TTree* fTTree; //! output TTree + AliPHOSDATreeEvent* fEvent; //! Contents of TTree + bool fEventClustered; //! Flag for + time_t fTime; // time + TH1I* fH1Time; // x:bin1=StartTime bin2=EndTime + TH1F* fH1DigitNum; // x:Number of digits + TH1F* fH1ClusterNum; // x:Number of clusters + TH2F* fH2EneDigitId; // x:DigitId[0-3583] y:Digit Energy + TH2F* fH2MipDigitId; // x:DigitId[0-3583] y:Cluster Energy + TH2F* fH2Pi0DigitId; // x:DigitId[0-3583] y:Cluster Pair Mass + TH3F* fH3Pi0AsymPt; // x:asym y:pT(GeV/c) z:Cluster Pair Mass + + ClassDef(AliPHOSDApi0mip,1) +}; +#endif diff --git a/PHOS/AliPHOSPreprocessorPHYS.cxx b/PHOS/AliPHOSPreprocessorPHYS.cxx new file mode 100644 index 00000000000..e480bea50b3 --- /dev/null +++ b/PHOS/AliPHOSPreprocessorPHYS.cxx @@ -0,0 +1,110 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +/////////////////////////////////////////////////////////////////////////////// +// PHOS Preprocessor class. It runs by Shuttle at the end of the run, +// calculates calibration coefficients and dead/bad channels +// to be posted in OCDB +// +// Author: Hisayuki Torii, 11 August 2009 +/////////////////////////////////////////////////////////////////////////////// +#include "AliLog.h" +#include "TFile.h" +#include "TKey.h" +#include "TList.h" +#include "TString.h" +#include "TObjString.h" + +#include "AliPHOSPreprocessorPHYS.h" + +ClassImp(AliPHOSPreprocessorPHYS) + +//_______________________________________________________________________________________ +AliPHOSPreprocessorPHYS::AliPHOSPreprocessorPHYS() : AliPreprocessor("PHS",0) { + //default constructor +} +//_______________________________________________________________________________________ +AliPHOSPreprocessorPHYS::AliPHOSPreprocessorPHYS(AliShuttleInterface* shuttle): AliPreprocessor("PHS",shuttle) { + // Constructor + + AddRunType("PHYSICS"); + AddRunType("COSMICS"); // Does this exist?? +} +//_______________________________________________________________________________________ +UInt_t AliPHOSPreprocessorPHYS::Process(TMap* /*valueSet*/){ + // process data retrieved by the Shuttle + // + + TString runType = GetRunType(); + Log(Form("Run type: %s",runType.Data())); + + if( runType=="PHYSICS" || runType=="COSMICS") { + + Bool_t calibEmc_OK = CalibratePhys(); + if(calibEmc_OK) return 0; + else return 1; + } + Log(Form("Unknown run type %s. Do nothing and return OK.",runType.Data())); + return 0; +} +//_______________________________________________________________________________________ +Bool_t AliPHOSPreprocessorPHYS::CalibratePhys(){ + //process PHYSICS event retrieved by the Shuttle + // + + TList* list = 0; + list = GetFileSources(kDAQ,"PHOSDApi0mip"); + if(!list) { + Log(Form("DAQ sources list not found!")); + return kFALSE; + } + if(!list->GetEntries()) { + Log(Form("Got empty sources list. It seems PHYS DA did not produce any files!")); + return kFALSE; + } + + TIter iter(list); + TObjString *source; + + while ((source = dynamic_cast (iter.Next()))) { + AliInfo(Form("found source %s", source->String().Data())); + + TString fileName = GetFile(kDAQ, "PHOSDApi0mip", source->GetName()); + AliInfo(Form("Got filename: %s",fileName.Data())); + + TFile* file = TFile::Open(fileName); + if(!file) { + Log(Form("File %s is not opened, something goes wrong!",fileName.Data())); + return kFALSE; + } + + file->ls(); + //TList * keylist = file->GetListOfKeys(); + Int_t nkeys = file->GetNkeys(); + if(nkeys==0){ + Log(Form("Not enough (%d) for calibration.",nkeys)); + return 1; // it's not fatal! May be short run.. + } + + file->Close(); + delete file; + + + } + return kTRUE; +} +//_______________________________________________________________________________________ + + diff --git a/PHOS/AliPHOSPreprocessorPHYS.h b/PHOS/AliPHOSPreprocessorPHYS.h new file mode 100644 index 00000000000..35cd7796582 --- /dev/null +++ b/PHOS/AliPHOSPreprocessorPHYS.h @@ -0,0 +1,28 @@ +#ifndef ALIPHOSPREPROCESSORPHYS_H +#define ALIPHOSPREPROCESSORPHYS_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/////////////////////////////////////////////////////////////////////////////// +// Class AliPHOSPreprocessorPHYS +/////////////////////////////////////////////////////////////////////////////// + + +#include "AliPreprocessor.h" + +class AliPHOSPreprocessorPHYS : public AliPreprocessor { + public: + + AliPHOSPreprocessorPHYS(); + AliPHOSPreprocessorPHYS(AliShuttleInterface* shuttle); + + protected: + + virtual UInt_t Process(TMap* valueSet); + Bool_t CalibratePhys(); + + private: + + ClassDef(AliPHOSPreprocessorPHYS,1) // Preprocessor for PHYS event +}; +#endif diff --git a/PHOS/PHOSPHYSda.cxx b/PHOS/PHOSPHYSda.cxx new file mode 100644 index 00000000000..abe7e2055bd --- /dev/null +++ b/PHOS/PHOSPHYSda.cxx @@ -0,0 +1,193 @@ +/* + +PHOSPHYSda.cxx + +This program reads the DAQ data files passed as argument using the monitoring library. +It stores the required event information into an output file. +Messages on stdout are exported to DAQ log system. + +contact: Hisayuki.Torii@cern.ch + +*/ + +#include "event.h" +#include "monitor.h" +#include "daqDA.h" + +#include +#include + +#include "AliRawReader.h" +#include "AliRawReaderDate.h" +#include "AliPHOSRawDecoder.h" +#include "AliCaloAltroMapping.h" +#include "AliPHOSDApi0mip.h" + +/* Main routine + Arguments: list of DATE raw data files +*/ +int main(int argc, char **argv) { + + int status; + + /* log start of process */ + printf("PHOSPHYSda program started\n"); + + /* Time out counters*/ + int timeout = 180; // 3mins + time_t start_time = time(0); + + /* check that we got some arguments = list of files */ + if (argc<2) { + printf("Wrong number of arguments\n"); + return -1; + } + + + /* Retrieve mapping files from DAQ DB */ + const char* mapFiles[4] = {"RCU0.data","RCU1.data","RCU2.data","RCU3.data"}; + for(Int_t iFile=0; iFile<4; iFile++) { + int failed = daqDA_DB_getFile(mapFiles[iFile], mapFiles[iFile]); + if(failed) { + printf("Cannot retrieve file %s from DAQ DB. Exit.\n",mapFiles[iFile]); + return -1; + } + } + + /* Open mapping files */ + AliAltroMapping *mapping[4]; + TString path = "./"; + path += "RCU"; + TString path2; + for(Int_t i = 0; i < 4; i++) { + path2 = path; + path2 += i; + path2 += ".data"; + mapping[i] = new AliCaloAltroMapping(path2.Data()); + } + + /* report progress */ + daqDA_progressReport(10); + + /* init some counters */ + int nevents_physics=0; + int nevents_total=0; + + /* init PHOS reader*/ + AliRawReader *rawReader = NULL; + + /* Prepare DA algorithm */ + int iMod = 0; + AliPHOSDApi0mip* dapi0mip = new AliPHOSDApi0mip(iMod); + + /* read the data files */ + for(int n=1;neventTimestamp; + dapi0mip->SetTime(t); + + /* Applying into DA algorithm */ + int ix, iz, igain; + eventT=event->eventType; + if (eventT==PHYSICS_EVENT || eventT==CALIBRATION_EVENT) { + // --------------------------------------------------------------- + // User Defined Function + // --------------------------------------------------------------- + dapi0mip->NewEvent(); + rawReader = new AliRawReaderDate((void*)event); + //AliPHOSRawDecoderv1 dc(rawReader,mapping); + AliPHOSRawDecoder dc(rawReader,mapping); + dc.SubtractPedestals(false); + while(dc.NextDigit()) { + //printf("."); + ix = dc.GetRow() - 1; + iz = dc.GetColumn() - 1; + if(dc.IsLowGain()) igain = 0; else igain = 1; + if( igain == 1 ){ + dapi0mip->FillDigit(dc.GetEnergy(),ix,iz); + } + //dc.GetTime(); + } + delete rawReader; + //dapi0mip->FillHist(); // This is not necesarry in this DA algorithm + dapi0mip->FillTree(); + //if( nevents_total%1000 == 0 ){ // Dump for debugging + //dapi0mip->Print(); + //} + // --------------------------------------------------------------- + nevents_physics++; + } + nevents_total++; + + /* Check the time out */ + if( nevents_total%1000 == 0 ){ + time_t current_time = time(0); + if( current_time - start_time > timeout ){ + free(event); + printf(" Warning: Exit due to the processing time exceed the limitation of %d sec\n",timeout); + n = argc; + break; + } + } + + /* free resources */ + free(event); + } + + } + + /* report progress */ + daqDA_progressReport(90); + + /* Delete DA algorithm for saving output into a file */ + delete dapi0mip; + + /* Store output files to the File Exchange Server */ + char localfile[1024]; + char remotefile[1024]; + sprintf(localfile,"AliPHOSDApi0mip_mod%d.root",iMod); + sprintf(remotefile,"PHOSDApi0mip",iMod); + daqDA_FES_storeFile(localfile,remotefile); + + /* report progress */ + daqDA_progressReport(100); + + return status; +} + -- 2.39.3