From e2a4d72cfff407434bf9f6fdd9a29477c143e79e Mon Sep 17 00:00:00 2001 From: marian Date: Fri, 24 Oct 2008 16:52:23 +0000 Subject: [PATCH] Adding cut on isolations Adding the cluster shape (Adam Matyja) --- TPC/AliTPCclusterKr.cxx | 418 +++++--- TPC/AliTPCclusterKr.h | 177 ++-- TPC/AliTPCclustererKr.cxx | 1904 +++++++++++++++++++------------------ TPC/AliTPCclustererKr.h | 325 +++---- TPC/FindKrClusters.C | 273 +++--- TPC/FindKrClustersRaw.C | 330 +++---- 6 files changed, 1803 insertions(+), 1624 deletions(-) diff --git a/TPC/AliTPCclusterKr.cxx b/TPC/AliTPCclusterKr.cxx index 921c538f78d..33e0e292b92 100644 --- a/TPC/AliTPCclusterKr.cxx +++ b/TPC/AliTPCclusterKr.cxx @@ -1,141 +1,277 @@ -/************************************************************************** - * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * * - * Author: The ALICE Off-line Project. * - * Contributors are mentioned in the code where appropriate. * - * * - * Permission to use, copy, modify and distribute this software and its * - * documentation strictly for non-commercial purposes is hereby granted * - * without fee, provided that the above copyright notice appears in all * - * copies and that both the copyright notice and this permission notice * - * appear in the supporting documentation. The authors make no claims * - * about the suitability of this software for any purpose. It is * - * provided "as is" without express or implied warranty. * - **************************************************************************/ - -/* $Id: AliTPCclusterKr.cxx,v 1.7 2008/01/22 17:24:53 matyja Exp $ */ - -//----------------------------------------------------------------- -// Implementation of the TPC Kr cluster class -// -// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl -//----------------------------------------------------------------- - -#include "AliTPCclusterKr.h" -#include "AliCluster.h" -#include "AliTPCvtpr.h" -#include "TObjArray.h" - -ClassImp(AliTPCclusterKr) - - -AliTPCclusterKr::AliTPCclusterKr() -:AliCluster(), - fMax(), - fADCcluster(0), - fSec(0), - fNPads(0), - fNRows(0), - fSize(0), - fCenterX(0), - fCenterY(0), - fCenterT(0), - fCluster(0) -{ -// -// default constructor -// - fCluster=new TObjArray(); -} - -AliTPCclusterKr::AliTPCclusterKr(const AliTPCclusterKr ¶m) -:AliCluster(param), - fMax(), - fADCcluster(0), - fSec(0), - fNPads(0), - fNRows(0), - fSize(0), - fCenterX(0), - fCenterY(0), - fCenterT(0), - fCluster(0) -{ -// -// copy constructor -// - fADCcluster = param.fADCcluster; - fSec = param.fSec ; - fNPads = param.fNPads; - fNRows = param.fNRows; - fMax = param.fMax; - // fCluster = param.fCluster; - fCenterX = param.fCenterX; - fCenterY = param.fCenterY; - fCenterT = param.fCenterT; - fCluster=new TObjArray(*(param.fCluster)); - fSize = param.fSize; -} - -AliTPCclusterKr &AliTPCclusterKr::operator = (const AliTPCclusterKr & param) -{ - // - // assignment operator - // - (AliCluster&)(*this) = (AliCluster&)param; - fADCcluster = param.fADCcluster; - fSec = param.fSec ; - fNPads = param.fNPads; - fNRows = param.fNRows; - fMax = param.fMax; - // fCluster=param.fCluster; - fCenterX = param.fCenterX; - fCenterY = param.fCenterY; - fCenterT = param.fCenterT; - delete fCluster; - fCluster=new TObjArray(*(param.fCluster)); - fSize=param.fSize; - return (*this); -} - -AliTPCclusterKr::~AliTPCclusterKr() -{ - // - // destructor - // - if(fCluster) { - fCluster->SetOwner(kTRUE); - fCluster->Delete(); - delete fCluster; - } - fCluster=0; -} - -////____________________________________________________________________________ -void AliTPCclusterKr::SetCenter(){ - // - // calculate geometrical center of the cluster - // - Double_t rX=0; - Double_t rY=0; - Double_t rT=0; - - Short_t adc; - fADCcluster=0; - for(Int_t iter = 0; iter < fCluster->GetEntriesFast(); ++iter) { - AliTPCvtpr *iclus=(AliTPCvtpr *)fCluster->At(iter); - - //for( std::vector::iterator iclus = fCluster.begin(); - //iclus != fCluster.end(); ++iclus ) { - adc = (iclus)->GetAdc(); - fADCcluster+=adc; - rX += ((iclus)->GetX() * adc); - rY += ((iclus)->GetY() * adc); - rT += ((iclus)->GetT() * adc); - } - fCenterX=rX/fADCcluster; - fCenterY=rY/fADCcluster; - fCenterT=rT/fADCcluster; - - return; -} +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id: AliTPCclusterKr.cxx,v 1.7 2008/01/22 17:24:53 matyja Exp $ */ + +//----------------------------------------------------------------- +// Implementation of the TPC Kr cluster class +// +// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl +//----------------------------------------------------------------- + +#include "AliTPCclusterKr.h" +#include "AliCluster.h" +#include "AliTPCvtpr.h" +#include "TObjArray.h" +//#include "TH1F.h" +#include "TMath.h" +#include "TArrayI.h" + +ClassImp(AliTPCclusterKr) + + +AliTPCclusterKr::AliTPCclusterKr() +:AliCluster(), + fMax(), + fADCcluster(0), + fSec(0), + fNPads(0), + fNRows(0), + fTimebins1D(0), + fPads1D(0), + fPadRMS(0), + fRowRMS(0), + fTimebinRMS(0), + fSize(0), + fCenterX(0), + fCenterY(0), + fCenterT(0), + fCluster(0) +{ +// +// default constructor +// + fCluster=new TObjArray(); +} + +AliTPCclusterKr::AliTPCclusterKr(const AliTPCclusterKr ¶m) +:AliCluster(param), + fMax(), + fADCcluster(0), + fSec(0), + fNPads(0), + fNRows(0), + fTimebins1D(0), + fPads1D(0), + fPadRMS(0), + fRowRMS(0), + fTimebinRMS(0), + fSize(0), + fCenterX(0), + fCenterY(0), + fCenterT(0), + fCluster(0) +{ +// +// copy constructor +// + fADCcluster = param.fADCcluster; + fSec = param.fSec ; + fNPads = param.fNPads; + fNRows = param.fNRows; + fMax = param.fMax; + // fCluster = param.fCluster; + fCenterX = param.fCenterX; + fCenterY = param.fCenterY; + fCenterT = param.fCenterT; + fCluster=new TObjArray(*(param.fCluster)); + fSize = param.fSize; + fTimebins1D = param.fTimebins1D; + fPads1D = param.fPads1D; + fPadRMS = param.fPadRMS; + fRowRMS = param.fRowRMS; + fTimebinRMS = param.fTimebinRMS; +} + +AliTPCclusterKr &AliTPCclusterKr::operator = (const AliTPCclusterKr & param) +{ + // + // assignment operator + // + (AliCluster&)(*this) = (AliCluster&)param; + fADCcluster = param.fADCcluster; + fSec = param.fSec ; + fNPads = param.fNPads; + fNRows = param.fNRows; + fMax = param.fMax; + // fCluster=param.fCluster; + fCenterX = param.fCenterX; + fCenterY = param.fCenterY; + fCenterT = param.fCenterT; + delete fCluster; + fCluster=new TObjArray(*(param.fCluster)); + fSize=param.fSize; + fTimebins1D = param.fTimebins1D; + fPads1D = param.fPads1D; + fPadRMS = param.fPadRMS; + fRowRMS = param.fRowRMS; + fTimebinRMS = param.fTimebinRMS; + return (*this); +} + +AliTPCclusterKr::~AliTPCclusterKr() +{ + // + // destructor + // + if(fCluster) { + fCluster->SetOwner(kTRUE); + fCluster->Delete(); + delete fCluster; + } + fCluster=0; +} + +////____________________________________________________________________________ +void AliTPCclusterKr::SetCenter(){ + // + // calculate geometrical center of the cluster + // + Double_t rX=0; + Double_t rY=0; + Double_t rT=0; + + Short_t adc; + fADCcluster=0; + for(Int_t iter = 0; iter < fCluster->GetEntriesFast(); ++iter) { + AliTPCvtpr *iclus=(AliTPCvtpr *)fCluster->At(iter); + + //for( std::vector::iterator iclus = fCluster.begin(); + //iclus != fCluster.end(); ++iclus ) { + adc = (iclus)->GetAdc(); + fADCcluster+=adc; + rX += ((iclus)->GetX() * adc); + rY += ((iclus)->GetY() * adc); + rT += ((iclus)->GetT() * adc); + } + fCenterX=rX/fADCcluster; + fCenterY=rY/fADCcluster; + fCenterT=rT/fADCcluster; + + return; +} + +void AliTPCclusterKr::SetPadRMS(){ + // + // calculate RMS in pad direction + // + // TH1F *histo= new TH1F("","",200,0,200); + TArrayI *array= new TArrayI(fCluster->GetEntriesFast()); + for(Int_t i=0;iGetEntriesFast();i++) + { + array->SetAt(((AliTPCvtpr *)(fCluster->At(i)))->GetPad(),i); + //histo->Fill( ((AliTPCvtpr *)(fCluster->At(i)))->GetPad() ); + } + // fPadRMS=histo->GetRMS(); + fPadRMS=TMath::RMS(array->GetSize(),array->GetArray()); + // delete histo; + delete array; + return; +} + +void AliTPCclusterKr::SetRowRMS(){ + // + // calculate RMS in row direction + // + TArrayI *array= new TArrayI(fCluster->GetEntriesFast()); + // TH1F *histo= new TH1F("","",120,0,120); + for(Int_t i=0;iGetEntriesFast();i++) + { + array->SetAt(((AliTPCvtpr *)(fCluster->At(i)))->GetRow(),i); + // histo->Fill( ((AliTPCvtpr *)(fCluster->At(i)))->GetRow() ); + } + // fRowRMS=histo->GetRMS(); + fRowRMS=TMath::RMS(array->GetSize(),array->GetArray()); + // delete histo; + delete array; + return; +} + +void AliTPCclusterKr::SetTimebinRMS(){ + // + // calculate RMS in timebin direction + // + TArrayI *array= new TArrayI(fCluster->GetEntriesFast()); + // TH1F *histo= new TH1F("","",1000,0,1000); + for(Int_t i=0;iGetEntriesFast();i++) + { + array->SetAt(((AliTPCvtpr *)(fCluster->At(i)))->GetTime(),i); + // histo->Fill( ((AliTPCvtpr *)(fCluster->At(i)))->GetTime() ); + } + fTimebinRMS=TMath::RMS(array->GetSize(),array->GetArray()); + //histo->GetRMS(); + // delete histo; + delete array; + return; +} + +void AliTPCclusterKr::SetRMS(){ + // + // calculate RMS in pad,row,timebin direction + // + TArrayI *arrayPad = new TArrayI(fCluster->GetEntriesFast()); + TArrayI *arrayRow = new TArrayI(fCluster->GetEntriesFast()); + TArrayI *arrayTime= new TArrayI(fCluster->GetEntriesFast()); + // TH1F *histoPad= new TH1F("p","p",200,0,200); + // TH1F *histoRow= new TH1F("r","r",120,0,120); + // TH1F *histoTime= new TH1F("t","t",1000,0,1000); + for(Int_t i=0;iGetEntriesFast();i++) + { + arrayPad->SetAt(((AliTPCvtpr *)(fCluster->At(i)))->GetPad(),i); + arrayRow->SetAt(((AliTPCvtpr *)(fCluster->At(i)))->GetRow(),i); + arrayTime->SetAt(((AliTPCvtpr *)(fCluster->At(i)))->GetTime(),i); + + //histoPad->Fill( ((AliTPCvtpr *)(fCluster->At(i)))->GetPad() ); + //histoRow->Fill( ((AliTPCvtpr *)(fCluster->At(i)))->GetRow() ); + //histoTime->Fill( ((AliTPCvtpr *)(fCluster->At(i)))->GetTime() ); + } + // fPadRMS=histoPad->GetRMS(); + fPadRMS=TMath::RMS(arrayPad->GetSize(),arrayPad->GetArray()); + fRowRMS=TMath::RMS(arrayRow->GetSize(),arrayRow->GetArray()); + //histoRow->GetRMS(); + fTimebinRMS=TMath::RMS(arrayTime->GetSize(),arrayTime->GetArray()); + //histoTime->GetRMS(); + + delete arrayPad; + delete arrayRow; + delete arrayTime; + // delete histoPad; + // delete histoRow; + // delete histoTime; + + return; +} + + +void AliTPCclusterKr::Set1D(){ + // + // + // + Short_t maxTime=0; + Short_t minTime=1000; + Short_t maxPad=0; + Short_t minPad=1000; + + for(Int_t i=0;iGetEntriesFast();i++) + { + if(((AliTPCvtpr *)(fCluster->At(i)))->GetPad()>maxPad)maxPad =((AliTPCvtpr *)(fCluster->At(i)))->GetPad(); + if(((AliTPCvtpr *)(fCluster->At(i)))->GetPad()At(i)))->GetPad(); + if(((AliTPCvtpr *)(fCluster->At(i)))->GetTime()>maxTime)maxTime=((AliTPCvtpr *)(fCluster->At(i)))->GetTime(); + if(((AliTPCvtpr *)(fCluster->At(i)))->GetTime()At(i)))->GetTime(); + } + fPads1D=maxPad-minPad+1; + fTimebins1D=maxTime-minTime+1; + return; +} diff --git a/TPC/AliTPCclusterKr.h b/TPC/AliTPCclusterKr.h index c5945f80493..281ba71f031 100644 --- a/TPC/AliTPCclusterKr.h +++ b/TPC/AliTPCclusterKr.h @@ -1,75 +1,102 @@ -#ifndef ALITPCCLUSTERKR_H -#define ALITPCCLUSTERKR_H -/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -/* $Id: AliTPCclusterKr.h,v 1.8 2008/01/22 16:07:15 matyja Exp $ */ - -//------------------------------------------------------- -// TPC Kr Cluster Class -// -// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl -//------------------------------------------------------- - -#include "AliCluster.h" -#include "TObjArray.h" -#include "AliTPCvtpr.h" - - -//_____________________________________________________________________________ -class AliTPCclusterKr: public AliCluster{ -public: - AliTPCclusterKr(); - AliTPCclusterKr(const AliTPCclusterKr & param);//copy constructor - AliTPCclusterKr &operator = (const AliTPCclusterKr & param); - virtual ~AliTPCclusterKr(); - - virtual void SetCenter();//set center of the cluster weighted by charge - - virtual void SetMax(AliTPCvtpr q){fMax=q;}//set values of max. in cluster - virtual void SetADCcluster(Int_t q){fADCcluster=q;} - virtual void SetSec(Short_t q){fSec=q;} - virtual void SetNPads(Short_t q){fNPads=q;} - virtual void SetNRows(Short_t q){fNRows=q;} - virtual void SetSize(){fSize=fCluster->GetEntriesFast();} - virtual void SetCenterX(Double_t q){fCenterX=q;} - virtual void SetCenterY(Double_t q){fCenterY=q;} - virtual void SetCenterT(Double_t q){fCenterT=q;} - //void AddDigitToCluster(AliTPCvtpr *q){fCluster.push_back(q);} - virtual void AddDigitToCluster(AliTPCvtpr *q){ - fCluster->AddLast(q); - //fCluster->Compress(); - } - - AliTPCvtpr GetMax() const {return fMax;} - Int_t GetADCcluster() const {return fADCcluster;} - Short_t GetSec() const {return fSec;} - Short_t GetNPads() const {return fNPads;} - Short_t GetNRows() const {return fNRows;} - Short_t GetSize() const {return fSize;} - Double_t GetCenterX() const {return fCenterX;} - Double_t GetCenterY() const {return fCenterY;} - Double_t GetCenterT() const {return fCenterT;} - AliTPCvtpr *GetDigitFromCluster(Int_t i) const {return (AliTPCvtpr*)fCluster->At(i);} - -private: - AliTPCvtpr fMax;//max (ADC,timebin,pad,row) in cluster - Int_t fADCcluster; //ADC of cluster - Short_t fSec; //sector of the cluster - Short_t fNPads; //number of pads in cluster - Short_t fNRows; //number of rows in cluster - Short_t fSize; //size of vector - Double_t fCenterX;// X coordinate of the cluster center in cm - Double_t fCenterY;// Y coordinate of the cluster center in cm - Double_t fCenterT;// time coordinate of the cluster center in timebins - //std::vector< AliTPCvtpr*> fCluster;//cluster contents(adc,nt,np,nr) - TObjArray *fCluster;//cluster contents(adc,nt,np,nr) - - - ClassDef(AliTPCclusterKr,4) // Time Projection Chamber Kr clusters -}; - - -#endif - - +#ifndef ALITPCCLUSTERKR_H +#define ALITPCCLUSTERKR_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id: AliTPCclusterKr.h,v 1.8 2008/01/22 16:07:15 matyja Exp $ */ + +//------------------------------------------------------- +// TPC Kr Cluster Class +// +// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl +//------------------------------------------------------- + +#include "AliCluster.h" +#include "TObjArray.h" +#include "AliTPCvtpr.h" +//#include "TH1F.h" +#include "TMath.h" +#include "TArrayI.h" + +//_____________________________________________________________________________ +class AliTPCclusterKr: public AliCluster{ +public: + AliTPCclusterKr(); + AliTPCclusterKr(const AliTPCclusterKr & param);//copy constructor + AliTPCclusterKr &operator = (const AliTPCclusterKr & param); + virtual ~AliTPCclusterKr(); + + virtual void SetCenter();//set center of the cluster weighted by charge + + virtual void SetMax(AliTPCvtpr q){fMax=q;}//set values of max. in cluster + virtual void SetADCcluster(Int_t q){fADCcluster=q;} + virtual void SetSec(Short_t q){fSec=q;} + virtual void SetNPads(Short_t q){fNPads=q;} + virtual void SetNRows(Short_t q){fNRows=q;} + virtual void SetSize(){fSize=fCluster->GetEntriesFast();} + virtual void SetCenterX(Double_t q){fCenterX=q;} + virtual void SetCenterY(Double_t q){fCenterY=q;} + virtual void SetCenterT(Double_t q){fCenterT=q;} + + virtual void SetTimebins1D(Short_t q){fTimebins1D=q;} + virtual void SetPads1D(Short_t q){fPads1D=q;} + virtual void Set1D(); + virtual void SetPadRMS(Double_t q){fPadRMS=q;} + virtual void SetRowRMS(Double_t q){fRowRMS=q;} + virtual void SetTimebinRMS(Double_t q){fTimebinRMS=q;} + virtual void SetPadRMS(); + virtual void SetRowRMS(); + virtual void SetTimebinRMS(); + virtual void SetRMS(); + //void AddDigitToCluster(AliTPCvtpr *q){fCluster.push_back(q);} + virtual void AddDigitToCluster(AliTPCvtpr *q){ + fCluster->AddLast(q); + //fCluster->Compress(); + } + + AliTPCvtpr GetMax() const {return fMax;} + Int_t GetADCcluster() const {return fADCcluster;} + Short_t GetSec() const {return fSec;} + Short_t GetNPads() const {return fNPads;} + Short_t GetNRows() const {return fNRows;} + Short_t GetSize() const {return fSize;} + + Short_t GetTimebins1D(){return fTimebins1D;} + Short_t GetPads1D(){return fPads1D;} + Double_t GetPadRMS(){return fPadRMS;} + Double_t GetRowRMS(){return fRowRMS;} + Double_t GetTimebinRMS(){return fTimebinRMS;} + + Double_t GetCenterX() const {return fCenterX;} + Double_t GetCenterY() const {return fCenterY;} + Double_t GetCenterT() const {return fCenterT;} + AliTPCvtpr *GetDigitFromCluster(Int_t i) const {return (AliTPCvtpr*)fCluster->At(i);} + +private: + AliTPCvtpr fMax;//max (ADC,timebin,pad,row) in cluster + Int_t fADCcluster; //ADC of cluster + Short_t fSec; //sector of the cluster + Short_t fNPads; //number of pads in cluster + Short_t fNRows; //number of rows in cluster or row max - min + + Short_t fTimebins1D; //Timebin max - min + Short_t fPads1D; //Pad max - min + Double_t fPadRMS; //Pad RMS + Double_t fRowRMS; //Row RMS + Double_t fTimebinRMS; //Timebin RMS + + Short_t fSize; //size of vector + Double_t fCenterX;// X coordinate of the cluster center in cm + Double_t fCenterY;// Y coordinate of the cluster center in cm + Double_t fCenterT;// time coordinate of the cluster center in timebins + //std::vector< AliTPCvtpr*> fCluster;//cluster contents(adc,nt,np,nr) + TObjArray *fCluster;//cluster contents(adc,nt,np,nr) + + + ClassDef(AliTPCclusterKr,5) // Time Projection Chamber Kr clusters +}; + + +#endif + + diff --git a/TPC/AliTPCclustererKr.cxx b/TPC/AliTPCclustererKr.cxx index bac6dd9a0ea..b093b6a79a0 100644 --- a/TPC/AliTPCclustererKr.cxx +++ b/TPC/AliTPCclustererKr.cxx @@ -1,945 +1,959 @@ -/************************************************************************** - * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * * - * Author: The ALICE Off-line Project. * - * Contributors are mentioned in the code where appropriate. * - * * - * Permission to use, copy, modify and distribute this software and its * - * documentation strictly for non-commercial purposes is hereby granted * - * without fee, provided that the above copyright notice appears in all * - * copies and that both the copyright notice and this permission notice * - * appear in the supporting documentation. The authors make no claims * - * about the suitability of this software for any purpose. It is * - * provided "as is" without express or implied warranty. * - **************************************************************************/ - -/* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */ - -//----------------------------------------------------------------- -// Implementation of the TPC Kr cluster class -// -// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl -//----------------------------------------------------------------- - -/* -Instruction - how to use that: -There are two macros prepared. One is for preparing clusters from MC -samples: FindKrClusters.C. The output is kept in TPC.RecPoints.root. -The other macro is prepared for data analysis: FindKrClustersRaw.C. -The output is created for each processed file in root file named adc.root. -For each data subsample the same named file is created. So be careful -do not overwrite them. - -Additional selection criteria to select the GOLD cluster -Example: -// open file with clusters -TFile f("Krypton.root"); -TTree * tree = (TTree*)f.Get("Kr") -TCut cutR0("cutR0","fADCcluster/fSize<100"); // adjust it according v seetings - -TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal -TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.2"); // digital noise removal -TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal -TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks -TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1; -This values are typical values to be applied in selectors - - -* -**** MC **** -* - -To run clusterizaton for MC type: -.x FindKrClusters.C - -If you don't want to use the standard selection criteria then you -have to do following: - -// load the standard setup -AliRunLoader* rl = AliRunLoader::Open("galice.root"); -AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader"); -tpcl->LoadDigits(); -rl->LoadgAlice(); -gAlice=rl->GetAliRun(); -TDirectory *cwd = gDirectory; -AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC"); -Int_t ver = tpc->IsVersion(); -rl->CdGAFile(); -AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60"); -AliTPCDigitsArray *digarr=new AliTPCDigitsArray; -digarr->Setup(param); -cwd->cd(); - -//loop over events -Int_t nevmax=rl->GetNumberOfEvents(); -for(Int_t nev=0;nevGetEvent(nev); - TTree* input_tree= tpcl->TreeD();//load tree with digits - digarr->ConnectTree(input_tree); - TTree *output_tree =tpcl->TreeR();//load output tree - - AliTPCclustererKr *clusters = new AliTPCclustererKr(); - clusters->SetParam(param); - clusters->SetInput(input_tree); - clusters->SetOutput(output_tree); - clusters->SetDigArr(digarr); - -//If you want to change the cluster finder parameters for MC there are -//several of them: - -//1. signal threshold (everything below the given number is treated as 0) - clusters->SetMinAdc(3); - -//2. number of neighbouring timebins to be considered - clusters->SetMinTimeBins(2); - -//3. distance of the cluster center to the center of a pad in pad-padrow plane -//(in cm). Remenber that this is still quantified by pad size. - clusters->SetMaxPadRangeCm(2.5); - -//4. distance of the cluster center to the center of a padrow in pad-padrow -//plane (in cm). Remenber that this is still quantified by pad size. - clusters->SetMaxRowRangeCm(3.5); - -//5. distance of the cluster center to the max time bin on a pad (in tackts) -//ie. fabs(centerT - time)<7 - clusters->SetMaxTimeRange(7); - -//6. cut reduce peak at 0. There are noises which appear mostly as two -//timebins on one pad. - clusters->SetValueToSize(3.5); - - - clusters->finderIO(); - tpcl->WriteRecPoints("OVERWRITE"); -} -delete rl;//cleans everything - -* -********* DATA ********* -* - -To run clusterizaton for DATA for file named raw_data.root type: -.x FindKrClustersRaw.C("raw_data.root") - -If you want to change some criteria do the following: - -// -// remove Altro warnings -// -AliLog::SetClassDebugLevel("AliTPCRawStream",-5); -AliLog::SetClassDebugLevel("AliRawReaderDate",-5); -AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5); -AliLog::SetModuleDebugLevel("RAW",-5); - -// -// Get database with noises -// -// char *ocdbpath = gSystem->Getenv("OCDB_PATH"); -char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB"; -if (ocdbpath==0){ -ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/"; -} -AliCDBManager * man = AliCDBManager::Instance(); -man->SetDefaultStorage(ocdbpath); -man->SetRun(0); -AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); -AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); - -//define tree -TFile *hfile=new TFile("adc.root","RECREATE","ADC file"); -// Create a ROOT Tree -TTree *mytree = new TTree("Kr","Krypton cluster tree"); - -//define infput file -const char *fileName="data.root"; -AliRawReader *reader = new AliRawReaderRoot(fileName); -//AliRawReader *reader = new AliRawReaderDate(fileName); -reader->Reset(); -AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader); -stream->SelectRawData("TPC"); - -//one general output -AliTPCclustererKr *clusters = new AliTPCclustererKr(); -clusters->SetOutput(mytree); -clusters->SetRecoParam(0);//standard reco parameters -AliTPCParamSR *param=new AliTPCParamSR(); -clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.) - -//set cluster finder parameters (from data): -//1. zero suppression parameter - clusters->SetZeroSup(param->GetZeroSup()); - -//2. first bin - clusters->SetFirstBin(60); - -//3. last bin - clusters->SetLastBin(950); - -//4. maximal noise - clusters->SetMaxNoiseAbs(2); - -//5. maximal amount of sigma of noise - clusters->SetMaxNoiseSigma(3); - -//The remaining parameters are the same paramters as for MC (see MC section -//points 1-6) - clusters->SetMinAdc(3); - clusters->SetMinTimeBins(2); - clusters->SetMaxPadRangeCm(2.5); - clusters->SetMaxRowRangeCm(3.5); - clusters->SetMaxTimeRange(7); - clusters->SetValueToSize(3.5); - -while (reader->NextEvent()) { - clusters->FinderIO(reader); -} - -hfile->Write(); -hfile->Close(); -delete stream; - - -*/ - -#include "AliTPCclustererKr.h" -#include "AliTPCclusterKr.h" -//#include -#include -#include "TObject.h" -#include "AliPadMax.h" -#include "AliSimDigits.h" -#include "AliTPCv4.h" -#include "AliTPCParam.h" -#include "AliTPCDigitsArray.h" -#include "AliTPCvtpr.h" -#include "AliTPCClustersRow.h" -#include "TTree.h" -#include "TH1F.h" -#include "TH2F.h" -#include "TTreeStream.h" - -#include "AliTPCTransform.h" - -//used in raw data finder -#include "AliTPCROC.h" -#include "AliTPCCalPad.h" -#include "AliTPCAltroMapping.h" -#include "AliTPCcalibDB.h" -#include "AliTPCRawStream.h" -#include "AliTPCRecoParam.h" -#include "AliTPCReconstructor.h" -#include "AliRawReader.h" -#include "AliTPCCalROC.h" - -ClassImp(AliTPCclustererKr) - - -AliTPCclustererKr::AliTPCclustererKr() - :TObject(), - fRawData(kFALSE), - fInput(0), - fOutput(0), - fParam(0), - fDigarr(0), - fRecoParam(0), - fZeroSup(2), - fFirstBin(60), - fLastBin(950), - fMaxNoiseAbs(2), - fMaxNoiseSigma(3), - fMinAdc(3), - fMinTimeBins(2), -// fMaxPadRange(4), -// fMaxRowRange(3), - fMaxTimeRange(7), - fValueToSize(3.5), - fMaxPadRangeCm(2.5), - fMaxRowRangeCm(3.5), - fDebugLevel(-1), - fHistoRow(0), - fHistoPad(0), - fHistoTime(0), - fHistoRowPad(0) -{ -// -// default constructor -// -} - -AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m) - :TObject(), - fRawData(kFALSE), - fInput(0), - fOutput(0), - fParam(0), - fDigarr(0), - fRecoParam(0), - fZeroSup(2), - fFirstBin(60), - fLastBin(950), - fMaxNoiseAbs(2), - fMaxNoiseSigma(3), - fMinAdc(3), - fMinTimeBins(2), -// fMaxPadRange(4), -// fMaxRowRange(3), - fMaxTimeRange(7), - fValueToSize(3.5), - fMaxPadRangeCm(2.5), - fMaxRowRangeCm(3.5), - fDebugLevel(-1), - fHistoRow(0), - fHistoPad(0), - fHistoTime(0), - fHistoRowPad(0) -{ -// -// copy constructor -// - fParam = param.fParam; - fRecoParam = param.fRecoParam; - fRawData = param.fRawData; - fInput = param.fInput ; - fOutput = param.fOutput; - fDigarr = param.fDigarr; - fZeroSup = param.fZeroSup ; - fFirstBin = param.fFirstBin ; - fLastBin = param.fLastBin ; - fMaxNoiseAbs = param.fMaxNoiseAbs ; - fMaxNoiseSigma = param.fMaxNoiseSigma ; - fMinAdc = param.fMinAdc; - fMinTimeBins = param.fMinTimeBins; -// fMaxPadRange = param.fMaxPadRange ; -// fMaxRowRange = param.fMaxRowRange ; - fMaxTimeRange = param.fMaxTimeRange; - fValueToSize = param.fValueToSize; - fMaxPadRangeCm = param.fMaxPadRangeCm; - fMaxRowRangeCm = param.fMaxRowRangeCm; - fDebugLevel = param.fDebugLevel; - fHistoRow = param.fHistoRow ; - fHistoPad = param.fHistoPad ; - fHistoTime = param.fHistoTime; - fHistoRowPad = param.fHistoRowPad; - -} - -AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param) -{ - // - // assignment operator - // - fParam = param.fParam; - fRecoParam = param.fRecoParam; - fRawData = param.fRawData; - fInput = param.fInput ; - fOutput = param.fOutput; - fDigarr = param.fDigarr; - fZeroSup = param.fZeroSup ; - fFirstBin = param.fFirstBin ; - fLastBin = param.fLastBin ; - fMaxNoiseAbs = param.fMaxNoiseAbs ; - fMaxNoiseSigma = param.fMaxNoiseSigma ; - fMinAdc = param.fMinAdc; - fMinTimeBins = param.fMinTimeBins; -// fMaxPadRange = param.fMaxPadRange ; -// fMaxRowRange = param.fMaxRowRange ; - fMaxTimeRange = param.fMaxTimeRange; - fValueToSize = param.fValueToSize; - fMaxPadRangeCm = param.fMaxPadRangeCm; - fMaxRowRangeCm = param.fMaxRowRangeCm; - fDebugLevel = param.fDebugLevel; - fHistoRow = param.fHistoRow ; - fHistoPad = param.fHistoPad ; - fHistoTime = param.fHistoTime; - fHistoRowPad = param.fHistoRowPad; - return (*this); -} - -AliTPCclustererKr::~AliTPCclustererKr() -{ - // - // destructor - // - delete fOutput; -} - -void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam) -{ - // - // set reconstruction parameters - // - if (recoParam) { - fRecoParam = recoParam; - }else{ - //set default parameters if not specified - fRecoParam = AliTPCReconstructor::GetRecoParam(); - if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam(); - } - return; -} - - -////____________________________________________________________________________ -//// I/O -void AliTPCclustererKr::SetInput(TTree * tree) -{ - // - // set input tree with digits - // - fInput = tree; - if (!fInput->GetBranch("Segment")){ - cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n"; - fInput=0; - return; - } -} - -void AliTPCclustererKr::SetOutput(TTree * /*tree*/) -{ - // - //dummy - // - fOutput = new TTreeSRedirector("Krypton.root"); -} - -////____________________________________________________________________________ -//// with new I/O -Int_t AliTPCclustererKr::FinderIO() -{ - // Krypton cluster finder for simulated events from MC - - if (!fInput) { - Error("Digits2Clusters", "input tree not initialised"); - return 10; - } - - if (!fOutput) { - Error("Digits2Clusters", "output tree not initialised"); - return 11; - } - - FindClusterKrIO(); - return 0; -} - -Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader) -{ - // Krypton cluster finder for the TPC raw data - // - // fParam must be defined before - - if(rawReader)fRawData=kTRUE; //set flag to data - - if (!fOutput) { - Error("Digits2Clusters", "output tree not initialised"); - return 11; - } - - fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param - // used later for memory allocation - - Bool_t isAltro=kFALSE; - - AliTPCROC * roc = AliTPCROC::Instance(); - AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); - AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); - // - AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping); - - const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors - const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors - const Int_t kNS = kNIS + kNOS;//all sectors - - - //crate TPC view - AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim - digarr->Setup(fParam);//as usually parameters - - // - // Loop over sectors - // - for(Int_t iSec = 0; iSec < kNS; iSec++) { - AliTPCCalROC * noiseROC; - if(noiseTPC==0x0){ - noiseROC =new AliTPCCalROC(iSec);//noise=0 - } - else{ - noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector - } - Int_t nRows = 0; //number of rows in sector - Int_t nDDLs = 0; //number of DDLs - Int_t indexDDL = 0; //DDL index - if (iSec < kNIS) { - nRows = fParam->GetNRowLow(); - nDDLs = 2; - indexDDL = iSec * 2; - }else { - nRows = fParam->GetNRowUp(); - nDDLs = 4; - indexDDL = (iSec-kNIS) * 4 + kNIS * 2; - } - - // - // Load the raw data for corresponding DDLs - // - rawReader->Reset(); - rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); - - if(input.Next()) { - isAltro=kTRUE; - // Allocate memory for rows in sector (pads(depends on row) x timebins) - for(Int_t iRow = 0; iRow < nRows; iRow++) { - digarr->CreateRow(iSec,iRow); - }//end loop over rows - } - rawReader->Reset(); - rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); - - // - // Begin loop over altro data - // - while (input.Next()) { - - //check sector consistency - if (input.GetSector() != iSec) - AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector())); - - Int_t iRow = input.GetRow(); - Int_t iPad = input.GetPad(); - Int_t iTimeBin = input.GetTime(); - - // - if(fDebugLevel==72){ - fHistoRow->Fill(iRow); - fHistoPad->Fill(iPad); - fHistoTime->Fill(iTimeBin); - fHistoRowPad->Fill(iPad,iRow); - }else if(fDebugLevel>=0&&fDebugLevel<72){ - if(iSec==fDebugLevel){ - fHistoRow->Fill(iRow); - fHistoPad->Fill(iPad); - fHistoTime->Fill(iTimeBin); - fHistoRowPad->Fill(iPad,iRow); - } - }else if(fDebugLevel==73){ - if(iSec<36){ - fHistoRow->Fill(iRow); - fHistoPad->Fill(iPad); - fHistoTime->Fill(iTimeBin); - fHistoRowPad->Fill(iPad,iRow); - } - }else if(fDebugLevel==74){ - if(iSec>=36){ - fHistoRow->Fill(iRow); - fHistoPad->Fill(iPad); - fHistoTime->Fill(iTimeBin); - fHistoRowPad->Fill(iPad,iRow); - } - } - - //check row consistency - if (iRow < 0 ) continue; - if (iRow < 0 || iRow >= nRows){ - AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !", - iRow, 0, nRows -1)); - continue; - } - - //check pad consistency - if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) { - AliError(Form("Pad index (%d) outside the range (%d -> %d) !", - iPad, 0, roc->GetNPads(iSec,iRow) )); - continue; - } - - //check time consistency - if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){ - //cout< %d) !", - iTimeBin, 0, fRecoParam->GetLastBin() -1)); - } - - //signal - Int_t signal = input.GetSignal(); - if (signal <= fZeroSup || - iTimeBin < fFirstBin || - iTimeBin > fLastBin - ) { - digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad); - continue; - } - if (!noiseROC) continue; - Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector - if (noiseOnPad > fMaxNoiseAbs){ - digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad); - continue; // consider noisy pad as dead - } - if(signal <= fMaxNoiseSigma * noiseOnPad){ - digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad); - continue; - } - digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad); - }//end of loop over altro data - }//end of loop over sectors - - SetDigArr(digarr); - if(isAltro) FindClusterKrIO(); - delete digarr; - - return 0; -} - - -Int_t AliTPCclustererKr::CleanSector(Int_t sector){ - // - // clean isolated digits - // - const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector - for(Int_t iRow=0; iRowGetRow(sector,iRow);//real data - }else{ - digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC - } - if(!digrow) continue; - digrow->ExpandBuffer(); //decrunch - const Int_t kNPads = digrow->GetNCols(); // number of pads - const Int_t kNTime = digrow->GetNRows(); // number of timebins - for(Int_t iPad=1;iPadGetDigitsColumn(iPad); - - for(Int_t iTimeBin=1;iTimeBinGetNSector();//number of sectors - for(Int_t iSec=0; iSec maximaInSector; - TObjArray *maximaInSector=new TObjArray();//to store AliPadMax* - - // - // looking for the maxima on the pad - // - - - - const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector - for(Int_t iRow=0; iRowGetRow(iSec,iRow);//real data - }else{ - digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC - } - if(digrow){//if pointer exist - digrow->ExpandBuffer(); //decrunch - const Int_t kNPads = digrow->GetNCols(); // number of pads - const Int_t kNTime = digrow->GetNRows(); // number of timebins - for(Int_t iPad=0;iPadGetDigitsColumn(iPad); - for(Int_t iTimeBin=1;iTimeBinAddLast(oneMaximum); - - timeBinMax=-1; - valueMaximum=-1; - increaseBegin=-1; - sumAdc=0; - ifIncreaseBegin=true; - ifMaximum=false; - } - continue; - } - - if(ifIncreaseBegin){ - ifIncreaseBegin=false; - increaseBegin=iTimeBin; - } - - if(adc>valueMaximum){ - timeBinMax=iTimeBin; - valueMaximum=adc; - ifMaximum=true; - } - sumAdc+=adc; - if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge - //at least 3 timebins - //insert maximum, default values and set flags - //Double_t xCord,yCord; - //GetXY(iSec,iRow,iPad,xCord,yCord); - Double_t x[]={iRow,iPad,iTimeBin}; - Int_t i[]={iSec}; - AliTPCTransform trafo; - trafo.Transform(x,i,0,1); - AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum, - timeBinMax, - iPad, - iRow, - x[0],//xCord, - x[1],//yCord, - x[2]/*timeBinMax*/), - increaseBegin, - iTimeBin-1, - sumAdc); - maximaInSector->AddLast(oneMaximum); - - timeBinMax=-1; - valueMaximum=-1; - increaseBegin=-1; - sumAdc=0; - ifIncreaseBegin=true; - ifMaximum=false; - continue; - } - - }//end loop over timebins - }//end loop over pads -// }else{ -// cout<<"Pointer does not exist!!"<SetOwner(kTRUE); - maximaInSector->Delete(); - delete maximaInSector; - }//end sector for - cout<<"Number of clusters in event: "<GetEntriesFast(); - for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) { - - AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1); - if (!mp1) continue; - AliTPCclusterKr clusterKr; - - Int_t nUsedPads=1; - Int_t clusterValue=0; - clusterValue+=(mp1)->GetSum(); - list nUsedRows; - nUsedRows.push_back((mp1)->GetRow()); - - maxDig =(mp1)->GetAdc() ; - maxSumAdc =(mp1)->GetSum() ; - maxTimeBin =(mp1)->GetTime(); - maxPad =(mp1)->GetPad() ; - maxRow =(mp1)->GetRow() ; - maxX =(mp1)->GetX(); - maxY =(mp1)->GetY(); - maxT =(mp1)->GetT(); - - AliSimDigits *digrowTmp; - if(fRawData){ - digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow()); - }else{ - digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow()); - } - - digrowTmp->ExpandBuffer(); //decrunch - - for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){ - Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad()); - AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT()); - clusterKr.AddDigitToCluster(vtpr); - } - clusterKr.SetCenter();//set centr of the cluster - - for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) { - AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2); - if (!mp2) continue; - if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX())>fMaxPadRangeCm) continue; - if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue; - if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) >fMaxTimeRange) continue; - - { - clusterValue+=(mp2)->GetSum(); - - nUsedPads++; - nUsedRows.push_back((mp2)->GetRow()); - - AliSimDigits *digrowTmp1; - if(fRawData){ - digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow()); - }else{ - digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow()); - } - digrowTmp1->ExpandBuffer(); //decrunch - - for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){ - Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad()); - AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT()); - clusterKr.AddDigitToCluster(vtpr); - } - - clusterKr.SetCenter();//set center of the cluster - - //which one is bigger - if( (mp2)->GetAdc() > maxDig ){ - maxDig =(mp2)->GetAdc() ; - maxSumAdc =(mp2)->GetSum() ; - maxTimeBin =(mp2)->GetTime(); - maxPad =(mp2)->GetPad() ; - maxRow =(mp2)->GetRow() ; - maxX =(mp2)->GetX() ; - maxY =(mp2)->GetY() ; - maxT =(mp2)->GetT() ; - } else if ( (mp2)->GetAdc() == maxDig ){ - if( (mp2)->GetSum() > maxSumAdc){ - maxDig =(mp2)->GetAdc() ; - maxSumAdc =(mp2)->GetSum() ; - maxTimeBin =(mp2)->GetTime(); - maxPad =(mp2)->GetPad() ; - maxRow =(mp2)->GetRow() ; - maxX =(mp2)->GetX() ; - maxY =(mp2)->GetY() ; - maxT =(mp2)->GetT() ; - } - } - delete maximaInSector->RemoveAt(it2); - } - }//inside loop - delete maximaInSector->RemoveAt(it1); - clusterKr.SetSize(); - //through out clusters on the edge and noise - //if(clusterValue/clusterKr.fCluster.size()GetPadRowRadii(sec,row);//radius of row in sector in cm - - Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row - Float_t padXSize; - if(secGetNInnerSector())padXSize=0.4; - else padXSize=0.6; - Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad - - Float_t sin,cos; - fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector - - Double_t xGlob1 = xLocal * cos + yLocal * sin; - Double_t yGlob1 = -xLocal * sin + yLocal * cos; - - - Double_t rot=0; - rot=TMath::Pi()/2.; - - xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot); - yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot); - - yGlob=-1*yGlob; - if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob; - - - return; -} +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id: AliTPCclustererKr.cxx,v 1.7 2008/02/06 17:24:53 matyja Exp $ */ + +//----------------------------------------------------------------- +// Implementation of the TPC Kr cluster class +// +// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl +//----------------------------------------------------------------- + +/* +Instruction - how to use that: +There are two macros prepared. One is for preparing clusters from MC +samples: FindKrClusters.C. The output is kept in TPC.RecPoints.root. +The other macro is prepared for data analysis: FindKrClustersRaw.C. +The output is created for each processed file in root file named adc.root. +For each data subsample the same named file is created. So be careful +do not overwrite them. + +Additional selection criteria to select the GOLD cluster +Example: +// open file with clusters +TFile f("Krypton.root"); +TTree * tree = (TTree*)f.Get("Kr") +TCut cutR0("cutR0","fADCcluster/fSize<100"); // adjust it according v seetings - +TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal +TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.2"); // digital noise removal +TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal +TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks +TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutS1; +This values are typical values to be applied in selectors + + +* +**** MC **** +* + +To run clusterizaton for MC type: +.x FindKrClusters.C + +If you don't want to use the standard selection criteria then you +have to do following: + +// load the standard setup +AliRunLoader* rl = AliRunLoader::Open("galice.root"); +AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader"); +tpcl->LoadDigits(); +rl->LoadgAlice(); +gAlice=rl->GetAliRun(); +TDirectory *cwd = gDirectory; +AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC"); +Int_t ver = tpc->IsVersion(); +rl->CdGAFile(); +AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60"); +AliTPCDigitsArray *digarr=new AliTPCDigitsArray; +digarr->Setup(param); +cwd->cd(); + +//loop over events +Int_t nevmax=rl->GetNumberOfEvents(); +for(Int_t nev=0;nevGetEvent(nev); + TTree* input_tree= tpcl->TreeD();//load tree with digits + digarr->ConnectTree(input_tree); + TTree *output_tree =tpcl->TreeR();//load output tree + + AliTPCclustererKr *clusters = new AliTPCclustererKr(); + clusters->SetParam(param); + clusters->SetInput(input_tree); + clusters->SetOutput(output_tree); + clusters->SetDigArr(digarr); + +//If you want to change the cluster finder parameters for MC there are +//several of them: + +//1. signal threshold (everything below the given number is treated as 0) + clusters->SetMinAdc(3); + +//2. number of neighbouring timebins to be considered + clusters->SetMinTimeBins(2); + +//3. distance of the cluster center to the center of a pad in pad-padrow plane +//(in cm). Remenber that this is still quantified by pad size. + clusters->SetMaxPadRangeCm(2.5); + +//4. distance of the cluster center to the center of a padrow in pad-padrow +//plane (in cm). Remenber that this is still quantified by pad size. + clusters->SetMaxRowRangeCm(3.5); + +//5. distance of the cluster center to the max time bin on a pad (in tackts) +//ie. fabs(centerT - time)<7 + clusters->SetMaxTimeRange(7); + +//6. cut reduce peak at 0. There are noises which appear mostly as two +//timebins on one pad. + clusters->SetValueToSize(3.5); + + + clusters->finderIO(); + tpcl->WriteRecPoints("OVERWRITE"); +} +delete rl;//cleans everything + +* +********* DATA ********* +* + +To run clusterizaton for DATA for file named raw_data.root type: +.x FindKrClustersRaw.C("raw_data.root") + +If you want to change some criteria do the following: + +// +// remove Altro warnings +// +AliLog::SetClassDebugLevel("AliTPCRawStream",-5); +AliLog::SetClassDebugLevel("AliRawReaderDate",-5); +AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5); +AliLog::SetModuleDebugLevel("RAW",-5); + +// +// Get database with noises +// +// char *ocdbpath = gSystem->Getenv("OCDB_PATH"); +char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB"; +if (ocdbpath==0){ +ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/"; +} +AliCDBManager * man = AliCDBManager::Instance(); +man->SetDefaultStorage(ocdbpath); +man->SetRun(0); +AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); +AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); + +//define tree +TFile *hfile=new TFile("adc.root","RECREATE","ADC file"); +// Create a ROOT Tree +TTree *mytree = new TTree("Kr","Krypton cluster tree"); + +//define infput file +const char *fileName="data.root"; +AliRawReader *reader = new AliRawReaderRoot(fileName); +//AliRawReader *reader = new AliRawReaderDate(fileName); +reader->Reset(); +AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader); +stream->SelectRawData("TPC"); + +//one general output +AliTPCclustererKr *clusters = new AliTPCclustererKr(); +clusters->SetOutput(mytree); +clusters->SetRecoParam(0);//standard reco parameters +AliTPCParamSR *param=new AliTPCParamSR(); +clusters->SetParam(param);//TPC parameters(sectors, timebins, etc.) + +//set cluster finder parameters (from data): +//1. zero suppression parameter + clusters->SetZeroSup(param->GetZeroSup()); + +//2. first bin + clusters->SetFirstBin(60); + +//3. last bin + clusters->SetLastBin(950); + +//4. maximal noise + clusters->SetMaxNoiseAbs(2); + +//5. maximal amount of sigma of noise + clusters->SetMaxNoiseSigma(3); + +//The remaining parameters are the same paramters as for MC (see MC section +//points 1-6) + clusters->SetMinAdc(3); + clusters->SetMinTimeBins(2); + clusters->SetMaxPadRangeCm(2.5); + clusters->SetMaxRowRangeCm(3.5); + clusters->SetMaxTimeRange(7); + clusters->SetValueToSize(3.5); + +while (reader->NextEvent()) { + clusters->FinderIO(reader); +} + +hfile->Write(); +hfile->Close(); +delete stream; + + +*/ + +#include "AliTPCclustererKr.h" +#include "AliTPCclusterKr.h" +//#include +#include +#include "TObject.h" +#include "AliPadMax.h" +#include "AliSimDigits.h" +#include "AliTPCv4.h" +#include "AliTPCParam.h" +#include "AliTPCDigitsArray.h" +#include "AliTPCvtpr.h" +#include "AliTPCClustersRow.h" +#include "TTree.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TTreeStream.h" + +#include "AliTPCTransform.h" + +//used in raw data finder +#include "AliTPCROC.h" +#include "AliTPCCalPad.h" +#include "AliTPCAltroMapping.h" +#include "AliTPCcalibDB.h" +#include "AliTPCRawStream.h" +#include "AliTPCRecoParam.h" +#include "AliTPCReconstructor.h" +#include "AliRawReader.h" +#include "AliTPCCalROC.h" + +ClassImp(AliTPCclustererKr) + + +AliTPCclustererKr::AliTPCclustererKr() + :TObject(), + fRawData(kFALSE), + fInput(0), + fOutput(0), + fParam(0), + fDigarr(0), + fRecoParam(0), + fZeroSup(2), + fFirstBin(60), + fLastBin(950), + fMaxNoiseAbs(2), + fMaxNoiseSigma(3), + fMinAdc(3), + fMinTimeBins(2), +// fMaxPadRange(4), +// fMaxRowRange(3), + fMaxTimeRange(7), + fValueToSize(3.5), + fMaxPadRangeCm(2.5), + fMaxRowRangeCm(3.5), + fIsolCut(3), + fDebugLevel(-1), + fHistoRow(0), + fHistoPad(0), + fHistoTime(0), + fHistoRowPad(0) +{ +// +// default constructor +// +} + +AliTPCclustererKr::AliTPCclustererKr(const AliTPCclustererKr ¶m) + :TObject(), + fRawData(kFALSE), + fInput(0), + fOutput(0), + fParam(0), + fDigarr(0), + fRecoParam(0), + fZeroSup(2), + fFirstBin(60), + fLastBin(950), + fMaxNoiseAbs(2), + fMaxNoiseSigma(3), + fMinAdc(3), + fMinTimeBins(2), +// fMaxPadRange(4), +// fMaxRowRange(3), + fMaxTimeRange(7), + fValueToSize(3.5), + fMaxPadRangeCm(2.5), + fMaxRowRangeCm(3.5), + fIsolCut(3), + fDebugLevel(-1), + fHistoRow(0), + fHistoPad(0), + fHistoTime(0), + fHistoRowPad(0) +{ +// +// copy constructor +// + fParam = param.fParam; + fRecoParam = param.fRecoParam; + fRawData = param.fRawData; + fInput = param.fInput ; + fOutput = param.fOutput; + fDigarr = param.fDigarr; + fZeroSup = param.fZeroSup ; + fFirstBin = param.fFirstBin ; + fLastBin = param.fLastBin ; + fMaxNoiseAbs = param.fMaxNoiseAbs ; + fMaxNoiseSigma = param.fMaxNoiseSigma ; + fMinAdc = param.fMinAdc; + fMinTimeBins = param.fMinTimeBins; +// fMaxPadRange = param.fMaxPadRange ; +// fMaxRowRange = param.fMaxRowRange ; + fMaxTimeRange = param.fMaxTimeRange; + fValueToSize = param.fValueToSize; + fMaxPadRangeCm = param.fMaxPadRangeCm; + fMaxRowRangeCm = param.fMaxRowRangeCm; + fIsolCut = param.fIsolCut; + fDebugLevel = param.fDebugLevel; + fHistoRow = param.fHistoRow ; + fHistoPad = param.fHistoPad ; + fHistoTime = param.fHistoTime; + fHistoRowPad = param.fHistoRowPad; + + +} + +AliTPCclustererKr & AliTPCclustererKr::operator = (const AliTPCclustererKr & param) +{ + // + // assignment operator + // + fParam = param.fParam; + fRecoParam = param.fRecoParam; + fRawData = param.fRawData; + fInput = param.fInput ; + fOutput = param.fOutput; + fDigarr = param.fDigarr; + fZeroSup = param.fZeroSup ; + fFirstBin = param.fFirstBin ; + fLastBin = param.fLastBin ; + fMaxNoiseAbs = param.fMaxNoiseAbs ; + fMaxNoiseSigma = param.fMaxNoiseSigma ; + fMinAdc = param.fMinAdc; + fMinTimeBins = param.fMinTimeBins; +// fMaxPadRange = param.fMaxPadRange ; +// fMaxRowRange = param.fMaxRowRange ; + fMaxTimeRange = param.fMaxTimeRange; + fValueToSize = param.fValueToSize; + fMaxPadRangeCm = param.fMaxPadRangeCm; + fMaxRowRangeCm = param.fMaxRowRangeCm; + fIsolCut = param.fIsolCut; + fDebugLevel = param.fDebugLevel; + fHistoRow = param.fHistoRow ; + fHistoPad = param.fHistoPad ; + fHistoTime = param.fHistoTime; + fHistoRowPad = param.fHistoRowPad; + return (*this); +} + +AliTPCclustererKr::~AliTPCclustererKr() +{ + // + // destructor + // + delete fOutput; +} + +void AliTPCclustererKr::SetRecoParam(AliTPCRecoParam *recoParam) +{ + // + // set reconstruction parameters + // + if (recoParam) { + fRecoParam = recoParam; + }else{ + //set default parameters if not specified + fRecoParam = AliTPCReconstructor::GetRecoParam(); + if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam(); + } + return; +} + + +////____________________________________________________________________________ +//// I/O +void AliTPCclustererKr::SetInput(TTree * tree) +{ + // + // set input tree with digits + // + fInput = tree; + if (!fInput->GetBranch("Segment")){ + cerr<<"AliTPCclusterKr::FindClusterKr(): no proper input tree !\n"; + fInput=0; + return; + } +} + +void AliTPCclustererKr::SetOutput(TTree * /*tree*/) +{ + // + //dummy + // + fOutput = new TTreeSRedirector("Krypton.root"); +} + +////____________________________________________________________________________ +//// with new I/O +Int_t AliTPCclustererKr::FinderIO() +{ + // Krypton cluster finder for simulated events from MC + + if (!fInput) { + Error("Digits2Clusters", "input tree not initialised"); + return 10; + } + + if (!fOutput) { + Error("Digits2Clusters", "output tree not initialised"); + return 11; + } + + FindClusterKrIO(); + return 0; +} + +Int_t AliTPCclustererKr::FinderIO(AliRawReader* rawReader) +{ + // Krypton cluster finder for the TPC raw data + // + // fParam must be defined before + + if(rawReader)fRawData=kTRUE; //set flag to data + + if (!fOutput) { + Error("Digits2Clusters", "output tree not initialised"); + return 11; + } + + fParam->SetMaxTBin(fRecoParam->GetLastBin());//set number of timebins from reco -> param + // used later for memory allocation + + Bool_t isAltro=kFALSE; + + AliTPCROC * roc = AliTPCROC::Instance(); + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); + AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); + // + AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping); + + const Int_t kNIS = fParam->GetNInnerSector();//number of inner sectors + const Int_t kNOS = fParam->GetNOuterSector();//number of outer sectors + const Int_t kNS = kNIS + kNOS;//all sectors + + + //crate TPC view + AliTPCDigitsArray *digarr=new AliTPCDigitsArray(kFALSE);//data not sim + digarr->Setup(fParam);//as usually parameters + + // + // Loop over sectors + // + for(Int_t iSec = 0; iSec < kNS; iSec++) { + AliTPCCalROC * noiseROC; + if(noiseTPC==0x0){ + noiseROC =new AliTPCCalROC(iSec);//noise=0 + } + else{ + noiseROC = noiseTPC->GetCalROC(iSec); // noise per given sector + } + Int_t nRows = 0; //number of rows in sector + Int_t nDDLs = 0; //number of DDLs + Int_t indexDDL = 0; //DDL index + if (iSec < kNIS) { + nRows = fParam->GetNRowLow(); + nDDLs = 2; + indexDDL = iSec * 2; + }else { + nRows = fParam->GetNRowUp(); + nDDLs = 4; + indexDDL = (iSec-kNIS) * 4 + kNIS * 2; + } + + // + // Load the raw data for corresponding DDLs + // + rawReader->Reset(); + rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); + + if(input.Next()) { + isAltro=kTRUE; + // Allocate memory for rows in sector (pads(depends on row) x timebins) + for(Int_t iRow = 0; iRow < nRows; iRow++) { + digarr->CreateRow(iSec,iRow); + }//end loop over rows + } + rawReader->Reset(); + rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); + + // + // Begin loop over altro data + // + while (input.Next()) { + + //check sector consistency + if (input.GetSector() != iSec) + AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSec,input.GetSector())); + + Int_t iRow = input.GetRow(); + Int_t iPad = input.GetPad(); + Int_t iTimeBin = input.GetTime(); + + // + if(fDebugLevel==72){ + fHistoRow->Fill(iRow); + fHistoPad->Fill(iPad); + fHistoTime->Fill(iTimeBin); + fHistoRowPad->Fill(iPad,iRow); + }else if(fDebugLevel>=0&&fDebugLevel<72){ + if(iSec==fDebugLevel){ + fHistoRow->Fill(iRow); + fHistoPad->Fill(iPad); + fHistoTime->Fill(iTimeBin); + fHistoRowPad->Fill(iPad,iRow); + } + }else if(fDebugLevel==73){ + if(iSec<36){ + fHistoRow->Fill(iRow); + fHistoPad->Fill(iPad); + fHistoTime->Fill(iTimeBin); + fHistoRowPad->Fill(iPad,iRow); + } + }else if(fDebugLevel==74){ + if(iSec>=36){ + fHistoRow->Fill(iRow); + fHistoPad->Fill(iPad); + fHistoTime->Fill(iTimeBin); + fHistoRowPad->Fill(iPad,iRow); + } + } + + //check row consistency + if (iRow < 0 ) continue; + if (iRow < 0 || iRow >= nRows){ + AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !", + iRow, 0, nRows -1)); + continue; + } + + //check pad consistency + if (iPad < 0 || iPad >= (Int_t)(roc->GetNPads(iSec,iRow))) { + AliError(Form("Pad index (%d) outside the range (%d -> %d) !", + iPad, 0, roc->GetNPads(iSec,iRow) )); + continue; + } + + //check time consistency + if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){ + //cout< %d) !", + iTimeBin, 0, fRecoParam->GetLastBin() -1)); + } + + //signal + Int_t signal = input.GetSignal(); + if (signal <= fZeroSup || + iTimeBin < fFirstBin || + iTimeBin > fLastBin + ) { + digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad); + continue; + } + if (!noiseROC) continue; + Double_t noiseOnPad = noiseROC->GetValue(iRow,iPad);//noise on given pad and row in sector + if (noiseOnPad > fMaxNoiseAbs){ + digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad); + continue; // consider noisy pad as dead + } + if(signal <= fMaxNoiseSigma * noiseOnPad){ + digarr->GetRow(iSec,iRow)->SetDigitFast(0,iTimeBin,iPad); + continue; + } + digarr->GetRow(iSec,iRow)->SetDigitFast(signal,iTimeBin,iPad); + }//end of loop over altro data + }//end of loop over sectors + + SetDigArr(digarr); + if(isAltro) FindClusterKrIO(); + delete digarr; + + return 0; +} + +void AliTPCclustererKr::CleanSector(Int_t sector){ + // + // clean isolated digits + // + const Int_t kNRows=fParam->GetNRow(sector);//number of rows in sector + for(Int_t iRow=0; iRowGetRow(sector,iRow);//real data + }else{ + digrow = (AliSimDigits*)fDigarr->LoadRow(sector,iRow);//MC + } + if(!digrow) continue; + digrow->ExpandBuffer(); //decrunch + const Int_t kNPads = digrow->GetNCols(); // number of pads + const Int_t kNTime = digrow->GetNRows(); // number of timebins + for(Int_t iPad=1;iPadGetDigitsColumn(iPad); + + for(Int_t iTimeBin=1;iTimeBinGetNSector();//number of sectors + for(Int_t iSec=0; iSec maximaInSector; + TObjArray *maximaInSector=new TObjArray();//to store AliPadMax* + + // + // looking for the maxima on the pad + // + + const Int_t kNRows=fParam->GetNRow(iSec);//number of rows in sector + for(Int_t iRow=0; iRowGetRow(iSec,iRow);//real data + }else{ + digrow = (AliSimDigits*)fDigarr->LoadRow(iSec,iRow);//MC + } + if(digrow){//if pointer exist + digrow->ExpandBuffer(); //decrunch + const Int_t kNPads = digrow->GetNCols(); // number of pads + const Int_t kNTime = digrow->GetNRows(); // number of timebins + for(Int_t iPad=0;iPadGetDigitsColumn(iPad); + for(Int_t iTimeBin=1;iTimeBinAddLast(oneMaximum); + + timeBinMax=-1; + valueMaximum=-1; + increaseBegin=-1; + sumAdc=0; + ifIncreaseBegin=true; + ifMaximum=false; + } + continue; + } + + + + + + + if(ifIncreaseBegin){ + ifIncreaseBegin=false; + increaseBegin=iTimeBin; + } + + if(adc>valueMaximum){ + timeBinMax=iTimeBin; + valueMaximum=adc; + ifMaximum=true; + } + sumAdc+=adc; + if(iTimeBin==kNTime-1 && ifMaximum && kNTime-increaseBegin>fMinTimeBins){//on the edge + //at least 3 timebins + //insert maximum, default values and set flags + //Double_t xCord,yCord; + //GetXY(iSec,iRow,iPad,xCord,yCord); + Double_t x[]={iRow,iPad,iTimeBin}; + Int_t i[]={iSec}; + AliTPCTransform trafo; + trafo.Transform(x,i,0,1); + AliPadMax *oneMaximum = new AliPadMax(AliTPCvtpr(valueMaximum, + timeBinMax, + iPad, + iRow, + x[0],//xCord, + x[1],//yCord, + x[2]/*timeBinMax*/), + increaseBegin, + iTimeBin-1, + sumAdc); + maximaInSector->AddLast(oneMaximum); + + timeBinMax=-1; + valueMaximum=-1; + increaseBegin=-1; + sumAdc=0; + ifIncreaseBegin=true; + ifMaximum=false; + continue; + } + + }//end loop over timebins + }//end loop over pads +// }else{ +// cout<<"Pointer does not exist!!"<SetOwner(kTRUE); + maximaInSector->Delete(); + delete maximaInSector; + }//end sector for + cout<<"Number of clusters in event: "<GetEntriesFast(); + for(Int_t it1 = 0; it1 < entriesArr; ++it1 ) { + + AliPadMax *mp1=(AliPadMax *)maximaInSector->UncheckedAt(it1); + if (!mp1) continue; + AliTPCclusterKr clusterKr; + + Int_t nUsedPads=1; + Int_t clusterValue=0; + clusterValue+=(mp1)->GetSum(); + list nUsedRows; + nUsedRows.push_back((mp1)->GetRow()); + + maxDig =(mp1)->GetAdc() ; + maxSumAdc =(mp1)->GetSum() ; + maxTimeBin =(mp1)->GetTime(); + maxPad =(mp1)->GetPad() ; + maxRow =(mp1)->GetRow() ; + maxX =(mp1)->GetX(); + maxY =(mp1)->GetY(); + maxT =(mp1)->GetT(); + + AliSimDigits *digrowTmp; + if(fRawData){ + digrowTmp = (AliSimDigits*)fDigarr->GetRow(iSec,(mp1)->GetRow()); + }else{ + digrowTmp = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp1)->GetRow()); + } + + digrowTmp->ExpandBuffer(); //decrunch + + for(Int_t itb=(mp1)->GetBegin(); itb<((mp1)->GetEnd())+1; itb++){ + Int_t adcTmp = digrowTmp->GetDigitUnchecked(itb,(mp1)->GetPad()); + AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp1)->GetPad(),(mp1)->GetRow(),(mp1)->GetX(),(mp1)->GetY(),(mp1)->GetT()); + clusterKr.AddDigitToCluster(vtpr); + } + clusterKr.SetCenter();//set centr of the cluster + + for(Int_t it2 = it1+1; it2 < entriesArr; ++it2 ) { + AliPadMax *mp2=(AliPadMax *)maximaInSector->UncheckedAt(it2); + if (!mp2) continue; + if (TMath::Abs(clusterKr.GetCenterX() - (mp2)->GetX()) > fMaxPadRangeCm) continue; + if (TMath::Abs(clusterKr.GetCenterY() - (mp2)->GetY()) > fMaxRowRangeCm) continue; + if (TMath::Abs(clusterKr.GetCenterT() - (mp2)->GetT()) > fMaxTimeRange) continue; + + { + clusterValue+=(mp2)->GetSum(); + + nUsedPads++; + nUsedRows.push_back((mp2)->GetRow()); + + AliSimDigits *digrowTmp1; + if(fRawData){ + digrowTmp1 = (AliSimDigits*)fDigarr->GetRow(iSec,(mp2)->GetRow()); + }else{ + digrowTmp1 = (AliSimDigits*)fDigarr->LoadRow(iSec,(mp2)->GetRow()); + } + digrowTmp1->ExpandBuffer(); //decrunch + + for(Int_t itb=(mp2)->GetBegin(); itb<(mp2)->GetEnd()+1; itb++){ + Int_t adcTmp = digrowTmp1->GetDigitUnchecked(itb,(mp2)->GetPad()); + AliTPCvtpr *vtpr=new AliTPCvtpr(adcTmp,itb,(mp2)->GetPad(),(mp2)->GetRow(),(mp2)->GetX(),(mp2)->GetY(),(mp2)->GetT()); + clusterKr.AddDigitToCluster(vtpr); + } + + clusterKr.SetCenter();//set center of the cluster + + //which one is bigger + if( (mp2)->GetAdc() > maxDig ){ + maxDig =(mp2)->GetAdc() ; + maxSumAdc =(mp2)->GetSum() ; + maxTimeBin =(mp2)->GetTime(); + maxPad =(mp2)->GetPad() ; + maxRow =(mp2)->GetRow() ; + maxX =(mp2)->GetX() ; + maxY =(mp2)->GetY() ; + maxT =(mp2)->GetT() ; + } else if ( (mp2)->GetAdc() == maxDig ){ + if( (mp2)->GetSum() > maxSumAdc){ + maxDig =(mp2)->GetAdc() ; + maxSumAdc =(mp2)->GetSum() ; + maxTimeBin =(mp2)->GetTime(); + maxPad =(mp2)->GetPad() ; + maxRow =(mp2)->GetRow() ; + maxX =(mp2)->GetX() ; + maxY =(mp2)->GetY() ; + maxT =(mp2)->GetT() ; + } + } + delete maximaInSector->RemoveAt(it2); + } + }//inside loop + delete maximaInSector->RemoveAt(it1); + clusterKr.SetSize(); + //through out clusters on the edge and noise + //if(clusterValue/clusterKr.fCluster.size()GetPadRowRadii(sec,row);//radius of row in sector in cm + + Int_t padmax = fParam->GetNPads(sec,row);//number of pads in a given row + Float_t padXSize; + if(secGetNInnerSector())padXSize=0.4; + else padXSize=0.6; + Double_t xLocal=(pad+0.5-padmax/2.)*padXSize;//x-value of the center of pad + + Float_t sin,cos; + fParam->AdjustCosSin((Int_t)sec,cos,sin);//return sinus and cosinus of the sector + + Double_t xGlob1 = xLocal * cos + yLocal * sin; + Double_t yGlob1 = -xLocal * sin + yLocal * cos; + + + Double_t rot=0; + rot=TMath::Pi()/2.; + + xGlob = xGlob1 * TMath::Cos(rot) + yGlob1 * TMath::Sin(rot); + yGlob = -xGlob1 * TMath::Sin(rot) + yGlob1 * TMath::Cos(rot); + + yGlob=-1*yGlob; + if(sec<18||(sec>=36&&sec<54)) xGlob =-1*xGlob; + + + return; +} diff --git a/TPC/AliTPCclustererKr.h b/TPC/AliTPCclustererKr.h index 86fce7809b0..46a01d10e04 100644 --- a/TPC/AliTPCclustererKr.h +++ b/TPC/AliTPCclustererKr.h @@ -1,160 +1,165 @@ -#ifndef ALITPCCLUSTERERKR_H -#define ALITPCCLUSTERERKR_H -/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -/* $Id: AliTPCclusterKr.h,v 1.8 2008/02/07 16:07:15 matyja Exp $ */ - -//------------------------------------------------------- -// TPC Kr Clusterer Class -// -// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl -//------------------------------------------------------- - -#include "TObject.h" - -class AliTPCclusterKr; -class AliPadMax; -class AliSimDigits; -class AliTPCv4; -class AliTPCParam; -class AliTPCDigitsArray; -class AliTPCvtpr; -class AliTPCClustersRow; -class TTree; -class TH1F; -class TH2F; - -class AliTPCTransform; - -//used in raw data finder -class AliTPCROC; -class AliTPCCalPad; -class AliTPCAltroMapping; -class AliTPCcalibDB; -class AliTPCRawStream; -class AliTPCRecoParam; -class AliTPCReconstructor; -class AliRawReader; -class AliTPCCalROC; -class TTreeSRedirector; -//_____________________________________________________________________________ -class AliTPCclustererKr: public TObject{ -public: - AliTPCclustererKr(); - AliTPCclustererKr(const AliTPCclustererKr ¶m);//copy constructor - AliTPCclustererKr &operator = (const AliTPCclustererKr & param); - virtual ~AliTPCclustererKr(); - - //finders - virtual Int_t FinderIO();//for MC - virtual Int_t FinderIO(AliRawReader* rawReader);//for data - virtual Int_t FindClusterKrIO();//main routine for finding clusters - virtual Int_t CleanSector(Int_t sector); // clean isolated digits - - - //other - void GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob);//give XY coordinate of the pad - - virtual void SetInput(TTree * tree); //set input tree with digits - virtual void SetOutput(TTree * tree); //set output tree with clusters - virtual void SetParam(AliTPCParam *param){fParam=param;}//set TPC parameters - virtual void SetDigArr(AliTPCDigitsArray *digarr){fDigarr=digarr;}//set current array of digits - virtual void SetRecoParam(AliTPCRecoParam *recoParam=0);//set reconstruction parameters - - - - //setters for cluster finder parameters - virtual void SetZeroSup(Int_t v){fZeroSup=v;}//set zero suppresion parameter - virtual void SetFirstBin(Int_t v){fFirstBin=v;}//set first considered timebin - virtual void SetLastBin(Int_t v){fLastBin=v;}//set last considered timebin - virtual void SetMaxNoiseAbs(Float_t v){fMaxNoiseAbs=v;}//set maximal noise value - virtual void SetMaxNoiseSigma(Float_t v){fMaxNoiseSigma=v;}//set maximal noise sigma - - virtual void SetMinAdc(Int_t v){v<=0?fMinAdc=1:fMinAdc=v;}//set fMinAdc - virtual void SetMinTimeBins(Int_t v){fMinTimeBins=v;}//set fMinTimeBins -// virtual void SetMaxPadRange(Int_t v){fMaxPadRange=v;}//set fMaxPadRange -// virtual void SetMaxRowRange(Int_t v){fMaxRowRange=v;}//set fMaxRowRange - virtual void SetMaxTimeRange(Int_t v){fMaxTimeRange=v;}//set fMaxTimeRange - virtual void SetValueToSize(Float_t v){fValueToSize=v;}//set fValueToSize - - virtual void SetMaxPadRangeCm(Double_t v){fMaxPadRangeCm=v;}//set fMaxPadRangeCm - virtual void SetMaxRowRangeCm(Double_t v){fMaxRowRangeCm=v;}//set fMaxRowRangeCm - - virtual void SetDebugLevel(Int_t debug){fDebugLevel=debug;} - //debug = 0 to 71 -sector number to print - // = 72 - all sectors - // = 73 - inners - // = 74 - outers - - virtual void SetHistoRow(TH1F *histo) {fHistoRow =histo;} - virtual void SetHistoPad(TH1F *histo) {fHistoPad =histo;} - virtual void SetHistoTime(TH1F *histo) {fHistoTime =histo;} - virtual void SetHistoRowPad(TH2F *histo){fHistoRowPad=histo;} - - //getters for cluster finder parameters - Int_t GetZeroSup() const {return fZeroSup;}//get zero suppresion parameter - Int_t GetFirstBin() const {return fFirstBin;}//get first considered timebin - Int_t GetLastBin() const {return fLastBin;}//get last considered timebin - Float_t GetMaxNoiseAbs() const {return fMaxNoiseAbs;}//get maximal noise value - Float_t GetMaxNoiseSigma() const {return fMaxNoiseSigma;}//get maximal noise sigma - - Int_t GetMinAdc() const {return fMinAdc;}//get fMinAdc - Int_t GetMinTimeBins() const {return fMinTimeBins;}//get fMinTimeBins -// Int_t GetMaxPadRange() const {return fMaxPadRange;}//get fMaxPadRange -// Int_t GetMaxRowRange() const {return fMaxRowRange;}//get fMaxRowRange - Int_t GetMaxTimeRange() const {return fMaxTimeRange;}//get fMaxTimeRange - Float_t GetValueToSize() const {return fValueToSize;}//get fValueToSize - - Double_t GetMaxPadRangeCm() const {return fMaxPadRangeCm;}//get fMaxPadRangeCm - Double_t GetMaxRowRangeCm() const {return fMaxRowRangeCm;}//get fMaxRowRangeCm - - Int_t GetDebugLevel() const {return fDebugLevel;} - TH1F * GetHistoRow(){return fHistoRow;} - TH1F * GetHistoPad(){return fHistoPad;} - TH1F * GetHistoTime(){return fHistoTime;} - TH2F * GetHistoRowPad(){return fHistoRowPad;} - -private: - void MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter); - Bool_t fRawData; //flag =0 for MC =1 for real data - - TTree * fInput; //!input tree with digits - object not owner - TTreeSRedirector * fOutput; //!output tree with clusters - object not owner - - AliTPCParam * fParam;//!TPC parameters - AliTPCDigitsArray *fDigarr;//! pointer to current array if digits - - //only for raw data :) - const AliTPCRecoParam * fRecoParam; //! reconstruction parameters - - //cluster finder parameters - Int_t fZeroSup;//zero suppresion parameter = 2 def. - Int_t fFirstBin;//first considered time bin used by cluster finder = 60 def. - Int_t fLastBin;//last considered time bin used by cluster finder = 950 def. - Float_t fMaxNoiseAbs;// maximal noise value on pad used in cluster finder = 2 def. - Float_t fMaxNoiseSigma;// maximal noise sigma on pad used in cluster finder = 3 def. - - Int_t fMinAdc;//minimal value of acd count in each timebin = 3 def. - Int_t fMinTimeBins;//minimal value of time bins one after each other = 2 def. -// Int_t fMaxPadRange;//maximal pad range from maximum = 4 def. -// Int_t fMaxRowRange;//maximal row range from maximum = 3 def. - Int_t fMaxTimeRange;//maximal time bin range from maximum = 7 def. - Float_t fValueToSize;//ratio cluster value to cluster size = 3.5 def. - - Double_t fMaxPadRangeCm;//maximal pad range in cm from maximum = 2.5cm def. - Double_t fMaxRowRangeCm;//maximal row range in cm from maximum = 3.5cm def. - - Int_t fDebugLevel;//! debug level variable - TH1F *fHistoRow;//!debug histo for rows - TH1F *fHistoPad;//!debug histo for pads - TH1F *fHistoTime;//!debug histo for timebins - TH2F *fHistoRowPad;//!debug histo for rows and pads - - ClassDef(AliTPCclustererKr,6) // Time Projection Chamber Kr clusters -}; - - -#endif - - +#ifndef ALITPCCLUSTERERKR_H +#define ALITPCCLUSTERERKR_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id: AliTPCclusterKr.h,v 1.8 2008/02/07 16:07:15 matyja Exp $ */ + +//------------------------------------------------------- +// TPC Kr Clusterer Class +// +// Origin: Adam Matyja, INP PAN, adam.matyja@ifj.edu.pl +//------------------------------------------------------- + +#include "TObject.h" + +class AliTPCclusterKr; +class AliPadMax; +class AliSimDigits; +class AliTPCv4; +class AliTPCParam; +class AliTPCDigitsArray; +class AliTPCvtpr; +class AliTPCClustersRow; +class TTree; +class TH1F; +class TH2F; + +class AliTPCTransform; + +//used in raw data finder +class AliTPCROC; +class AliTPCCalPad; +class AliTPCAltroMapping; +class AliTPCcalibDB; +class AliTPCRawStream; +class AliTPCRecoParam; +class AliTPCReconstructor; +class AliRawReader; +class AliTPCCalROC; +class TTreeSRedirector; +//_____________________________________________________________________________ +class AliTPCclustererKr: public TObject{ +public: + AliTPCclustererKr(); + AliTPCclustererKr(const AliTPCclustererKr ¶m);//copy constructor + AliTPCclustererKr &operator = (const AliTPCclustererKr & param); + virtual ~AliTPCclustererKr(); + + //finders + virtual Int_t FinderIO();//for MC + virtual Int_t FinderIO(AliRawReader* rawReader);//for data + virtual Int_t FindClusterKrIO();//main routine for finding clusters + virtual void CleanSector(Int_t sector); // clean isolated digits + + //other + void GetXY(Int_t sec,Int_t row,Int_t pad,Double_t& xGlob,Double_t& yGlob);//give XY coordinate of the pad + + virtual void SetInput(TTree * tree); //set input tree with digits + virtual void SetOutput(TTree * tree); //set output tree with clusters + virtual void SetParam(AliTPCParam *param){fParam=param;}//set TPC parameters + virtual void SetDigArr(AliTPCDigitsArray *digarr){fDigarr=digarr;}//set current array of digits + virtual void SetRecoParam(AliTPCRecoParam *recoParam=0);//set reconstruction parameters + + + + //setters for cluster finder parameters + virtual void SetZeroSup(Int_t v){fZeroSup=v;}//set zero suppresion parameter + virtual void SetFirstBin(Int_t v){fFirstBin=v;}//set first considered timebin + virtual void SetLastBin(Int_t v){fLastBin=v;}//set last considered timebin + virtual void SetMaxNoiseAbs(Float_t v){fMaxNoiseAbs=v;}//set maximal noise value + virtual void SetMaxNoiseSigma(Float_t v){fMaxNoiseSigma=v;}//set maximal noise sigma + + virtual void SetMinAdc(Int_t v){v<=0?fMinAdc=1:fMinAdc=v;}//set fMinAdc + virtual void SetMinTimeBins(Int_t v){fMinTimeBins=v;}//set fMinTimeBins +// virtual void SetMaxPadRange(Int_t v){fMaxPadRange=v;}//set fMaxPadRange +// virtual void SetMaxRowRange(Int_t v){fMaxRowRange=v;}//set fMaxRowRange + virtual void SetMaxTimeRange(Int_t v){fMaxTimeRange=v;}//set fMaxTimeRange + virtual void SetValueToSize(Float_t v){fValueToSize=v;}//set fValueToSize + + virtual void SetMaxPadRangeCm(Double_t v){fMaxPadRangeCm=v;}//set fMaxPadRangeCm + virtual void SetMaxRowRangeCm(Double_t v){fMaxRowRangeCm=v;}//set fMaxRowRangeCm + + virtual void SetIsolCut(Short_t v){fIsolCut=v;} + + virtual void SetDebugLevel(Int_t debug){fDebugLevel=debug;} + //debug = 0 to 71 -sector number to print + // = 72 - all sectors + // = 73 - inners + // = 74 - outers + + virtual void SetHistoRow(TH1F *histo) {fHistoRow =histo;} + virtual void SetHistoPad(TH1F *histo) {fHistoPad =histo;} + virtual void SetHistoTime(TH1F *histo) {fHistoTime =histo;} + virtual void SetHistoRowPad(TH2F *histo){fHistoRowPad=histo;} + + //getters for cluster finder parameters + Int_t GetZeroSup() const {return fZeroSup;}//get zero suppresion parameter + Int_t GetFirstBin() const {return fFirstBin;}//get first considered timebin + Int_t GetLastBin() const {return fLastBin;}//get last considered timebin + Float_t GetMaxNoiseAbs() const {return fMaxNoiseAbs;}//get maximal noise value + Float_t GetMaxNoiseSigma() const {return fMaxNoiseSigma;}//get maximal noise sigma + + Int_t GetMinAdc() const {return fMinAdc;}//get fMinAdc + Int_t GetMinTimeBins() const {return fMinTimeBins;}//get fMinTimeBins +// Int_t GetMaxPadRange() const {return fMaxPadRange;}//get fMaxPadRange +// Int_t GetMaxRowRange() const {return fMaxRowRange;}//get fMaxRowRange + Int_t GetMaxTimeRange() const {return fMaxTimeRange;}//get fMaxTimeRange + Float_t GetValueToSize() const {return fValueToSize;}//get fValueToSize + + Double_t GetMaxPadRangeCm() const {return fMaxPadRangeCm;}//get fMaxPadRangeCm + Double_t GetMaxRowRangeCm() const {return fMaxRowRangeCm;}//get fMaxRowRangeCm + + Short_t GetIsolCut() const {return fIsolCut;} + + Int_t GetDebugLevel() const {return fDebugLevel;} + TH1F * GetHistoRow(){return fHistoRow;} + TH1F * GetHistoPad(){return fHistoPad;} + TH1F * GetHistoTime(){return fHistoTime;} + TH2F * GetHistoRowPad(){return fHistoRowPad;} + +private: + void MakeClusters(TObjArray * maximaInSector, Int_t iSec, Int_t &clusterCounter); + Bool_t fRawData; //flag =0 for MC =1 for real data + + TTree * fInput; //!input tree with digits - object not owner + TTreeSRedirector * fOutput; //!output tree with clusters - object not owner + + AliTPCParam * fParam;//!TPC parameters + AliTPCDigitsArray *fDigarr;//! pointer to current array if digits + + //only for raw data :) + const AliTPCRecoParam * fRecoParam; //! reconstruction parameters + + //cluster finder parameters + Int_t fZeroSup;//zero suppresion parameter = 2 def. + Int_t fFirstBin;//first considered time bin used by cluster finder = 60 def. + Int_t fLastBin;//last considered time bin used by cluster finder = 950 def. + Float_t fMaxNoiseAbs;// maximal noise value on pad used in cluster finder = 2 def. + Float_t fMaxNoiseSigma;// maximal noise sigma on pad used in cluster finder = 3 def. + + Int_t fMinAdc;//minimal value of acd count in each timebin = 3 def. + Int_t fMinTimeBins;//minimal value of time bins one after each other = 2 def. +// Int_t fMaxPadRange;//maximal pad range from maximum = 4 def. +// Int_t fMaxRowRange;//maximal row range from maximum = 3 def. + Int_t fMaxTimeRange;//maximal time bin range from maximum = 7 def. + Float_t fValueToSize;//ratio cluster value to cluster size = 3.5 def. + + Double_t fMaxPadRangeCm;//maximal pad range in cm from maximum = 2.5cm def. + Double_t fMaxRowRangeCm;//maximal row range in cm from maximum = 3.5cm def. + + Short_t fIsolCut;//isolation cut in 3D = 5 def. + + Int_t fDebugLevel;//! debug level variable + TH1F *fHistoRow;//!debug histo for rows + TH1F *fHistoPad;//!debug histo for pads + TH1F *fHistoTime;//!debug histo for timebins + TH2F *fHistoRowPad;//!debug histo for rows and pads + + ClassDef(AliTPCclustererKr,6) // Time Projection Chamber Kr clusters +}; + + +#endif + + diff --git a/TPC/FindKrClusters.C b/TPC/FindKrClusters.C index d68596967e6..d00b9ccc07f 100644 --- a/TPC/FindKrClusters.C +++ b/TPC/FindKrClusters.C @@ -1,138 +1,135 @@ -/**************************************************************************** - * Origin: A.Matyja amatyja@cern.ch * - ****************************************************************************/ - -/* - - macro to create array of clusters from TPC digits - input files - galice.root - digits.root - file with digits - usualy use link to galice.root - - in splitted mode - neccesary to create link to proper file - - output file - TPC.RecPoints.root - - -// Warning - if cluster file AliTPCclusters.root already exist - macro exit and don't produce anything - - -*/ - - -#ifndef __CINT__ -#include -#include "AliRun.h" -#include "AliTPCv4.h" -#include "AliTPCParam.h" -#include "AliTPCclusterKr.h" -#include "AliTPCclustererKr.h" -#include "TFile.h" -#include "TStopwatch.h" -#include "TTree.h" -#endif - -Int_t FindKrClusters(){ - - // - //Load DataBase - // - char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB"; - //char *ocdbpath ="local:///home/matyja/baza/OCDB"; - if (ocdbpath==0){ - ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/"; - } - printf("OCDB PATH = %s\n",ocdbpath); - AliCDBManager * man = AliCDBManager::Instance(); - man->SetDefaultStorage(ocdbpath); - man->SetRun(0); - - AliRunLoader* rl = AliRunLoader::Open("galice.root"); - if (rl == 0x0) { - cerr<<"Can not open session"<GetLoader("TPCLoader"); - if (tpcl == 0x0) { - cerr<<"Can not get TPC Loader"<LoadDigits()) { - cerr<<"Error occured while loading digits"<LoadgAlice()) { - cerr<<"Error occured while LoadgAlice"<GetAliRun(); - if (!gAlice) { - cerr<<"Can't get gAlice !\n"; - return 1; - } - - TDirectory *cwd = gDirectory; - - AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC"); - Int_t ver = tpc->IsVersion(); - cerr<<"TPC version "<CdGAFile(); - - AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60"); - if (!param) {cerr<<"TPC parameters have not been found !\n"; return 4;} - - AliTPCDigitsArray *digarr=new AliTPCDigitsArray; - digarr->Setup(param); - - cerr<<"It has begun"<cd(); - - Int_t nevmax=rl->GetNumberOfEvents();//number of events in run - for(Int_t nev=0;nevGetEvent(nev); - - TTree* input_tree= tpcl->TreeD();//tree with digits - if (input_tree == 0x0){ - cerr << "Can not get TreeD for event " <ConnectTree(input_tree); - - TTree *output_tree =tpcl->TreeR(); - if(output_tree==0x0){ - tpcl->MakeTree("R"); - output_tree = tpcl->TreeR(); - if (output_tree == 0x0){ - cerr << "Problems with output tree (TreeR) for event "<GetEntries(); - - AliTPCclustererKr *clusters = new AliTPCclustererKr(); - clusters->SetParam(param); - clusters->SetInput(input_tree); - clusters->SetOutput(output_tree); - clusters->SetDigArr(digarr); - clusters->FinderIO(); - - tpcl->WriteRecPoints("OVERWRITE"); - } - - - timer.Stop(); timer.Print(); - - delete rl;//cleans everything - - return 0; -} + +/**************************************************************************** + * Origin: A.Matyja amatyja@cern.ch * + ****************************************************************************/ + +/* + + macro to create array of clusters from TPC digits + input files - galice.root + digits.root - file with digits - usualy use link to galice.root + - in splitted mode - neccesary to create link to proper file + + output file - TPC.RecPoints.root + + +// Warning - if cluster file AliTPCclusters.root already exist - macro exit and don't produce anything + + +*/ + + +#ifndef __CINT__ +#include +#include "AliRun.h" +#include "AliTPCv4.h" +#include "AliTPCParam.h" +#include "AliTPCclusterKr.h" +#include "AliTPCclustererKr.h" +#include "TFile.h" +#include "TStopwatch.h" +#include "TTree.h" +#endif + +Int_t FindKrClusters(){ + + // + //Load DataBase + // + //char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB"; + //char *ocdbpath ="local:///home/matyja/baza/OCDB"; + char *ocdbpath ="local:///data/baza/OCDB"; + if (ocdbpath==0){ + ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/"; + } + printf("OCDB PATH = %s\n",ocdbpath); + AliCDBManager * man = AliCDBManager::Instance(); + man->SetDefaultStorage(ocdbpath); + man->SetRun(0); + + AliRunLoader* rl = AliRunLoader::Open("galice.root"); + if (rl == 0x0) { + cerr<<"Can not open session"<GetLoader("TPCLoader"); + if (tpcl == 0x0) { + cerr<<"Can not get TPC Loader"<LoadDigits()) { + cerr<<"Error occured while loading digits"<LoadgAlice()) { + cerr<<"Error occured while LoadgAlice"<GetAliRun(); + if (!gAlice) { + cerr<<"Can't get gAlice !\n"; + return 1; + } + + TDirectory *cwd = gDirectory; + + AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC"); + Int_t ver = tpc->IsVersion(); + cerr<<"TPC version "<CdGAFile(); + + AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60"); + if (!param) {cerr<<"TPC parameters have not been found !\n"; return 4;} + + AliTPCDigitsArray *digarr=new AliTPCDigitsArray; + digarr->Setup(param); + + cerr<<"It has begun"<cd(); + + TTree *output_tree; + + AliTPCclustererKr *clusters = new AliTPCclustererKr(); + clusters->SetParam(param); + clusters->SetOutput(output_tree); + + clusters->SetMinAdc(3);//signal threshold (everything below is treated as 0) + clusters->SetMinTimeBins(2);//number of neighbouring timebins + clusters->SetMaxPadRangeCm(5.);//distance of the cluster center to the center of a pad (in cm) + clusters->SetMaxRowRangeCm(5.);//distance of the cluster center to the center of a padrow (in cm) + clusters->SetMaxTimeRange(7.);//distance of the cluster center to the max time bin on a pad (in tackts) + //ie. fabs(centerT - time)<7 + + clusters->SetIsolCut(3);//set isolation cut threshold + clusters->SetValueToSize(3.1);//cut reduce peak at 0 + + Int_t nevmax=rl->GetNumberOfEvents();//number of events in run + for(Int_t nev=0;nevGetEvent(nev); + + TTree* input_tree= tpcl->TreeD();//tree with digits + if (input_tree == 0x0){ + cerr << "Can not get TreeD for event " <ConnectTree(input_tree); + clusters->SetInput(input_tree); + clusters->SetDigArr(digarr); + cout<<"Processing event "<FinderIO(); + + } + delete clusters; + timer.Stop(); timer.Print(); + + delete rl;//cleans everything + + return 0; +} diff --git a/TPC/FindKrClustersRaw.C b/TPC/FindKrClustersRaw.C index 0569bd9b703..0a00aeadc24 100644 --- a/TPC/FindKrClustersRaw.C +++ b/TPC/FindKrClustersRaw.C @@ -1,165 +1,165 @@ -// - -Int_t FindKrClusterCheck(const char *fileName="data.root"); - - -Int_t FindKrClustersRaw(const char *fileName="data.root"){ - - - - // - // remove Altro warnings - // - AliLog::SetClassDebugLevel("AliTPCRawStream",-5); - AliLog::SetClassDebugLevel("AliRawReaderDate",-5); - AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5); - AliLog::SetModuleDebugLevel("RAW",-5); - // - // Get calibration - // - char *ocdbpath = gSystem->Getenv("OCDB_PATH"); - //char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB"; - if (ocdbpath==0){ - ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/"; - } - printf("OCDB PATH = %s\n",ocdbpath); - AliCDBManager * man = AliCDBManager::Instance(); - man->SetDefaultStorage(ocdbpath); - man->SetRun(100000000); - - AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); - AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); - // - // - - - //define tree - TFile *hfile=new TFile("adc.root","RECREATE","ADC file"); - // Create a ROOT Tree - TTree *mytree = new TTree("Kr","Krypton cluster tree"); - - - Int_t debugLevel=1; - if(debugLevel>0){ - TH1F *histoRow =new TH1F("histoRow","rows",100,0.,100.); - TH1F *histoPad =new TH1F("histoPad","pads",150,0.,150.); - TH1F *histoTime =new TH1F("histoTime","timebins",100,0.,1000.); - TH2F *histoRowPad=new TH2F("histoRowPad","pads-vs-rows" ,150,0.,150.,100,0.,100.); - } - - - //one general output - AliTPCclustererKr *clusters = new AliTPCclustererKr(); - clusters->SetOutput(mytree); - clusters->SetRecoParam(0); - - if(debugLevel>0){ - clusters->SetDebugLevel(debugLevel); - clusters->SetHistoRow(histoRow ); - clusters->SetHistoPad(histoPad ); - clusters->SetHistoTime(histoTime ); - clusters->SetHistoRowPad(histoRowPad); - } - - - AliTPCParamSR *param=new AliTPCParamSR(); - //only for geometry parameters loading - temporarly -// AliRunLoader* rl = AliRunLoader::Open("galice.root"); -// AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60"); - //if (!param) {cerr<<"TPC parameters have not been found !\n"; return 4;} - - clusters->SetParam(param); - - //set cluster finder parameters (from data) - clusters->SetZeroSup(param->GetZeroSup());//zero suppression parameter - clusters->SetFirstBin(60);//first bin - clusters->SetLastBin(950);//last bin - clusters->SetMaxNoiseAbs(2);//maximal noise - clusters->SetMaxNoiseSigma(3);//maximal amount of sigma of noise - - //set cluster finder parameters (from MC) - clusters->SetMinAdc(3);//signal threshold (everything below is treated as 0) - clusters->SetMinTimeBins(2);//number of neighbouring timebins - clusters->SetMaxPadRangeCm(2.5);//distance of the cluster center to the center of a pad (in cm) - clusters->SetMaxRowRangeCm(3.5);//distance of the cluster center to the center of a padrow (in cm) - clusters->SetMaxTimeRange(7);//distance of the cluster center to the max time bin on a pad (in tackts) - //ie. fabs(centerT - time)<7 - clusters->SetValueToSize(6);//cut reduce peak at 0 - - - AliRawReader *reader = new AliRawReaderRoot(fileName); - reader->Reset(); - - TStopwatch timer; - timer.Start(); - - AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader); - stream->SelectRawData("TPC"); - - Int_t evtnr=0; - while (reader->NextEvent()) { - //output for each event - - //if(evtnr>4) break; - cout<<"Evt = "<FinderIO(reader); - evtnr++; - AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr); - } - - - mytree->Print();//print rootuple summary - // Save all objects in this file - hfile->Write(); - // Close the file - hfile->Close(); - - timer.Stop(); - timer.Print(); - printf("Deleting clusterer\n"); - delete clusters; - printf("Deleting stream\n"); - delete stream; - printf("Deleting raw reader\n"); - delete reader; - -// TCanvas *c2=new TCanvas("c2","title",800,800); -// c2->SetHighLightColor(2); -// c2->Range(-458.9552,-2948.238,3296.642,26856.6); -// c2->SetBorderSize(2); -// c2->SetLeftMargin(0.15); -// c2->SetRightMargin(0.06); -// c2->SetFrameFillColor(0); - -// gStyle->SetOptStat(111111); -// histoRow->Draw(); -// c2->Print("rows.ps"); -// histoPad->Draw(); -// c2->Print("pads.ps"); -// histoTime->Draw(); -// c2->Print("timebins.ps"); -// histoRowPad->Draw(); -// c2->Print("row-pad.ps"); - - return 0; -} - - -Int_t FindKrClusterCheck(const char *fileName){ - // - // - gSystem->Load("$ROOTSYS/lib/libGui.so"); - gSystem->Load("$ROOTSYS/lib/libTree.so"); - gSystem->Load("$MEMSTAT/libMemStat.so"); - { - TMemStat memstat(1000000,100000,kTRUE); - AliSysInfo::AddCallBack(TMemStatManager::GetInstance()->fStampCallBack); - FindKrClustersRaw(fileName); - } - // the output memstat.root file - TMemStat draw("memstat.root"); - // Print some information - // code info - draw.MakeReport(0,0); - -} +// + +Int_t FindKrClusterCheck(const char *fileName="data.root"); + + +Int_t FindKrClustersRaw(const char *fileName="data.root"){ + + + + // + // remove Altro warnings + // + AliLog::SetClassDebugLevel("AliTPCRawStream",-5); + AliLog::SetClassDebugLevel("AliRawReaderDate",-5); + AliLog::SetClassDebugLevel("AliTPCAltroMapping",-5); + AliLog::SetModuleDebugLevel("RAW",-5); + // + // Get calibration + // + char *ocdbpath = gSystem->Getenv("OCDB_PATH"); + //char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB"; + if (ocdbpath==0){ + ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/"; + } + printf("OCDB PATH = %s\n",ocdbpath); + AliCDBManager * man = AliCDBManager::Instance(); + man->SetDefaultStorage(ocdbpath); + man->SetRun(100000000); + + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); + AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); + // + // + + + //define tree + TFile *hfile=new TFile("adc.root","RECREATE","ADC file"); + // Create a ROOT Tree + TTree *mytree = new TTree("Kr","Krypton cluster tree"); + + + Int_t debugLevel=1; + if(debugLevel>0){ + TH1F *histoRow =new TH1F("histoRow","rows",100,0.,100.); + TH1F *histoPad =new TH1F("histoPad","pads",150,0.,150.); + TH1F *histoTime =new TH1F("histoTime","timebins",100,0.,1000.); + TH2F *histoRowPad=new TH2F("histoRowPad","pads-vs-rows" ,150,0.,150.,100,0.,100.); + } + + + //one general output + AliTPCclustererKr *clusters = new AliTPCclustererKr(); + clusters->SetOutput(mytree); + clusters->SetRecoParam(0); + + if(debugLevel>0){ + clusters->SetDebugLevel(debugLevel); + clusters->SetHistoRow(histoRow ); + clusters->SetHistoPad(histoPad ); + clusters->SetHistoTime(histoTime ); + clusters->SetHistoRowPad(histoRowPad); + } + + + AliTPCParamSR *param=new AliTPCParamSR(); + //only for geometry parameters loading - temporarly +// AliRunLoader* rl = AliRunLoader::Open("galice.root"); +// AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60"); + //if (!param) {cerr<<"TPC parameters have not been found !\n"; return 4;} + + clusters->SetParam(param); + + //set cluster finder parameters (from data) + clusters->SetZeroSup(param->GetZeroSup());//zero suppression parameter + clusters->SetFirstBin(60);//first bin + clusters->SetLastBin(950);//last bin + clusters->SetMaxNoiseAbs(2);//maximal noise + clusters->SetMaxNoiseSigma(3);//maximal amount of sigma of noise + + //set cluster finder parameters (from MC) + clusters->SetMinAdc(3);//signal threshold (everything below is treated as 0) + clusters->SetMinTimeBins(2);//number of neighbouring timebins + clusters->SetMaxPadRangeCm(5.);//distance of the cluster center to the center of a pad (in cm) + clusters->SetMaxRowRangeCm(5.);//distance of the cluster center to the center of a padrow (in cm) + clusters->SetMaxTimeRange(5.);//distance of the cluster center to the max time bin on a pad (in tackts) + //ie. fabs(centerT - time)<7 + clusters->SetValueToSize(7.);//cut reduce peak at 0 + clusters->SetIsolCut(3);//set isolation cut threshold + + AliRawReader *reader = new AliRawReaderRoot(fileName); + reader->Reset(); + + TStopwatch timer; + timer.Start(); + + AliAltroRawStreamFast* stream = new AliAltroRawStreamFast(reader); + stream->SelectRawData("TPC"); + + Int_t evtnr=0; + while (reader->NextEvent()) { + //output for each event + + //if(evtnr>4) break; + cout<<"Evt = "<FinderIO(reader); + evtnr++; + AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr); + } + + + mytree->Print();//print rootuple summary + // Save all objects in this file + hfile->Write(); + // Close the file + hfile->Close(); + + timer.Stop(); + timer.Print(); + printf("Deleting clusterer\n"); + delete clusters; + printf("Deleting stream\n"); + delete stream; + printf("Deleting raw reader\n"); + delete reader; + +// TCanvas *c2=new TCanvas("c2","title",800,800); +// c2->SetHighLightColor(2); +// c2->Range(-458.9552,-2948.238,3296.642,26856.6); +// c2->SetBorderSize(2); +// c2->SetLeftMargin(0.15); +// c2->SetRightMargin(0.06); +// c2->SetFrameFillColor(0); + +// gStyle->SetOptStat(111111); +// histoRow->Draw(); +// c2->Print("rows.ps"); +// histoPad->Draw(); +// c2->Print("pads.ps"); +// histoTime->Draw(); +// c2->Print("timebins.ps"); +// histoRowPad->Draw(); +// c2->Print("row-pad.ps"); + + return 0; +} + + +Int_t FindKrClusterCheck(const char *fileName){ + // + // + gSystem->Load("$ROOTSYS/lib/libGui.so"); + gSystem->Load("$ROOTSYS/lib/libTree.so"); + gSystem->Load("$MEMSTAT/libMemStat.so"); + { + TMemStat memstat(1000000,100000,kTRUE); + AliSysInfo::AddCallBack(TMemStatManager::GetInstance()->fStampCallBack); + FindKrClustersRaw(fileName); + } + // the output memstat.root file + TMemStat draw("memstat.root"); + // Print some information + // code info + draw.MakeReport(0,0); + +} -- 2.43.0