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