--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2007, 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. *
+ **************************************************************************/
+
+/* $Log$ */
+
+//_________________________________________________________________________
+// Calibration coefficients
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+#include "AliEMCALCalibCoefs.h"
+
+#include "AliRun.h"
+#include "AliEMCALCalibData.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBId.h"
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliEMCALGeometry.h"
+
+#include "TH1F.h"
+
+TableClassImpl(AliEMCALCalibCoefs,calibCoef)
+
+// Get initial Calib Data from DB
+AliEMCALCalibCoefs* AliEMCALCalibCoefs::GetCalibTableFromDb(const char *tn)
+{
+ // Initial cc with decalibration
+ char* dbString = "local:///data/r22b/ALICE/PROD/CALIBRATION_May_2007/PI0/PDSF/10GEV/DECALIB/DeCalibDB";
+ AliEMCALCalibData* caldata = (AliEMCALCalibData*)
+ (AliCDBManager::Instance()->GetStorage(dbString)
+ ->Get("EMCAL/Calib/Data",gAlice->GetRunNumber())->GetObject());
+ if(caldata == 0) return 0;
+
+ AliEMCALCalibCoefs *tab = new AliEMCALCalibCoefs(tn);
+ tab->SetCalibMethod(AliEMCALCalibCoefs::kMC);
+
+ AliEMCALGeometry *g = AliEMCALGeometry::GetInstance();
+ for(Int_t id=0; id<g->GetNCells(); id++){
+ Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0, iphiCell=0, ietaCell=0;
+ g->GetCellIndex(id, nSupMod, nModule, nIphi, nIeta);
+ g->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta, iphiCell, ietaCell);
+
+ calibCoef r;
+ r.absId = id;
+ r.cc = caldata->GetADCchannel(nSupMod, ietaCell, iphiCell);
+ r.eCc = 0.0;
+
+ tab->AddAt(&r);
+ }
+ tab->Purge();
+ tab->SetCalibMethod(AliEMCALCalibCoefs::kPI0);
+
+ return tab;
+}
+
+TH1F* AliEMCALCalibCoefs::GetHistOfCalibTableFromDb(const char *tn)
+{
+ // First SM only
+ AliEMCALGeometry *g = AliEMCALGeometry::GetInstance("");
+
+ AliEMCALCalibCoefs* tab = GetCalibTableFromDb(tn);
+ if(tab==0) return 0;
+
+ TH1F *h = new TH1F("hCCfirst", " cc first (in MeV)", 70, 12., 19.);
+ calibCoef *r;
+ for(Int_t i=0; i<tab->GetNRows(); i++){
+ r = tab->GetTable(i);
+ if(i>=1152) break;
+ h->Fill(r->cc*1000.);
+ }
+ delete tab;
+ return h;
+}
+
+calibCoef *AliEMCALCalibCoefs::GetRow(const int absId)
+{
+ calibCoef *r=0;
+ for(int id=0; id<GetNRows(); id++){
+ r = GetTable(id);
+ if(r->absId == absId) return r;
+ }
+ return 0;
+}
+
+void AliEMCALCalibCoefs::PrintTable()
+{
+ printf(" Table : %s : nrows %i \n", GetName(), int(GetNRows()));
+ for(int i=0; i<GetNRows(); i++) PrintTable(i);
+}
+
+void AliEMCALCalibCoefs::PrintTable(const Int_t i)
+{
+ if(i>=GetNRows()) return;
+ printf("row %i \n", i);
+ PrintRec(GetTable(i));
+}
+
+void AliEMCALCalibCoefs::PrintRec(calibCoef* r)
+{
+ if(r==0) return;
+ printf(" abs Id %5.5i cc %7.6f eCc %7.6f \n", r->absId, r->cc, r->eCc);
+}
--- /dev/null
+#ifndef ALIEMCALCALIBCOEFS_H
+#define ALIEMCALCALIBCOEFS_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+// Table of Calibration coefficients
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+// --- ROOT system ---
+#include <TTable.h>
+
+// unit is GeV
+struct calibCoef {
+ Int_t absId; // absolute id of cell
+ Double_t cc; // Calib. coef
+ Double_t eCc; // Calib. coef. error
+};
+
+class TH1F;
+
+class AliEMCALCalibCoefs : public TTable {
+ public:
+ enum EEmcalCalibType {kMC, kEQUALIZATION, kMIP, kPI0}; // type of EMCAL calibrations
+
+ void SetCalibMethod(Int_t var) {fCalibMethod=var;}
+ Int_t GetCalibMethod() {return fCalibMethod;}
+ calibCoef* GetRow(const int absId);
+ // Get initial Calib Data from DB
+ static AliEMCALCalibCoefs *GetCalibTableFromDb(const char *tn="CCIN");
+ static TH1F *GetHistOfCalibTableFromDb(const char *tn="CCTMP");
+ // Menu
+ void PrintTable(); // *MENU*
+ void PrintTable(const Int_t i); // *MENU*
+ void PrintRec(calibCoef *r);
+
+ protected:
+ Int_t fCalibMethod; // method of calibration - EEmcalCalibType
+
+ ClassDefTable(AliEMCALCalibCoefs , calibCoef)
+ ClassDef(AliEMCALCalibCoefs,1) // Table of Calibration coefficients
+};
+
+#endif // ALIEMCALCalibCoefs_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2007, 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:*/
+
+//_________________________________________________________________________
+// Top EMCAL folder which will keep all information about EMCAL itself,
+// super Modules (SM), modules, towers, set of hists and so on.
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+#include "AliEMCALCell.h"
+#include "AliEMCALHistoUtilities.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALFolder.h"
+#include "AliEMCALSuperModule.h"
+#include "AliEMCALCalibData.h"
+
+#include "AliEMCALCalibCoefs.h"
+
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TList.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TNtuple.h>
+#include <TObjectSet.h>
+
+typedef AliEMCALHistoUtilities u;
+
+ClassImp(AliEMCALCell)
+
+Double_t ADCCHANNELEC = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
+Double_t MPI0 = 0.13498; // mass of pi0
+Double_t MPI02 = MPI0*MPI0; // mass**2
+
+AliEMCALCell::AliEMCALCell() :
+TObjectSet(),
+fAbsId(0),fSupMod(0),fModule(0),fPhi(0),fEta(0),fPhiCell(0),fEtaCell(0),fCcIn(0),fCcOut(0),
+fFun(0)
+{
+}
+
+AliEMCALCell::AliEMCALCell(const Int_t absId, const char* title) :
+TObjectSet(Form("Cell%4.4i",absId)),
+fAbsId(absId),fSupMod(0),fModule(0),fPhi(0),fEta(0),fPhiCell(0),fEtaCell(0),fCcIn(0),fCcOut(0),
+fFun(0)
+{
+ SetTitle(title);
+
+ AliEMCALGeometry *g = AliEMCALGeometry::GetInstance();
+ g->GetCellIndex(fAbsId, fSupMod, fModule, fPhi, fEta);
+ g->GetCellPhiEtaIndexInSModule(fSupMod, fModule, fPhi, fEta, fPhiCell, fEtaCell);
+
+}
+
+AliEMCALCell::~AliEMCALCell()
+{
+ // dtor
+}
+
+void AliEMCALCell::SetCCfromDB(AliEMCALCalibData *ccDb)
+{
+ if(ccDb == 0) return;
+ // fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
+ // fetaCel-column; fPhiCell- row
+ fCcIn = ccDb->GetADCchannel(fSupMod, fEtaCell, fPhiCell); // in GeV
+
+ TH1* h = (TH1*)GetHists()->At(0);
+ u::AddToNameAndTitle(h, 0, Form(", cc %5.2f MeV", fCcIn*1.e+3));
+}
+
+void AliEMCALCell::SetCCfromCCTable(AliEMCALCalibCoefs *t)
+{
+ if(t == 0) return;
+ this->AddObject((TObject*)BookHists(), kTRUE);
+
+ calibCoef *r = t->GetTable(fAbsId);
+ if(r && r->absId == fAbsId) {
+ fCcIn = r->cc;
+ } else { // something wrong
+ if(r) printf(" fAbsId %i : r->absId %i \n", fAbsId, r->absId);
+ assert(0);
+ }
+
+ TH1* h = (TH1*)GetHists()->At(0);
+ u::AddToNameAndTitle(h, 0, Form(", cc %5.2f MeV", fCcIn*1.e+3));
+}
+
+void AliEMCALCell::FillEffMass(const Double_t mgg)
+{
+ u::FillH1(GetHists(), 0, mgg);
+}
+
+void AliEMCALCell::FillCellNtuple(TNtuple *nt)
+{
+ if(nt==0) return;
+ nt->Fill(fAbsId,fSupMod,fModule,fPhi,fEta,fPhiCell,fEtaCell,fCcIn,fCcOut);
+}
+
+void AliEMCALCell::FitHist(TH1* h, const char* name, const char* opt)
+{
+ TString optFit(""), OPT(opt);
+ OPT.ToUpper();
+ if(h==0) return;
+ printf("<I> AliEMCALCell::FitHist : h %p |%s| is started : opt %s\n", h, h->GetName(), opt);
+ TString tit(h->GetTitle());
+
+ TF1 *GausPol2 = 0, *g=0, *bg=0;
+ if(h->GetListOfFunctions()->GetSize() == 0) {
+ g = u::Gausi(name, 0.0, 0.4, h); // gaus estimation
+
+ g->SetParLimits(0, h->Integral()/20., h->Integral());
+ g->SetParLimits(1, 0.08, 0.20);
+ g->SetParLimits(2, 0.001, 0.02);
+
+ g->SetParameter(0, 1200.);
+ g->SetParameter(1, h->GetBinLowEdge(h->GetMaximumBin()));
+ g->SetParameter(2, 0.010);
+
+ g->SetLineColor(kRed);
+ g->SetLineWidth(2);
+
+ optFit = "0NQ";
+ h->Fit(g, optFit.Data(),"", 0.001, 0.3);
+
+ optFit = "0NQIME";
+ h->Fit(g, optFit.Data(),"", 0.001, 0.3);
+
+ bg = new TF1(Form("bg%s",name), "pol2", 0.0, 0.3);
+ optFit = "0NQ";
+ h->Fit(bg, optFit.Data(),"", 0.0, 0.3);
+
+ GausPol2 = u::GausiPol2(name, 0.00, 0.3, g, bg);
+ optFit = "0Q";
+ h->Fit(GausPol2, optFit.Data(),"", 0.03, 0.28);
+ // Clean up
+ delete g;
+ delete bg;
+ optFit = "0IME+"; // no drwaing at all
+ if(tit.Contains("SM") || OPT.Contains("DRAW")) optFit = "IME+";
+ } else {
+ GausPol2 = (TF1*)h->GetListOfFunctions()->At(0);
+ optFit = "IME+";
+ }
+ // optFit = "IME+";
+ h->Fit(GausPol2, optFit.Data(),"", 0.01, 0.28);
+
+ if(optFit.Contains("0") == 0) {
+ gStyle->SetOptFit(111);
+ u::DrawHist(h,2);
+ }
+ printf("<I> AliEMCALCell::FitHist : h %p |%s| is ended \n\n", h, h->GetName());
+}
+
+void AliEMCALCell::FitEffMassHist(const char* opt)
+{
+ AliEMCALFolder* EMCAL = (AliEMCALFolder*)(GetParent()->GetParent()->GetParent());
+ Int_t it = EMCAL->GetIterationNumber();
+
+ TH1* h = (TH1*)GetHists()->At(0);
+
+ FitHist(h, GetName(), opt);
+
+ fFun = (TF1*)h->GetListOfFunctions()->At(0);
+ if(fFun) {
+ Double_t mpi = fFun->GetParameter(1), mpi2 = mpi*mpi;
+ Double_t ccTmp = fCcIn * MPI02 / mpi2;
+ if(it<=6) {
+ fCcOut = ccTmp;
+ } else {
+ fCcOut = (ccTmp + fCcIn)/2.;
+ }
+ printf(" fFun %s | %s : iet %i\n", fFun->GetName(), fFun->GetTitle(), it);
+ }
+ printf(" %s | fCcIn %6.5f -> % 6.5f <- fCcOut \n", GetTitle(), fCcIn , fCcOut);
+}
+
+void AliEMCALCell::Print()
+{
+ TH1* h = (TH1*)GetHists()->At(0);
+ TF1 *f = (TF1*)h->GetListOfFunctions()->At(0);
+ printf(" %s %s \n", GetName(), GetTitle());
+ if(fFun) printf(" fFun : %s | %s \n", fFun->GetName(), fFun->GetTitle());
+ else fFun = f;
+ if(f) f->Dump();
+}
+
+TList* AliEMCALCell::BookHists()
+{
+ gROOT->cd();
+ TH1::AddDirectory(1);
+
+ AliEMCALFolder* EMCAL = (AliEMCALFolder*)(GetParent()->GetParent()->GetParent());
+ Int_t it = EMCAL->GetIterationNumber();
+
+ new TH1F("01_EffMass", "effective mass of #gamma,#gamma(m_{#pi^{0}}=134.98 MeV) ", 60,0.0,0.3);
+
+ TList *l = AliEMCALHistoUtilities::MoveHistsToList(Form("HistsOfCell%4.4i",fAbsId), kFALSE);
+ AliEMCALHistoUtilities::AddToNameAndTitleToList(l, Form("4.4i_It%i",fAbsId,it),
+ Form(" Cell %4.4i, iter. %i",fAbsId, it));
+
+ TH1::AddDirectory(0);
+ return l;
+}
--- /dev/null
+#ifndef ALIEMCALCELL_H
+#define ALIEMCALCELL_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+// EMCAL cell - keep everyrhing for calibration task
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+#include <TObjectSet.h>
+
+class TList;
+class TH1;
+class TF1;
+class TNtuple;
+
+class AliEMCALCalibData;
+class AliEMCALCalibCoefs;
+
+class AliEMCALCell : public TObjectSet {
+
+ public:
+
+ AliEMCALCell();
+ AliEMCALCell(const Int_t absId, const char* title="EMCAL cell");
+
+ virtual ~AliEMCALCell();
+
+ void SetCCfromDB(AliEMCALCalibData *ccDb); // obsolete
+ void SetCCfromCCTable(AliEMCALCalibCoefs *t);
+
+ TList* GetHists() {return (TList*)fObj;}
+ Int_t GetAbsId() const {return fAbsId;}
+ Int_t GetSupMod() const {return fSupMod;}
+ Int_t GetModule() const {return fModule;}
+
+ Double_t GetCcIn() {return fCcIn;}
+ Double_t GetCcOut() {return fCcOut;}
+ TF1* GetFunction() {return fFun;}
+
+ void FillEffMass(const Double_t mgg);
+ void FillCellNtuple(TNtuple *nt);
+
+ static void FitHist(TH1* h, const char* name="",const char* opt="");
+ // Menu
+ void FitEffMassHist(const char* opt=""); //*MENU*
+ void Print(); //*MENU*
+ protected:
+ Int_t fAbsId; // abs cell id
+ Int_t fSupMod; // super module number
+ Int_t fModule; // module number inside SM
+ Int_t fPhi; // phi number of cell inside module
+ Int_t fEta; // eta number of cell inside module
+ Int_t fPhiCell; // phi number of cell SM
+ Int_t fEtaCell; // eta number of cell SM
+ // CC staf
+ Double_t fCcIn; // input cc in GeV (from Db or table
+ Double_t fCcOut; // output cc in GeV (from fir now)
+
+ TF1* fFun; //! fitting function - gaus + pol2
+ //
+ TList* BookHists();
+
+ ClassDef(AliEMCALCell,1) // EMCAL cell
+
+};
+
+#endif // ALIEMCALCELL_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2007, 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$ */
+
+//_________________________________________________________________________
+// Top EMCAL folder which will keep all information about EMCAL itself,
+// super Modules (SM), modules, towers, set of hists and so on.
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+#include "AliEMCALFolder.h"
+#include "AliEMCALHistoUtilities.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALSuperModule.h"
+#include "AliEMCALCell.h"
+#include "AliESDCaloCluster.h"
+
+#include "AliRun.h"
+
+#include "AliEMCALCalibData.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBId.h"
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+
+#include "AliEMCALCalibCoefs.h"
+#include "AliEMCALDigit.h"
+#include "AliEMCALRecPoint.h"
+
+#include "AliEMCALPi0SelectionParam.h"
+
+#include <cassert>
+
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TList.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TFile.h>
+#include <TFileIter.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TKey.h>
+#include <TNtuple.h>
+#include <TLegend.h>
+#include <TArrayS.h>
+
+ const TString AliEMCALFolder::fgkBaseFolderName("EMCAL");
+ const TString AliEMCALFolder::fgkCCFirstName("CCFIRST");
+ const TString AliEMCALFolder::fgkCCinName("CCIN");
+ const TString AliEMCALFolder::fgkCCoutName("CCOUT");
+
+typedef AliEMCALHistoUtilities u;
+
+ClassImp(AliEMCALFolder)
+
+//AliEMCALGeometry* AliEMCALFolder::fGeometry = 0;
+
+ TList *lObj = 0;
+
+AliEMCALFolder::AliEMCALFolder() :
+ TObjectSet(),
+ fCounter(0), fGeometry(0), fNumOfCell(0), fLhists(0), fLofCells(0),fPi0SelPar(0),fCalibData(0),
+ fCellNtuple(0)
+{
+}
+
+AliEMCALFolder::AliEMCALFolder(const char* name, const char* title, Bool_t putToBrowser) :
+ TObjectSet(name),
+ fCounter(-1), fGeometry(0), fNumOfCell(0), fLhists(0), fLofCells(0),fPi0SelPar(0),fCalibData(0),
+ fCellNtuple(0)
+{
+ SetTitle(title);
+ Init(putToBrowser);
+}
+
+AliEMCALFolder::AliEMCALFolder(const Int_t it, const char* title, Bool_t putToBrowser) :
+ TObjectSet(Form("%s_%2.2i", AliEMCALFolder::fgkBaseFolderName.Data(),it)),
+ fCounter(it), fGeometry(0), fNumOfCell(0), fLhists(0), fLofCells(0),fPi0SelPar(0),fCalibData(0),
+ fCellNtuple(0)
+{
+ SetTitle(title);
+ Init(putToBrowser);
+}
+
+AliEMCALFolder::~AliEMCALFolder()
+{
+ // dtor
+}
+
+void AliEMCALFolder::Init(Bool_t putToBrowser)
+{
+ lObj = new TList;
+ lObj->SetName("Objects"); // name is good ?
+ this->AddObject((TObject*)lObj, kTRUE);
+ // Get default geometry - "SHISH_77_TRD1_2X2_FINAL_110DEG"; May 29, 2007
+ fGeometry = AliEMCALGeometry::GetInstance(); // should be define before
+ lObj->Add(fGeometry);
+
+ // Initial cc with decalibration
+ char* dbString = "local:///data/r22b/ALICE/PROD/CALIBRATION_May_2007/PI0/PDSF/10GEV/DECALIB/DeCalibDB"; // hard coded now
+ fCalibData = (AliEMCALCalibData*)
+ (AliCDBManager::Instance()->GetStorage(dbString)
+ ->Get("EMCAL/Calib/Data",gAlice->GetRunNumber())->GetObject());
+ lObj->Add(fCalibData);
+ // Jun 13, 2007
+ Add(AliEMCALCalibCoefs::GetCalibTableFromDb(fgkCCFirstName.Data()));
+ if(GetIterationNumber()==1) {
+ Add(AliEMCALCalibCoefs::GetCalibTableFromDb()); // first iteration only
+ }
+ // Selection Parameter
+ fPi0SelPar = AliEMCALPi0SelectionParam::Set1();
+ this->Add(fPi0SelPar);
+ //
+ fLhists = BookHists();
+ lObj->Add(fLhists);
+
+ // dimension should be get from geometry - 4*12*24*11);
+ fNumOfCell = fGeometry->GetNCells();
+ fLofCells = new AliEMCALCell*[fNumOfCell];
+ for(int i=0; i<fNumOfCell; i++) fLofCells[i] = 0;
+
+ printf("<I> Create AliEMCALFolder : it %i : name %s\n ", fCounter, GetName());
+
+ if(putToBrowser) gROOT->GetListOfBrowsables()->Add(this); // for testing purpuse
+}
+
+AliEMCALSuperModule* AliEMCALFolder::GetSuperModule(const Int_t nm)
+{
+ TDataSet *set = 0;
+ AliEMCALSuperModule *sm = 0;
+
+ set = FindByName(Form("SM%2.2i",nm));
+ if(set) sm = (AliEMCALSuperModule*)set;
+
+ return sm;
+}
+
+
+AliEMCALCell* AliEMCALFolder::GetCell(const Int_t absId)
+{ // May 30, 2007
+ if(absId<0 || absId >= fNumOfCell) return 0;
+ else return fLofCells[absId];
+}
+
+void AliEMCALFolder::SetCell(AliEMCALCell *cell, const Int_t absId)
+{
+ if(absId>=0 && absId < fNumOfCell) {
+ fLofCells[absId] = cell;
+ }
+}
+
+pi0SelectionParam* AliEMCALFolder::GetPi0SelectionParRow(Int_t nrow)
+{
+ pi0SelectionParam* r=0;
+ if(fPi0SelPar) {
+ r = fPi0SelPar->GetTable(nrow);
+ }
+ return r;
+}
+
+void AliEMCALFolder::FillPi0Candidate(const Double_t mgg, AliESDCaloCluster* cl1, AliESDCaloCluster* cl2)
+{
+ static Int_t absIdMax, nm1, nm2;
+
+ u::FillH1(fLhists, 0, 1.); // number entries
+ u::FillH1(fLhists, 1, mgg);
+
+ nm1 = GetSMNumber(cl1);
+ nm2 = GetSMNumber(cl2);
+
+ if(nm1==-1 || nm2==-1) assert(0);
+
+ if(nm1 != nm2) return; // Both cluster should be in the same SM
+
+ AliESDCaloCluster* cl = cl1;
+ if(cl1->E() < cl2->E()) cl = cl2; // Get cluster with highest energy
+
+ const Int_t nDigits = cl->GetNumberOfDigits();
+ // const UShort_t* adc = cl->GetDigitAmplitude(); // for future
+ const Short_t* absId = cl->GetDigitIndex()->GetArray();
+
+ int indMax = 0;
+ double emax = cl->GetDigitAmplitude()->At(0);
+ if(nDigits > 1) {
+ for(int i=1; i<nDigits; i++) {
+ if(emax < cl->GetDigitAmplitude()->At(i)) {
+ indMax = i;
+ emax = cl->GetDigitAmplitude()->At(i);
+ }
+ }
+ }
+ if(emax/cl->E() > 0.5) { // more than 50% of cluster energy
+ u::FillH1(fLhists, 0, 2.); // number "good" entries
+ absIdMax = Int_t(absId[indMax]);
+ FillPi0Candidate(mgg, absIdMax, nm1);
+ }
+}
+
+void AliEMCALFolder::FillPi0Candidate(const Double_t mgg, Int_t absIdMax, Int_t nm)
+{
+ //
+ // Jun 08;
+ //
+ static Int_t nSupModMax, nModuleMax, nIphiMax, nIetaMax, iphiCellMax, ietaCellMax;
+ static AliEMCALCell* cell;
+ static TDataSet *set;
+ static AliEMCALSuperModule *sm;
+
+ fGeometry->GetCellIndex(absIdMax, nSupModMax, nModuleMax, nIphiMax, nIetaMax);
+ if(nm != nSupModMax) assert(0);
+
+ fGeometry->GetCellPhiEtaIndexInSModule(nSupModMax, nModuleMax, nIphiMax, nIetaMax, iphiCellMax, ietaCellMax);
+
+ cell = 0;
+ set = 0;
+ sm = 0;
+ if(GetCell(absIdMax)==0) {
+ cell = new AliEMCALCell(absIdMax, Form("sm%2.2i:phi%2.2i:eta%2.2i(%4.4i)",nSupModMax,iphiCellMax,ietaCellMax,absIdMax));
+ SetCell(cell, absIdMax);
+ // For browser
+ set = FindByName(Form("SM%2.2i",nSupModMax));
+ if(set==0) {
+ sm = new AliEMCALSuperModule(nSupModMax);
+ Add(sm);
+ sm->Init();
+ } else {
+ sm = (AliEMCALSuperModule*)set;
+ }
+ sm->AddCellToEtaRow(cell, ietaCellMax);
+ //
+ cell->SetCCfromCCTable(GetCCIn());
+ } else {
+ cell = GetCell(absIdMax);
+ set = FindByName(Form("SM%2.2i",nm));
+ if(set) sm = (AliEMCALSuperModule*)set;
+ }
+ if(sm == 0) assert(0);
+ if(nm != sm->GetSMNumber()) assert(0);
+
+ u::FillH1(sm->GetHists(), 0, mgg);
+ cell->FillEffMass(mgg);
+}
+
+void AliEMCALFolder::FitAllSMs()
+{ // Jun 14, 2007
+ // Only first SM now - should be changed in the future
+ AliEMCALSuperModule *sm0 = GetSuperModule(0);
+ sm0->FitForAllCells();
+ // Get input calibration table
+ AliEMCALCalibCoefs *ccIn = GetCCIn();
+ if(ccIn == 0) {
+ printf("<E> no input cc \n");
+ return;
+ }
+ // New calibration table
+ AliEMCALCalibCoefs *ccOut = new AliEMCALCalibCoefs(fgkCCoutName.Data(), ccIn->GetNRows());
+ calibCoef *rIn=0, rOut;
+ for(Int_t i=0; i<ccIn->GetNRows(); i++){
+ rIn = ccIn->GetTable(i);
+ rOut = (*rIn);
+ AliEMCALCell* cell = GetCell(rIn->absId);
+ if(cell && cell->GetSupMod() == 0) { // only first module now
+ rOut.cc = cell->GetCcOut();
+ }
+ ccOut->AddAt(&rOut);
+ }
+ ccOut->Purge();
+ Add(ccOut);
+}
+
+AliEMCALCalibCoefs* AliEMCALFolder::GetCCTable(const char* name)
+{
+ TDataSet *ds = FindByName(name);
+ if(ds) return (AliEMCALCalibCoefs*)ds;
+ else return 0;
+}
+
+Int_t AliEMCALFolder::GetSMNumber(AliESDCaloCluster* cl)
+{
+ static Int_t absId, nSupMod, nModule, nIphi, nIeta;
+ nSupMod = -1; // if negative something wrong
+ if(cl) {
+ absId = Int_t(cl->GetDigitIndex()->At(0));
+ fGeometry->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
+ }
+ return nSupMod;
+}
+
+// Recalibration staf - Jun 18,2007
+AliEMCALRecPoint* AliEMCALFolder::GetRecPoint(AliESDCaloCluster *cl, AliEMCALCalibCoefs *tOld,AliEMCALCalibCoefs *tNew,
+TList *l)
+{
+ //
+ // Static function;
+ // Get recalibrated rec.point from ESD cluster
+ //
+ // printf(" AliEMCALFolder::GetRecPoint() : RECALIBRATION \n");
+ static Double_t ADCCHANNELEC = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
+ static Float_t ECAW0 = 4.5; // hard coded now - see AliEMCALClusterizerv1::InitParameters()
+ Int_t ampDigi=0;
+ Double_t eDigiNew=0.0;
+
+ AliEMCALRecPoint *rp=0;
+ calibCoef *rOld=0, *rNew=0;
+ if(cl && tOld && tNew){
+ // cl->PrintClusterInfo(1);
+ const Int_t ndg = cl->GetNumberOfDigits();
+ //const Int_t nprim = cl->GetNumberOfPrimaries();
+
+ const Short_t* prim = cl->GetLabels()->GetArray();
+ const Short_t* dgAbsId = cl->GetDigitIndex()->GetArray();
+ // const UShort_t* dgAmp = cl->GetDigitAmplitude()(); // This is energy - bad definition
+
+ rp = new AliEMCALRecPoint(""); // opt=""
+ rp->SetClusterType(AliESDCaloCluster::kClusterv1);
+ AliEMCALDigit* dg=0;
+ TClonesArray digits("AliEMCALDigit", ndg);
+ Int_t absId = 0;
+ for(Int_t i=0; i<ndg; i++){
+ // Save just abs id and amplitude of digits which will be used for recalculation of
+ // cluster energy and position
+ absId = Int_t(dgAbsId[i]);
+ rOld = tOld->GetTable(absId);
+ rNew = tNew->GetTable(absId);
+ ampDigi = Int_t(cl->GetDigitAmplitude()->At(i) / rOld->cc);
+
+ new(digits[i]) AliEMCALDigit(Int_t(prim[i]),0,absId, ampDigi, 0.0, i, 0.0);
+ dg = (AliEMCALDigit*)digits[i];
+ eDigiNew = ampDigi * rNew->cc;
+ //eDigiNew = ampDigi * ADCCHANNELEC;
+ //eDigiNew = Double_t(cl->GetDigitAmplitude()->At(i)); // Copy from ESD for checking
+ rp->AddDigit(*dg, Float_t(eDigiNew));
+ u::FillH1(l, 14, eDigiNew);
+ u::FillH1(l, 15, Double_t(absId));
+ //printf("<I> digit %i amp %i rOld->cc %6.5f GeV rNew->cc %6.5f GeV\n", i, ampDigi, rOld->cc, rNew->cc);
+ }
+ //printf("<I> recalibration of digits was done ! \n");
+ // rp->EvalAll(ECAW0, &digits);
+ rp->EvalLocalPosition(ECAW0, &digits); // I need just position
+ digits.Delete();
+ }
+ //rp->Print("print");
+ return rp;
+}
+
+void AliEMCALFolder::Save(const char *fn, const char *opt)
+{
+ //
+ // Jun 5, 2007; See TFileIter and StFMC.cxx
+ //
+ TString FN(fn);
+ if(FN.Contains(".root")==0) FN += ".root";
+ TFileIter f(FN.Data(),opt,"EMCAL object");
+ UInt_t eventNum = 0; // just one object
+ UInt_t runNumber = 0; // 0 now, - may statistics on selectorr
+ f.NextEventPut(this, eventNum, runNumber);
+ printf(" Save %s to file %s\n", GetName(), FN.Data());
+}
+
+AliEMCALFolder* AliEMCALFolder::Read(const char *fn, const char *opt)
+{
+ //
+ // Jun 27, 2007
+ //
+ AliEMCALFolder* EMCAL = 0;
+ TH1::AddDirectory(0); // this is obligatory
+
+ TFile f(fn,opt);
+ TList *l = f.GetListOfKeys();
+ printf(" The total number of the objects: %d in file %s\n", l->GetSize(), fn);
+
+ TKey *key = (TKey*)l->At(0);
+ EMCAL = (AliEMCALFolder*)key->ReadObj();
+ f.Close();
+ if(EMCAL) EMCAL->InitAfterRead();
+
+ return EMCAL;
+}
+
+
+void AliEMCALFolder::InitAfterRead()
+{
+ lObj = (TList*)fObj;
+ fLhists = (TList*)lObj->FindObject("HistsOfEmcal");
+}
+
+void AliEMCALFolder::DrawQA(const int nsm)
+{
+ //
+ // Jun 25, 2007
+ //
+
+ AliEMCALSuperModule* sm = GetSuperModule(nsm);
+ if(sm==0) return;
+ TList *l = sm-> GetHists();
+ Int_t nx=2, ny=2, wh=530, ww=750;
+
+ TCanvas *c = new TCanvas(Form("QA_%i",GetIterationNumber()), Form("QA_%i",GetIterationNumber()),
+ 10, 10, ww, wh);
+ c->Divide(nx,ny);
+
+ int ic=1;
+ c->cd(ic++);
+ TH1 *h1 = (TH1*)l->At(0);
+ sm->FitEffMassHist();
+ u::DrawHist(h1,2);
+ h1->SetAxisRange(0.03, 0.28);
+
+ c->cd(ic++);
+ sm->DrawCC(0);
+ TH1 *hccin = (TH1*)l->At(1);
+ hccin->SetAxisRange(14., 20.);
+
+ gStyle->SetOptStat(1111);
+ c->cd(ic++);
+ TH1 *hmass = (TH1*)l->At(3);
+ u::DrawHist(hmass, 2);
+ hmass->SetAxisRange(0.12, 0.16);
+
+ c->cd(ic++);
+ TH1 *hres = (TH1*)l->At(4);
+ u::DrawHist(hres, 2);
+ hres->SetAxisRange(0.05, 0.120);
+
+ if(ic<nx*ny) {
+ c->cd(ic++);
+ u::DrawHist((TH1*)l->At(5), 2);
+
+ c->cd(ic++);
+ u::DrawHist((TH1*)l->At(6), 2);
+ }
+ c->Update();
+}
+
+TList* AliEMCALFolder::BookHists()
+{
+ gROOT->cd();
+ TH1::AddDirectory(1);
+
+ new TH1F("00_HStat", "hist of common EMCAL statistics", 100, 0.5, 100.5);
+ new TH1F("01_EffMassAll", "effective mass of #gamma,#gamma(m_{#pi^{0}}=134.9766 MeV) - whole EMCAL", 250,0.0,0.5);
+
+ TList *l = AliEMCALHistoUtilities::MoveHistsToList("HistsOfEmcal", kFALSE);
+
+ TH1::AddDirectory(0);
+ return l;
+}
+
+void AliEMCALFolder::CreateCellNtuple()
+{
+ // Jun 28, 2007
+ if(fCellNtuple) { // Already exist
+ fCellNtuple->Print();
+ return;
+ }
+ // Create ntuple
+ Int_t bsize = int(1.e+5);
+ fCellNtuple = new TNtuple("cells","Cells Ntuple for quick analysis",
+ "fAbsId:fSupMod:fModule:fPhi:fEta:fPhiCell:fEtaCell:fCcIn:fCcOut", bsize);
+ AliEMCALCell *cell=0;
+ // Fill ntuple
+ AliEMCALSuperModule* sm = GetSuperModule(0);
+ if(sm) printf(" TNtuple was created ! sm0 %s \n", sm->GetName());
+ for(int eta=0; eta<48; eta++) { // eta row
+ TDataSet *setEta = sm->FindByName(Form("ETA%2.2i",eta));
+ if(setEta) {
+ printf(" ***** eta row %s ******\n", setEta->GetName());
+ for(int phi=0; phi<setEta->GetListSize(); phi++) { // cycle on cells (phi directions)
+ cell = (AliEMCALCell*)setEta->At(phi);
+ if(cell) {
+ cell->FillCellNtuple(fCellNtuple);
+ //printf(" fill cell %s : %s \n", cell->GetName(), cell->GetTitle());
+ }
+ }
+ }
+ }
+ fCellNtuple->Print();
+ lObj->Add(fCellNtuple);
+ // --
+ CreateAndFillAdditionalHists();
+}
+
+void AliEMCALFolder::CreateAndFillAdditionalHists()
+{
+ TH1::AddDirectory(0);
+ fLhists->Add(new TH1F("02_CCoutOnEdge", "cc out on edge of calorimeter (in MeV)", 70, 12., 19.));
+ fLhists->Add(new TH1F("03_CCoutInside", "cc out inside of calorimeter (in MeV)", 70, 12., 19.));
+ // Fill
+ Float_t* args;
+ for(Int_t i=0; i<(Int_t)fCellNtuple->GetEntries(); i++){
+ fCellNtuple->GetEvent(i);
+ args = fCellNtuple->GetArgs();
+ Int_t phi = (Int_t)args[5];
+ Int_t eta = (Int_t)args[6];
+ Double_t cc = (Double_t)args[8]*1000.;
+ if((phi==0||phi==23) || (eta==0||eta==47)) u::FillH1(fLhists, 2, cc);
+ else u::FillH1(fLhists, 3, cc);
+ }
+ // Drawing
+ Int_t wh=530, ww=750;
+ TCanvas *c = new TCanvas("c_edge","CEDGE", 10, 10, ww, wh);
+
+ gStyle->SetOptStat(1100);
+ gStyle->SetOptFit(111);
+ TH1 *h1 = (TH1*)fLhists->At(3);
+ TF1 *g = u::Gausi("ccInside", 14.7, 16.4, h1);
+ g->SetLineColor(kRed);
+ h1->Fit(g,"Q+","", 14.7, 16.4);
+ u::DrawHist(h1,2);
+ h1->SetTitle("CC distribution after #pi^{0} calibration");
+ h1->SetXTitle(" MeV ");
+ h1->SetYTitle(" N ");
+ TLatex *lat1 = u::lat(Form("rel.width = %4.2f%%",
+ 100.*h1->GetRMS()/ h1->GetMean()), 16.5, 100., 12, 0.045);
+ TLatex *lat2 = u::lat(Form("rel.width = %4.2f%% (from fit)",
+ 100.*g->GetParameter(2)/ g->GetParameter(1)), 16.5, 70., 12, 0.045);
+
+ if(0) {
+ TH1 *h2 = (TH1*)fLhists->At(2);
+ u::DrawHist(h2,2,1,"same",2);
+ }
+
+ TH1F *hccFirst = AliEMCALCalibCoefs::GetHistOfCalibTableFromDb("ccTmp");
+ u::DrawHist(hccFirst,2,1,"same",3);
+
+
+ TLegend *leg = new TLegend(0.1,0.6, 0.45,0.85);
+ leg->AddEntry(hccFirst, "Initial cc ", "L");
+ leg->AddEntry(h1, "Final cc", "L");
+ leg->Draw();
+
+ c->Update();
+}
+
+void AliEMCALFolder::TestSMStruct()
+{
+ // testing May 22, 2007
+ for(int m=0; m<12; m++) {
+ AliEMCALSuperModule *sm = new AliEMCALSuperModule(m);
+ Add(sm);
+ }
+}
+
--- /dev/null
+#ifndef ALIEMCALFOLDER_H
+#define ALIEMCALFOLDER_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $log:$ */
+
+//_________________________________________________________________________
+// Top EMCAL folder - keep everyrhing for calibration task
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+// --- ROOT system ---
+
+#include <TObjectSet.h>
+#include <TString.h>
+
+class AliEMCALGeometry;
+class AliEMCALSuperModule;
+class AliEMCALCell;
+class AliESDCaloCluster;
+class AliEMCALPi0SelectionParam;
+class pi0SelectionParam;
+class AliEMCALCalibData;
+class AliEMCALCalibCoefs;
+class AliEMCALRecPoint;
+
+class TList;
+class TNtuple;
+
+class AliEMCALFolder : public TObjectSet {
+
+ public:
+
+ AliEMCALFolder();
+ AliEMCALFolder(const char* name, const char* title="Top EMCAL folder", Bool_t putToBrowser=kFALSE);
+ AliEMCALFolder(const Int_t it, const char* title="Top EMCAL folder", Bool_t putToBrowser=kFALSE);
+
+ virtual ~AliEMCALFolder();
+
+ void Init(Bool_t putToBrowser=kFALSE);
+ // Get methods
+ Int_t GetIterationNumber() const {return fCounter;}
+ AliEMCALSuperModule* GetSuperModule(const Int_t nm);
+ TList* GetHists() {return fLhists;}
+ AliEMCALCell** GetListOfCells() {return fLofCells;}
+ AliEMCALCell* GetCell(const Int_t absId);
+ void SetCell(AliEMCALCell *cell, const Int_t absId);
+ AliEMCALPi0SelectionParam* GetPi0SelectionPar() {return fPi0SelPar;}
+ pi0SelectionParam* GetPi0SelectionParRow(Int_t nrow);
+
+ void FillPi0Candidate(const Double_t mgg, AliESDCaloCluster* cl1, AliESDCaloCluster* cl2);
+ void FillPi0Candidate(const Double_t mgg, Int_t absIdMax, Int_t nm);
+ // Define CC
+ void FitAllSMs(); // SM0 now
+ // Service routine
+ AliEMCALCalibCoefs* GetCCTable(const char* name);
+ AliEMCALCalibCoefs* GetCCFirst() {return GetCCTable(fgkCCFirstName.Data());}
+ AliEMCALCalibCoefs* GetCCIn() {return GetCCTable(fgkCCinName.Data());}
+ AliEMCALCalibCoefs* GetCCOut(){return GetCCTable(fgkCCoutName.Data());}
+ Int_t GetSMNumber(AliESDCaloCluster* cl);
+ // Recalibration staf - Jun 18,2007
+ static AliEMCALRecPoint *GetRecPoint(AliESDCaloCluster *cl,AliEMCALCalibCoefs *tOld,AliEMCALCalibCoefs *tNew, TList *l=0);
+ // MENU
+ void Save(const char *fn = "EMCALFOLDER.root", const char *opt="RECREATE"); // *MENU*
+ static AliEMCALFolder* Read(const char *fn = "EMCALFOLDER.root", const char *opt="READ");
+ void InitAfterRead(); // *MENU*
+ void DrawQA(const int nsm=0); // *MENU*
+ void CreateCellNtuple(); // *MENU*
+ void CreateAndFillAdditionalHists(); // *MENU*
+
+ protected:
+ TList* BookHists();
+ Int_t fCounter; // Coonter of iteration
+ //
+ AliEMCALGeometry *fGeometry; //
+ //
+ Int_t fNumOfCell; // number of cells as in geometry
+
+ TList* fLhists; //! for speed
+ AliEMCALCell** fLofCells; //! unifrom array of cells for fast access; invisible from browser
+ AliEMCALPi0SelectionParam* fPi0SelPar; // pi0 selection parameters
+ AliEMCALCalibData *fCalibData; //!
+ //
+ TNtuple *fCellNtuple; //! for quick cell anaylsis
+
+ public:
+ static const TString fgkBaseFolderName; // base name of EMCAL Folder
+ static const TString fgkCCFirstName; // name of first calib.table
+ static const TString fgkCCinName; // name of initial calib.coefs. table
+ static const TString fgkCCoutName; // name of out calib.coefs. table
+
+ void TestSMStruct();
+
+ ClassDef(AliEMCALFolder,2) // EMCAL folder
+
+};
+
+#endif // ALIEMCALFOLDER_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2007, 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. *
+ **************************************************************************/
+
+/* $Log$ */
+
+//_________________________________________________________________________
+// Set of parameters for pi0 selection
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+#include "AliEMCALPi0SelectionParam.h"
+
+TableClassImpl(AliEMCALPi0SelectionParam,pi0SelectionParam)
+
+void AliEMCALPi0SelectionParam::PrintTable()
+{
+ printf(" Table : %s : nrows %i \n", GetName(), int(GetNRows()));
+ for(int i=0; i<GetNRows(); i++) PrintTable(i);
+}
+
+void AliEMCALPi0SelectionParam::PrintTable(const Int_t i)
+{
+ if(i>=GetNRows()) return;
+ printf("row %i \n", i);
+ PrintRec(GetTable(i));
+}
+
+void AliEMCALPi0SelectionParam::PrintRec(pi0SelectionParam* r)
+{
+ if(r==0) return;
+ printf(" cluster energy window %7.2f -> %7.2f \n", r->eOfRpMin, r->eOfRpMax);
+ printf(" gamma,gamma mass window %7.2f -> %7.2f \n", r->massGGMin, r->massGGMax);
+ printf(" pi0 momentum window %7.2f -> %7.2f \n", r->momPi0Min, r->momPi0Max);
+}
+// Set 1;
+AliEMCALPi0SelectionParam* AliEMCALPi0SelectionParam::Set1()
+{
+ pi0SelectionParam r;
+ r.eOfRpMin = 0.3;
+ r.eOfRpMax = 30.;
+ r.massGGMin = 0.03;
+ r.massGGMax = 0.28;
+ r.momPi0Min = 1.8;
+ r.momPi0Max = 12.0;
+ AliEMCALPi0SelectionParam *t = new AliEMCALPi0SelectionParam("Pi0Set1",1);
+ t->AddAt(&r);
+ return t;
+}
--- /dev/null
+#ifndef ALIEMCALPI0SELECTIONPARAM_H
+#define ALIEMCALPI0SELECTIONPARAM_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+// Set of parameters for pi0 selection
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+// --- ROOT system ---
+#include <TTable.h>
+
+// unit is GeV
+struct pi0SelectionParam {
+ double eOfRpMin; // minimal energy of em.cluster (rec point)
+ double eOfRpMax; // maximal energy of em.cluster (rec point)
+ double massGGMin; // minimal mass of gamma,gamma
+ double massGGMax; // maximal mass of gamma,gamma
+ double momPi0Min; // minimal pi0 momentum
+ double momPi0Max; // maximal pi0 momentum
+};
+
+class AliEMCALPi0SelectionParam : public TTable {
+ public:
+ ClassDefTable(AliEMCALPi0SelectionParam , pi0SelectionParam)
+
+ // Menu
+ void PrintTable(); // *MENU*
+ void PrintTable(const Int_t i); // *MENU*
+ void PrintRec(pi0SelectionParam *r);
+
+ // Set of parameter(s)
+ static AliEMCALPi0SelectionParam* Set1();
+
+ ClassDef(AliEMCALPi0SelectionParam,1) // Set of Parameters For Pi0 Selection
+};
+
+#endif // ALIEMCALPI0SELECTIONPARAM_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2002, 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. *
+ **************************************************************************/
+
+/* $Log$ */
+
+//_________________________________________________________________________
+// This is selector for making a set of histogramms from ESD for rec.points
+//*-- Authors: Aleksei Pavlinov (WSU)
+//*-- Feb 17, 2007 - May 23, 2007
+//*-- Pi0 calibration
+//_________________________________________________________________________
+
+#include "AliEMCALRecPointsQaESDSelector.h"
+#include "AliEMCALFolder.h"
+#include "AliEMCALSuperModule.h"
+#include "AliEMCALCell.h"
+#include "AliEMCALPi0SelectionParam.h"
+#include "AliEMCALHistoUtilities.h"
+#include "AliEMCALCalibCoefs.h"
+#include "AliEMCALGeometry.h"
+#include "AliLog.h"
+#include "AliESD.h"
+#include "AliEMCALRecPoint.h"
+
+#include <assert.h>
+#include <map>
+#include <string>
+
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TSystem.h>
+#include <TCanvas.h>
+#include <TVector3.h>
+#include <TLorentzVector.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <TArrayS.h>
+#include <TBrowser.h>
+#include <TDataSet.h>
+
+using namespace std;
+
+typedef AliEMCALHistoUtilities u;
+
+ClassImp(AliEMCALRecPointsQaESDSelector)
+
+Int_t nmaxCell = 4*12*24*11;
+
+AliEMCALRecPointsQaESDSelector::AliEMCALRecPointsQaESDSelector() :
+ AliSelector(),
+ fChain(0),
+ fLofHistsPC(0),
+ fLofHistsRP(0),
+ fEmcalPool(0),
+ fEMCAL(0),
+ fEMCALOld(0)
+{
+ //
+ // Constructor. Initialization of pointers
+ //
+}
+
+AliEMCALRecPointsQaESDSelector::~AliEMCALRecPointsQaESDSelector()
+{
+ //
+ // Destructor
+ //
+
+ // histograms are in the output list and deleted when the output
+ // list is deleted by the TSelector dtor
+}
+
+void AliEMCALRecPointsQaESDSelector::Begin(TTree* tree)
+{
+ // Begin function
+
+ AliSelector::Begin(tree);
+}
+
+void AliEMCALRecPointsQaESDSelector::Init(TTree* tree)
+{
+ // Init function
+
+ AliSelector::Init(tree);
+}
+
+void AliEMCALRecPointsQaESDSelector::SlaveBegin(TTree* tree)
+{
+ // The SlaveBegin() function is called after the Begin() function.
+ // When running with PROOF SlaveBegin() is called on each slave server.
+ // The tree argument is deprecated (on PROOF 0 is passed).
+
+ printf("**** <I> Slave Begin \n **** ");
+
+ AliSelector::SlaveBegin(tree);
+}
+
+void AliEMCALRecPointsQaESDSelector::InitStructure(Int_t it)
+{
+ AliEMCALGeometry::GetInstance("SHISH_TRD1_CURRENT_2X2"); // initialize geometry just once
+ //AliEMCALGeometry::GetInstance(""); // initialize geometry just once
+
+ fLofHistsPC = DefineHistsOfRP("PseudoCl");
+ fLofHistsRP = DefineHistsOfRP("RP");
+
+ fEmcalPool = new TObjectSet("PoolOfEMCAL");
+ if(it == 1) {
+ fEMCAL = new AliEMCALFolder(it); // folder for first itteration
+ fEmcalPool->Add(fEMCAL);
+ }
+}
+
+Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
+{
+ // The Process() function is called for each entry in the tree (or possibly
+ // keyed object in the case of PROOF) to be processed. The entry argument
+ // specifies which entry in the currently loaded tree is to be processed.
+ // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
+ // to read either all or the required parts of the data. When processing
+ // keyed objects with PROOF, the object is already loaded and is available
+ // via the fObject pointer.
+ //
+ // This function should contain the "body" of the analysis. It can contain
+ // simple or elaborate selection criteria, run algorithms on the data
+ // of the event and typically fill histograms.
+
+ // WARNING when a selector is used with a TChain, you must use
+ // the pointer to the current TTree to call GetEntry(entry).
+ // The entry is always the local entry number in the current tree.
+ // Assuming that fTree is the pointer to the TChain being processed,
+ // use fTree->GetTree()->GetEntry(entry).
+
+ if (AliSelector::Process(entry) == kFALSE)
+ return kFALSE;
+
+ // Check prerequisites
+ if (!fESD)
+ {
+ AliDebug(AliLog::kError, "ESD branch not available");
+ return kFALSE;
+ }
+ static pi0SelectionParam* rPar = GetEmcalFolder()->GetPi0SelectionParRow(0);
+
+ static Int_t nEmcalClusters, indOfFirstEmcalRP, nEmcalRP,nEmcalPseudoClusters;
+ nEmcalClusters = fESD->GetNumberOfEMCALClusters();
+ indOfFirstEmcalRP = fESD->GetFirstEMCALCluster();
+ u::FillH1(fLofHistsRP, 1, double(indOfFirstEmcalRP));
+
+ static AliESDCaloCluster *cl = 0;
+ nEmcalRP = nEmcalPseudoClusters = 0;
+ TList *l=0;
+ double eDigi=0;
+
+ TClonesArray lvM1("TLorentzVector", 100); // for convenience; 100 is enough now
+ TArrayI indLv(100); // index of RP
+
+ static TLorentzVector v, vcl;
+ int nrp = 0; // # of RP for gg analysis
+ for(int i=indOfFirstEmcalRP; i<indOfFirstEmcalRP+nEmcalClusters; i++) {
+ cl = fESD->GetCaloCluster(i);
+ if(cl->GetClusterType() == AliESDCaloCluster::kPseudoCluster) {
+ nEmcalPseudoClusters++;
+ l = fLofHistsPC;
+ } else if(cl->GetClusterType() == AliESDCaloCluster::kClusterv1){
+ nEmcalRP++;
+ if(fEMCAL->GetIterationNumber()>1) {
+ AliEMCALRecPoint *rp = AliEMCALFolder::GetRecPoint(cl, fEMCAL->GetCCFirst(), fEMCAL->GetCCIn(), fLofHistsRP);
+ // if(rp->GetPointEnergy()>=rPar->eOfRpMin && u::GetLorentzVectorFromRecPoint(v, rp)) {
+ if(u::GetLorentzVectorFromRecPoint(v, rp)) { // comparing with RP
+ new(lvM1[nrp]) TLorentzVector(v);
+ indLv[nrp] = i;
+ nrp++;
+ // Conroling of recalibration
+ u::GetLorentzVectorFromESDCluster(vcl, cl);
+ u::FillH1(fLofHistsRP, 11, vcl.P()-v.P());
+ u::FillH1(fLofHistsRP, 12, TMath::RadToDeg()*(vcl.Theta()-v.Theta()));
+ u::FillH1(fLofHistsRP, 13, TMath::RadToDeg()*(vcl.Phi()-v.Phi()));
+ l = 0; // no filling
+ }
+ delete rp;
+ } else { // first iteration
+ // if(cl->E()>=rPar->eOfRpMin && u::GetLorentzVectorFromESDCluster(v, cl)) {
+ if(u::GetLorentzVectorFromESDCluster(v, cl)) { // comparing with RP
+ // cut 0.4 GeV may be high !
+ new(lvM1[nrp]) TLorentzVector(v);
+ indLv[nrp] = i;
+ nrp++;
+ l = fLofHistsRP;
+ }
+ }
+
+ } else {
+ printf(" wrong cluster type : %i\n", cl->GetClusterType());
+ assert(0);
+ }
+ u::FillH1(l, 2, double(cl->GetClusterType()));
+
+ u::FillH1(l, 3, double(cl->GetNumberOfDigits()));
+ u::FillH1(l, 4, double(cl->E()));
+ // Cycle on digits (towers)
+ Short_t *digiAmpl = cl->GetDigitAmplitude()->GetArray();
+ Short_t *digiTime = cl->GetDigitTime()->GetArray();
+ Short_t *digiAbsId = cl->GetDigitIndex()->GetArray();
+ for(int id=0; id<cl->GetNumberOfDigits(); id++) {
+ eDigi = double(digiAmpl[id]) / 500.; // See AliEMCALClusterizerv1
+ // if(eDigi <= 0.0) { // sometimes it is happen
+ //if(eDigi > 10.0 && cl->GetClusterType() == AliESDCaloCluster::kClusterv1) {
+ // printf(" %i digiAmpl %i : %f \n", id, int(digiAmpl[id]), eDigi);
+ //}
+ u::FillH1(l, 5, eDigi);
+ u::FillH1(l, 6, double(digiTime[id]));
+ u::FillH1(l, 7, double(digiAbsId[id]));
+ if(int(digiAbsId[id]) >= nmaxCell) {
+ printf(" id %i : digiAbsId[id] %i (%i) : %s \n",
+ id, int(digiAbsId[id]), nmaxCell, l->GetName());
+ }
+ }
+ }
+ u::FillH1(fLofHistsRP, 0, double(nEmcalRP));
+ u::FillH1(fLofHistsPC, 0, double(nEmcalPseudoClusters));
+
+ static TLorentzVector *lv1=0, *lv2=0, lgg;
+ static Double_t mgg, pgg;
+ mgg = pgg = 0.;
+ nrp = lvM1.GetEntriesFast();
+ if(nrp >= 2) {
+ // eff.mass analysis
+ for(int i1=0; i1<nrp-1; i1++){
+ lv1 = (TLorentzVector*)lvM1.At(i1);
+ for(int i2=i1+1; i2<nrp; i2++){
+ lv2 = (TLorentzVector*)lvM1.At(i2);
+ lgg = (*lv1) + (*lv2);
+ mgg = lgg.M(); // eff.mass
+ pgg = lgg.P(); // momentum
+ u::FillH1(fLofHistsRP, 8, mgg);
+
+ if((mgg>=rPar->massGGMin && mgg<=rPar->massGGMax)) {// pi0 candidates
+ if((pgg>=rPar->momPi0Min && pgg>=rPar->momPi0Min)) {
+ fEMCAL->FillPi0Candidate(mgg,fESD->GetCaloCluster(indLv[i1]),fESD->GetCaloCluster(indLv[i2]));
+ u::FillH1(fLofHistsRP, 9, pgg);
+ u::FillH1(fLofHistsRP,10, lv1->P());
+ u::FillH1(fLofHistsRP,10, lv2->P());
+ }
+ }
+ }
+ }
+ }
+
+ lvM1.Delete();
+
+ if(nEmcalClusters != (nEmcalRP+nEmcalPseudoClusters))
+ Info("Process","nEmcalClusters %i : RP %i + PC %i ",nEmcalClusters, nEmcalRP, nEmcalPseudoClusters);
+
+ return kTRUE;
+}
+
+void AliEMCALRecPointsQaESDSelector::SlaveTerminate()
+{
+ // The SlaveTerminate() function is called after all entries or objects
+ // have been processed. When running with PROOF SlaveTerminate() is called
+ // on each slave server.
+
+ AliSelector::SlaveTerminate();
+
+ // Add the histograms to the output on each slave server
+ if (!fOutput)
+ {
+ AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
+ return;
+ }
+
+ // fOutput->Add(fMultiplicity);
+}
+
+void AliEMCALRecPointsQaESDSelector::Terminate()
+{
+ // The Terminate() function is the last function to be called during
+ // a query. It always runs on the client, it can be used to present
+ // the results graphically or save the results to file.
+
+ AliSelector::Terminate();
+ /*
+ fMultiplicity = dynamic_cast<TH1F*> (fOutput->FindObject("fMultiplicity"));
+
+ if (!fMultiplicity)
+ {
+ AliDebug(AliLog::kError, Form("ERROR: Histogram not available %p", (void*) fMultiplicity));
+ return;
+ }
+
+ new TCanvas;
+ fMultiplicity->DrawCopy();
+ */
+ PrintInfo();
+}
+//
+TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfRP(const char *name)
+{
+ double ADCchannelEC = 0.00305; // ~3mev per adc count
+
+ gROOT->cd();
+ TH1::AddDirectory(1);
+ new TH1F("00_EmcalMultiplicity", "multiplicity of EMCAL ", 201, -0.5, 200.5); // real and pseudo
+ new TH1F("01_IndexOfFirstEmcal", "index of first emcal rec.points ", 201, -0.5, 200.5);
+
+ new TH1F("02_NumberOf", "number of ", 6, -0.5, 5.5);
+ new TH1F("03_NumberOfDigitsIn", "number of digits(towers) in rec.points ", 101, -0.5, 100.5);
+ new TH1F("04_EnergyOf", "energy of ", 1000, 0.0, 100.);
+
+ int nmax=10000;
+ double xmi = ADCchannelEC/2., xma = xmi + ADCchannelEC*nmax;
+ // All energy(momentum) unit is GeV if don't notice
+ new TH1F("05_DigitEnergyIn", "digit energy in ", nmax, xmi, xma);
+ new TH1F("06_DigitTimeIn", "digit time in 10ps(0.01ns) ", 1000, 0.0, 3.e+3); // ns/100 = 10 ps
+ new TH1F("07_DigitAbsIdIn", "digit abs id in ", nmaxCell, -0.5, double(nmaxCell)-0.5);
+ new TH1F("08_EffMass", "effective mass of #gamma,#gamma(m_{#pi^{0}}=134.9766 MeV)", 100, 0.0, 0.5);
+ new TH1F("09_MomOfPi0Candidate", "momentum of #pi^{0} candidates (0.085 <mgg<0.185)", 600, 0.0, 30.0);
+ new TH1F("10_MomOfRpPi0Candidate", "momentum of RP for #pi^{0} candidates (0.085 <mgg<0.185)", 600, 0.0, 30.0);
+ // Recalibration staf
+ new TH1F("11_MomClESD-RpRecalib", "difference of momentum cl(ESD) - rp(Recalib)", 200, -1., +1.);
+ new TH1F("12_ThetaClESD-RpRecalib", "difference of #theta cl(ESD) - rp(Recalib) in degree", 100, -0.05, +0.05);
+ new TH1F("13_PhiClESD-RpRecalib", "difference of #phi cl(ESD) - rp(Recalib) in degree ", 100, -0.05, +0.05);
+ // Digi
+ new TH1F("14_EDigiRecalib", "energy of digits after recalibration", 2000, 0.0, 20.);
+ AliEMCALGeometry* g = AliEMCALGeometry::GetInstance();
+ new TH1F("15_AbsIdRecalib", "abs Id of digits after recalibration", g->GetNCells(),-0.5,Double_t(g->GetNCells())-0.5);
+
+ TList *l = u::MoveHistsToList(Form("ListOfHists%s",name), kFALSE);
+ u::AddToNameAndTitleToList(l, name, name);
+
+ return l;
+}
+
+/* unused now
+TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfTowers(const char *name)
+{
+ //
+ // ESD: towers information was saved to pseudo clusters
+ //
+ gROOT->cd();
+ TH1::AddDirectory(1);
+ new TH1F("00_EmcalMultiplicity", "multiplicity of EMCAL ", 201, -0.5, 200.5); // number of pseudo RP
+
+ new TH1F("01_EnergyOf", "energy of ", 1000, 0.0, 100.);
+
+ TList *l = u::MoveHistsToList(Form("ListOfHists%s",name, " - ESD"), kFALSE);
+ u::AddToNameAndTitleToList(l, name, name);
+
+ return l;
+}
+*/
+
+void AliEMCALRecPointsQaESDSelector::FitEffMassHist()
+{
+ TH1* h = (TH1*)fLofHistsRP->At(8);
+ AliEMCALCell::FitHist(h, GetName(), "draw");
+}
+
+void AliEMCALRecPointsQaESDSelector::PrintInfo()
+{
+ // Service routine
+ TList *l[2] = {fLofHistsPC, fLofHistsRP};
+ printf("\n");
+ for(int i=0; i<2; i++){
+ TH1F *h = (TH1F*)l[i]->At(2);
+ printf(" %s \t: %i \n", h->GetTitle(), int(h->GetEntries()));
+ }
+}
+
+AliEMCALFolder* AliEMCALRecPointsQaESDSelector::CreateEmcalFolder(const Int_t it)
+{
+ AliEMCALFolder* newFolder = new AliEMCALFolder(it); // folder for iteration #it
+ if(it>1) {
+ fEMCALOld = fEMCAL;
+ AliEMCALCalibCoefs* tabOldOut = fEMCALOld->GetCCOut();
+ AliEMCALCalibCoefs* tabNewIn = new AliEMCALCalibCoefs(*tabOldOut);
+ tabNewIn->SetName(AliEMCALFolder::fgkCCinName.Data());
+ newFolder->Add(tabNewIn);
+ }
+ fEmcalPool->Add(newFolder);
+ fEMCAL = newFolder;
+
+ return fEMCAL;
+}
+
+AliEMCALFolder* AliEMCALRecPointsQaESDSelector::GetEmcalOldFolder(const Int_t nsm)
+{
+ AliEMCALFolder* folder=0;
+ if(fEmcalPool) folder = (AliEMCALFolder*)fEmcalPool->FindByName(Form("EMCAL_%2.2i",nsm));
+ return folder;
+}
+
+
+void AliEMCALRecPointsQaESDSelector::SetEmcalFolder(AliEMCALFolder* folder)
+{
+ fEMCAL = folder;
+ fEmcalPool->Add(fEMCAL);
+}
+
+void AliEMCALRecPointsQaESDSelector::SetEmcalOldFolder(AliEMCALFolder* folder)
+{
+ fEMCALOld = folder;
+ fEmcalPool->Add(fEMCALOld);
+}
+
+void AliEMCALRecPointsQaESDSelector::Browse(TBrowser* b)
+{
+ if(fESD) b->Add(fESD);
+ if(fLofHistsPC) b->Add(fLofHistsPC);
+ if(fLofHistsRP) b->Add(fLofHistsRP);
+ if(fChain) b->Add(fChain);
+ if(fEmcalPool) b->Add(fEmcalPool);
+ // if(u) b->Add(u);
+}
+
+Bool_t AliEMCALRecPointsQaESDSelector::IsFolder() const
+{
+ if(fESD || fLofHistsRP || fEmcalPool) return kTRUE;
+ return kFALSE;
+}
+
+void AliEMCALRecPointsQaESDSelector::ReadAllEmcalFolders()
+{
+ if(fEmcalPool==0) {
+ fEmcalPool = new TObjectSet("PoolOfEMCAL");
+ for(Int_t it=1; it<=10; it++){
+ AliEMCALFolder* fold = AliEMCALFolder::Read(Form("EMCALFOLDER_It%i_fit.root",it), "READ");
+ if(fold) fEmcalPool->Add(fold);
+ }
+ }
+}
+
+void AliEMCALRecPointsQaESDSelector::PictVsIterNumber(const Int_t ind, const Int_t nsm)
+{
+ // Jun 27, 2007 - unfinished; which picture is the best
+ if(ind<0 || ind>5) return;
+ gROOT->cd();
+ TH1::AddDirectory(1);
+
+ Int_t itMax = 10, it=0;
+ map <int, char*> indName;
+ indName[0] = "eff.mass";
+ indName[3] = "mass of #pi_{0}";
+ indName[4] = "resolution of #pi_{0}";
+ indName[5] = "chi^{2}/ndf";
+
+ TH1F *hout = new TH1F(indName[ind], indName[ind], itMax, 0.5, double(itMax)+0.5);
+
+ TH1::AddDirectory(0);
+ Double_t content, error;
+ for(Int_t i=0; i<fEmcalPool->GetListSize(); i++) { // cycle on EMCAL folders
+ AliEMCALFolder* folder = static_cast<AliEMCALFolder*>(fEmcalPool->At(i));
+ if(folder==0) continue;
+ it = folder->GetIterationNumber();
+
+ AliEMCALSuperModule* sm = folder->GetSuperModule(nsm);
+ if(sm==0) continue;
+
+ TList* l = sm->GetHists();
+ if(l==0) continue;
+
+ TH1F *hin = (TH1F*)l->At(ind);
+ if(hin==0) continue;
+
+ if(ind !=0 ) {
+ content = hin->GetMean();
+ error = hin->GetRMS();
+ } else {
+ sm->FitEffMassHist();
+ TF1 *f = (TF1*)hin->GetListOfFunctions()->At(0);
+ content = error = -1.;
+ if(f) {
+ // content = f->GetParameter(1);
+ //error = f->GetParameter(2);
+ content = f->GetParameter(2);
+ error = f->GetParError(2);
+ }
+ }
+
+ if(content > 0.0) {
+ hout->SetBinContent(it, content);
+ hout->SetBinError(it, error);
+ printf(" it %i content %f +/- %f \n", it, content, error);
+ }
+ }
+
+ u::DrawHist(hout,2);
+ hout->SetMinimum(0.0);
+
+}
--- /dev/null
+#ifndef ALIEMCALRECPOINTSQAESDSELECTOR_H
+#define ALIEMCALRECPOINTSQAESDSELECTOR_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+//*-- Authors: Aleksei Pavlinov (WSU)
+
+#include "AliSelector.h"
+
+class AliEMCALFolder;
+
+class TList;
+class TBrowser;
+class TChain;
+class TObjectSet;
+
+class AliEMCALRecPointsQaESDSelector : public AliSelector {
+ public:
+ AliEMCALRecPointsQaESDSelector();
+ virtual ~AliEMCALRecPointsQaESDSelector();
+
+ virtual void Begin(TTree* tree);
+ virtual void Init(TTree* tree);
+ virtual void SlaveBegin(TTree *tree);
+ virtual Bool_t Process(Long64_t entry);
+ virtual void SlaveTerminate();
+ virtual void Terminate();
+ //
+ void InitStructure(Int_t it);
+ static TList *DefineHistsOfRP(const char *name="RP");
+ // static TList *DefineHistsOfTowers(const char *name="towers");
+ //
+ void FitEffMassHist(); //*MENU*
+ void PrintInfo(); // *MENU*
+ //
+ void SetChain(TChain *chain) {fChain = chain;}
+ TChain* GetChain() {return fChain;}
+
+ AliEMCALFolder* CreateEmcalFolder(const Int_t it);
+ void SetEmcalFolder(AliEMCALFolder* folder);
+ void SetEmcalOldFolder(AliEMCALFolder* folder);
+ AliEMCALFolder* GetEmcalFolder() {return fEMCAL;}
+ AliEMCALFolder* GetEmcalOldFolder() {return fEMCALOld;}
+ AliEMCALFolder* GetEmcalOldFolder(const Int_t nsm);
+ //
+ virtual void Browse(TBrowser* b);
+ virtual Bool_t IsFolder() const;
+
+ //
+ //// Pictures staf - Jun 26, 2007
+ //
+ void ReadAllEmcalFolders();
+ void PictVsIterNumber(const Int_t ind=0, const Int_t nsm=0);
+
+ protected:
+ //
+ TChain* fChain; //! chain if ESD files
+ TList* fLofHistsPC; //! list of histograms of pseudo clusters
+ TList* fLofHistsRP; //! list of histograms of rec.points
+ //
+ TObjectSet* fEmcalPool;
+ AliEMCALFolder* fEMCAL; //! current EMCAL object
+ AliEMCALFolder* fEMCALOld; //! previous EMCAL object
+
+ private:
+ AliEMCALRecPointsQaESDSelector(const AliEMCALRecPointsQaESDSelector&);
+ AliEMCALRecPointsQaESDSelector& operator=(const AliEMCALRecPointsQaESDSelector&);
+
+ ClassDef(AliEMCALRecPointsQaESDSelector, 1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2007, 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$ */
+
+//_________________________________________________________________________
+// Top EMCAL folder which will keep all information about EMCAL itself,
+// super Modules (SM), modules, towers, set of hists and so on.
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+#include "AliEMCALSuperModule.h"
+#include "AliEMCALFolder.h"
+#include "AliEMCALCell.h"
+#include "AliEMCALHistoUtilities.h"
+
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TList.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+
+typedef AliEMCALHistoUtilities u;
+
+ClassImp(AliEMCALSuperModule)
+
+AliEMCALSuperModule::AliEMCALSuperModule() : TObjectSet(),
+fSMNumber(0)
+{
+}
+
+AliEMCALSuperModule::AliEMCALSuperModule(const Int_t m, const char* title) :
+TObjectSet(Form("SM%2.2i",m)),
+fSMNumber(m)
+{
+ SetTitle(title);
+}
+
+AliEMCALSuperModule::~AliEMCALSuperModule()
+{
+ // dtor
+}
+
+void AliEMCALSuperModule::Init()
+{
+ if(GetHists()==0) this->AddObject(BookHists(), kTRUE);
+}
+
+void AliEMCALSuperModule::AddCellToEtaRow(AliEMCALCell *cell, const Int_t etaRow)
+{
+ static TDataSet *set;
+ set = FindByName(Form("ETA%2.2i",etaRow));
+ if(set==0) {
+ set = new TDataSet(Form("ETA%2.2i",etaRow));
+ Add(set);
+ }
+ set->Add(cell);
+}
+
+void AliEMCALSuperModule::FitForAllCells()
+{
+ Int_t ncells=0;
+ for(int eta=0; eta<48; eta++) { // eta row
+ TDataSet *setEta = FindByName(Form("ETA%2.2i",eta));
+ if(setEta) {
+ printf(" eta %i : %s : cells %i \n", eta, setEta->GetName(), setEta->GetListSize());
+ for(int phi=0; phi<setEta->GetListSize(); phi++) { // cycle on cells (phi directions)
+ AliEMCALCell* cell = (AliEMCALCell*)setEta->At(phi);
+ cell->FitEffMassHist();
+
+ u::FillH1(GetHists(), 1, cell->GetCcIn()*1.e+3);
+ u::FillH1(GetHists(), 2, cell->GetCcOut()*1.e+3);
+
+ TF1 *f = cell->GetFunction();
+ u::FillH1(GetHists(), 3, f->GetParameter(1));
+ u::FillH1(GetHists(), 4, f->GetParameter(2));
+ u::FillH1(GetHists(), 5, f->GetChisquare()/f->GetNDF());
+ u::FillH1(GetHists(), 6, f->GetParameter(0));
+
+ ncells++;
+ }
+ }
+ }
+ printf(" <I> AliEMCALSuperModule::FitForAllCells() : ncells %i with fit \n", ncells);
+}
+
+void AliEMCALSuperModule::FitEffMassHist()
+{
+ TH1* h = (TH1*)GetHists()->At(0);
+ AliEMCALCell::FitHist(h, GetName());
+}
+
+
+void AliEMCALSuperModule::PrintInfo()
+{
+ printf(" Super Module : %s : %i \n", GetName(), fSMNumber);
+ printf(" # of active cells %i \n", GetNumberOfCells());
+ TH1* h = (TH1*)GetHists()->At(0);
+ printf(" # h->Integral() of eff.mass hist %i \n", int(h->Integral()));
+}
+
+void AliEMCALSuperModule::DrawCC(int iopt)
+{
+ TCanvas *c=0;
+ if(iopt==1) c = new TCanvas("ccInOut","ccInOut");
+
+ gStyle->SetOptStat(0);
+
+ TH1 *hCCIn = (TH1*)GetHists()->At(1);
+ TH1 *hCCOut = (TH1*)GetHists()->At(2);
+
+ hCCIn->SetStats(kFALSE);
+ hCCOut->SetStats(kFALSE);
+ hCCOut->SetTitle("CC in and out");
+ hCCOut->SetXTitle("cc in MeV");
+ hCCOut->SetYTitle(" N ");
+
+ u::DrawHist(hCCOut,2);
+ hCCOut->SetAxisRange(10., 24.);
+ u::DrawHist(hCCIn,2, kRed, "same");
+
+ TLegend *leg = new TLegend(0.5,0.36, 0.97,0.80);
+ TLegendEntry *leIn = leg->AddEntry(hCCIn, Form("input cc : %6.2f #pm %6.2f", hCCIn->GetMean(),hCCIn->GetRMS()), "L");
+ leIn->SetTextColor(hCCIn->GetLineColor());
+ leg->AddEntry(hCCOut, Form("output cc : %6.2f #pm %6.2f", hCCOut->GetMean(),hCCOut->GetRMS()), "L");
+ leg->Draw();
+
+ if(c) c->Update();
+}
+
+Int_t AliEMCALSuperModule::GetNumberOfCells()
+{
+ Int_t ncells=0;
+ for(int eta=0; eta<GetListSize(); eta++) { // cycle on eta row
+ TDataSet *set = At(eta);
+ ncells += set->GetListSize();
+ }
+ return ncells;
+}
+
+TList* AliEMCALSuperModule::BookHists()
+{
+ gROOT->cd();
+ TH1::AddDirectory(1);
+
+ AliEMCALFolder* EMCAL = (AliEMCALFolder*)GetParent();
+ Int_t it = EMCAL->GetIterationNumber();
+
+ new TH1F("00_EffMass", "effective mass of #gamma,#gamma(m_{#pi^{0}}=134.98 MeV) ", 250,0.0,0.5);
+ new TH1F("01_CCInput", "input CC dist.(MEV) ", 200, 5., 25.);
+ new TH1F("02_CCOutput", "output CC dist.(MEV) ", 200, 5., 25.);
+ new TH1F("03_MPI0", "mass of #pi_{0} dist. ", 170, 0.05, 0.22);
+ new TH1F("04_RESPI0", "resolution at #pi_{0} dist. ", 50, 0.0, 0.05);
+ new TH1F("05_XI2/NDF", "#chi^{2} / ndf", 50, 0.0, 5.0);
+ new TH1F("06_NPI0", "number of #pi_{0}", 150, 0.0, 1500.);
+
+ TList *l = AliEMCALHistoUtilities::MoveHistsToList(Form("HistsOfSM_%2.2i",fSMNumber), kFALSE);
+ AliEMCALHistoUtilities::AddToNameAndTitleToList(l, Form("_%2.2i_It%i",fSMNumber, it),
+ Form(" SM %2.2i, Iter %i",fSMNumber, it));
+
+ TH1::AddDirectory(0);
+ return l;
+}
--- /dev/null
+#ifndef ALIEMCALSUPERMODULE_H
+#define ALIEMCALSUPERMODULE_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+// Emcal Super Module
+//
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+
+// --- ROOT system ---
+
+#include <TObjectSet.h>
+
+class TList;
+class AliEMCALCell;
+
+class AliEMCALSuperModule : public TObjectSet {
+
+ public:
+
+ AliEMCALSuperModule();
+ AliEMCALSuperModule(const Int_t m, const char* title="Emcal Super Module");
+
+ virtual ~AliEMCALSuperModule();
+
+ void Init();
+ void AddCellToEtaRow(AliEMCALCell *cell, const Int_t etaRow);
+ TList* GetHists() {return (TList*)fObj;}
+ // MENU
+ void FitForAllCells(); //*MENU*
+ void FitEffMassHist(); //*MENU*
+ void PrintInfo(); //*MENU*
+ void DrawCC(int iopt=1); //*MENU*
+ //
+ Int_t GetNumberOfCells();
+ Int_t GetSMNumber() const {return fSMNumber;}
+ protected:
+ TList* BookHists();
+ //
+ Int_t fSMNumber;
+
+ ClassDef(AliEMCALSuperModule,1) // EMCAL SuperModule
+
+};
+
+#endif // ALIEMCALSUPERMODULE_H
#pragma link C++ class AliEMCALCalibHistoProducer+;
#pragma link C++ class AliEMCALPreprocessor+;
#pragma link C++ class AliEMCALRawUtils+;
-
+// Calibration staff
+#pragma link C++ class AliEMCALFolder;
+#pragma link C++ class AliEMCALSuperModule;
+#pragma link C++ class AliEMCALCell;
+// Tables
+#pragma link C++ class AliEMCALPi0SelectionParam-;
+#pragma link C++ class pi0SelectionParam+;
+#pragma link C++ class AliEMCALCalibCoefs-;
+#pragma link C++ class calibCoef+;
#endif
#pragma link C++ class AliEMCALTrack+;
#pragma link C++ class AliEMCALTracker+;
#pragma link C++ class AliEMCALPID+;
+#pragma link C++ class AliEMCALRecPointsQaESDSelector+;
#pragma link C++ class AliEMCALRecParam+;
#endif
AliEMCALCalibData.cxx \
AliEMCALCalibHistoProducer.cxx \
AliEMCALPreprocessor.cxx \
-AliEMCALRawUtils.cxx
+AliEMCALRawUtils.cxx \
+AliEMCALFolder.cxx \
+AliEMCALSuperModule.cxx \
+AliEMCALCell.cxx \
+AliEMCALPi0SelectionParam.cxx\
+AliEMCALCalibCoefs.cxx
HDRS= $(SRCS:.cxx=.h)
AliEMCALTrack.cxx \
AliEMCALTracker.cxx \
AliEMCALPID.cxx \
+AliEMCALRecPointsQaESDSelector.cxx \
AliEMCALRecParam.cxx
HDRS= $(SRCS:.cxx=.h)
-
// Script to create calibration parameters and store them into CDB
-// Two sets of calibration parameters can be created:
+// Thre sets of calibration parameters can be created:
// 1) equal parameters
// 2) randomly distributed parameters for decalibrated detector silumations
+// 3) gausian distributed parameters for decalibrated detector silumations
+//
// Modified from PHOS script for EMCAL by Gustavo Conesa
+// May 3, 2007: Modified by Aleksei Pavlinov
+
+//.x $ALICE_ROOT/EMCAL/macros/CalibrationDB/AliEMCALSetCDB.C
#if !defined(__CINT__)
-#include "TControlBar.h"
-#include "TString.h"
-#include "TRandom.h"
-#include "TH2F.h"
-#include "TCanvas.h"
+#include <TControlBar.h>
+#include <TString.h>
+#include <TRandom.h>
+#include <TStyle.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TCanvas.h>
#include "AliRun.h"
#include "AliEMCALCalibData.h"
"Explains how to use EMCAL CDS menus");
menu->AddButton("Equal CC","SetCC(0)",
"Set equal calibration coefficients");
- menu->AddButton("Decalibrate","SetCC(1)",
+ menu->AddButton("Random De-calibration","SetCC(1)",
"Set random decalibration calibration coefficients");
+ menu->AddButton("Gaussian De-calibration","SetCC(2)",
+ "Set gausian decalibration calibration coefficients");
menu->AddButton("Read equal CC","GetCC(0)",
"Read initial equal calibration coefficients");
menu->AddButton("Read random CC","GetCC(1)",
"Read random decalibration calibration coefficients");
+ menu->AddButton("Read gaussian CC","GetCC(2)",
+ "Read gausian decalibration calibration coefficients");
menu->Show();
}
// Arguments:
// flag=0: all calibration coefficients are equal
// flag=1: all calibration coefficients random (decalibration)
+ // flag=2: all calibration coefficients have Gaussian random distribution (decalibration)
// Author: Boris Polishchuk (Boris.Polichtchouk at cern.ch)
TString DBFolder;
- Int_t firstRun = 0;
+ Int_t firstRun = 0; // What is this
Int_t lastRun = 10;
Int_t beamPeriod = 1;
char* objFormat = "";
lastRun = 12;
objFormat = "EMCAL random pedestals and ADC gain factors (12x48x24)";
}
+ else if (flag == 2) {
+ DBFolder ="local://DeCalibDB"; // create directory DeCalibDB in current directory
+ firstRun = 0;
+ lastRun = 12; // Why 12 ?
+ objFormat = "EMCAL random pedestals and gausian ADC gain factors (12x48x24)";
+ }
AliEMCALCalibData *calibda=new AliEMCALCalibData("EMCAL");
- Float_t fADCpedestal = 0.0;
+ Float_t fADCpedestal = 0.009;
Float_t fADCchannel = 0.0153; // 250 GeV / (16*1024)
+ // Float_t fADCchannel = 0.00305;
+ Float_t rDecalibration = 0.1 * fADCchannel; // 10% decalibration - just a guess
+ Float_t cc=0, ped;
TRandom rn;
Int_t nSMod = 12;
if(supermodule >= 10)
nRow = nRow2;
for(Int_t row=0; row< nRow; row++) {
+ cc = fADCchannel;
+ ped = fADCpedestal;
if (flag == 1) {
// Decalibration:
- // Spread calibration coefficients uniformly around central value
- fADCchannel = rn.Uniform(0.00140,0.00160);
- fADCpedestal = 0;
+ // Spread calibration coefficients uniformly with
+ // Cmax/Cmin = 5, (Cmax-Cmin)/2 = 0.0015
+ // and pedestals 0.005 +-10%
+ // fADCchannel = rn.Uniform(0.00075,0.00375);
+ cc = rn.Uniform(0.00140,0.00160);
+ ped = 0;
+ } else if (flag == 2) { // Gaussian
+ cc = rn.Gaus(fADCchannel, rDecalibration);
+ ped = rn.Uniform(0.0045,0.0055);
}
- calibda->SetADCchannel (supermodule,column,row,fADCchannel);
- calibda->SetADCpedestal(supermodule,column,row,fADCpedestal);
+ calibda->SetADCchannel (supermodule,column,row, cc);
+ calibda->SetADCpedestal(supermodule,column,row, ped);
cout<<"Set SM: "<<supermodule<<" col "<<column<<" row "<<row
- <<" chan "<<fADCchannel<<" ped fADCpedestal"<<endl;
+ <<" cc "<< cc <<" ped "<<ped<<endl;
}
}
}
AliCDBMetaData md;
md.SetComment(objFormat);
md.SetBeamPeriod(beamPeriod);
- md.SetResponsible("Boris Polichtchouk");
+ md.SetResponsible("Aleksei Pavlinov");
- AliCDBId id("EMCAL/Calib/Data",firstRun,lastRun);
+ AliCDBId id("EMCAL/Calib/Data",firstRun,lastRun); // create in EMCAL/Calib/Data DBFolder
AliCDBManager* man = AliCDBManager::Instance();
AliCDBStorage* loc = man->GetStorage(DBFolder.Data());
// Author: Yuri.Kharlov at cern.ch
TString DBFolder;
+ int drawKey=1;
if (flag == 0) {
DBFolder ="local://InitCalibDB";
}
- else if (flag == 1) {
- DBFolder ="local://DeCalibDB";
+ else if (flag == 1 || flag == 2) {
+ // DBFolder ="local://DeCalibDB"; // Get DB in current folder
+ DBFolder ="local:///data/r22b/ALICE/PROD/CALIBRATION_May_2007/PI0/DeCalibDB"; // absolute name - Aug 31, 2007
+ TString HOST(gSystem->Getenv("HOST"));
+ if(HOST.Contains("pc")) { // pdsf; May 31, 2007
+ DBFolder ="local:///eliza5/alice/pavlinov/PROD/CALIBRATION_May_2007/PI0/10GEV/DECALIB/DeCalibDB";
+ }
}
AliEMCALCalibData* clb = (AliEMCALCalibData*)
static const Int_t nCol = 48;
Int_t nRow = 24;
Int_t nRow2 = 12; //Modules 11 and 12 are half modules
+ Int_t nCC = 0;
TH2F *hPed[nSMod], *hGain[nSMod];
- TCanvas *cPed = new TCanvas("cPed" ,"Pedestals Mod 1-6" , 10,10,400,800);
- TCanvas *cGain = new TCanvas("cGain","Gain factors Mod 1-6", 410,10,400,800);
- TCanvas *cPed2 = new TCanvas("cPed2","Pedestals SMod 7-12", 10,10,400,800);
- TCanvas *cGain2 = new TCanvas("cGain2","Gain factors SMod 7-12", 410,10,400,800);
- cPed ->Divide(2,3);
- cGain ->Divide(2,3);
- cPed2 ->Divide(2,3);
- cGain2 ->Divide(2,3);
+ TH1F *hCCSum = new TH1F("hCCSum"," CC summary (in MeV) ", 200, 0.0, 20.);
+ TH1F *hPedSum = new TH1F("hPedSum"," pedestal summary (in MeV) ", 100, 4., 6.);
+
+ TCanvas *cPed=0, *cGain=0, *cPed2=0, *cGain2=0;
+ if(drawKey>1) {
+ cPed = new TCanvas("cPed" ,"Pedestals Mod 0-5" , 10,10,400,800);
+ cGain = new TCanvas("cGain","Gain factors Mod 0-5", 410,10,400,800);
+ cPed2 = new TCanvas("cPed2","Pedestals SMod 6-11", 10,10,400,800);
+ cGain2 = new TCanvas("cGain2","Gain factors SMod 6-11", 410,10,400,800);
+ cPed ->Divide(2,3);
+ cGain ->Divide(2,3);
+ cPed2 ->Divide(2,3);
+ cGain2 ->Divide(2,3);
+ }
+ TCanvas *cSum = new TCanvas("cSum" ,"summary" , 10,10,600,800);
+ cSum->Divide(1,2);
+
+ cout<<endl;
for (Int_t supermodule=0; supermodule<nSMod; supermodule++) {
if(supermodule >= 10)
Float_t ped = clb->GetADCpedestal(supermodule,column,row);
Float_t gain = clb->GetADCchannel (supermodule,column,row);
//cout<<"Get SM: "<<supermodule<<" col "<<column<<" row "<<row
- // <<" chan "<<gain<<endl;
- hPed[supermodule] ->SetBinContent(column,row,ped);
- hGain[supermodule]->SetBinContent(column,row,gain);
+ //<<" chan "<<gain<<endl;
+ hPed[supermodule] ->SetBinContent(column+1,row+1,ped*1.e+3); // in mev
+ hGain[supermodule]->SetBinContent(column+1,row+1,gain*1.e+3); // in mev
+
+ hPedSum->Fill(ped*1.e+3);
+ hCCSum->Fill(gain*1.e+3);
+
+ nCC++;
}
}
- if(supermodule < 7){
- cPed ->cd(supermodule);
- hPed[supermodule] ->Draw("lego2");
- cGain->cd(supermodule);
- hGain[supermodule]->Draw("lego2");
- }
- else{
- cPed2 ->cd(supermodule-6);
- hPed[supermodule] ->Draw("lego2");
- cGain2->cd(supermodule-6);
- hGain[supermodule]->Draw("lego2");
+ cout<<" Fill cc for SM "<< supermodule << " nCC "<< nCC << endl;
+
+ if(drawKey>1) {
+ if(supermodule < 6){
+ cPed ->cd(supermodule+1);
+ hPed[supermodule]->Draw("lego2");
+ cGain->cd(supermodule+1);
+ hGain[supermodule]->Draw("lego2");
+ }
+ else{
+ cPed2 ->cd(supermodule-5);
+ hPed[supermodule]->Draw("lego2");
+ cGain2->cd(supermodule-5);
+ hGain[supermodule]->Draw("lego2");
+ }
}
}
- cPed ->Print("pedestals_SM_1_6.eps");
- cGain ->Print("gains_SM_1-6.eps");
- cPed2 ->Print("pedestals_SM_7-12.eps");
- cGain2 ->Print("gains_SM_7-12.eps");
+ cout << " Get "<<nCC<<" calibration coeffs"<<endl;
+
+ cSum->cd(1);
+ gStyle->SetOptFit(111);
+ hCCSum->Fit("gaus");
+ hCCSum->SetLineWidth(2);
+ hCCSum->GetFunction("gaus")->SetLineColor(2);
+
+ cSum->cd(2);
+ hPedSum->Draw();
+ hPedSum->SetLineWidth(2);
+
+ cSum->Update();
+
+ /*
+ cPed ->Print("pedestals_SM_0_6.eps");
+ cGain ->Print("gains_SM_0_5.eps");
+ cPed2 ->Print("pedestals_SM_6_11.eps");
+ cGain2 ->Print("gains_SM_6_11.eps");
+ */
}
# Libraries will be linked against SHLIB
# ROOT libraries
-ROOTCLIBS := $(shell root-config --glibs) -lThread -lMinuit -lHtml -lVMC -lEG -lGeom -lTreePlayer -lXMLIO -lProof -lProofPlayer
+ROOTCLIBS := $(shell root-config --glibs) -lThread -lMinuit -lHtml -lVMC -lEG -lGeom -lTreePlayer -lXMLIO -lProof -lProofPlayer -lTable
ROOTPLIBS := -lEGPythia6
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff