added pi0 calibration, linearity, shower profile
authorpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 11 Sep 2007 19:38:15 +0000 (19:38 +0000)
committerpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 11 Sep 2007 19:38:15 +0000 (19:38 +0000)
14 files changed:
EMCAL/AliEMCALCalibCoefs.cxx
EMCAL/AliEMCALCalibCoefs.h
EMCAL/AliEMCALCell.cxx
EMCAL/AliEMCALCell.h
EMCAL/AliEMCALFolder.cxx
EMCAL/AliEMCALFolder.h
EMCAL/AliEMCALHistoUtilities.cxx
EMCAL/AliEMCALHistoUtilities.h
EMCAL/AliEMCALPi0SelectionParam.cxx
EMCAL/AliEMCALPi0SelectionParam.h
EMCAL/AliEMCALRecPointsQaESDSelector.cxx
EMCAL/AliEMCALRecPointsQaESDSelector.h
EMCAL/AliEMCALSuperModule.cxx
EMCAL/AliEMCALSuperModule.h

index 7d171e3..9158222 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Log$ */
+/* $Log$co: warning: `/* $Log' is obsolescent; use ` * $Log'.
+
+ * Revision 1.1  2007/08/08 15:58:01  hristov
+ * */
+
 
 //_________________________________________________________________________
 //    Calibration coefficients
 #include "AliEMCALGeometry.h"
 
 #include "TH1F.h"
+#include <TString.h>
 
-TableClassImpl(AliEMCALCalibCoefs,calibCoef)
+ClassImp(calibCoef) 
 
-// Get initial Calib Data from DB
-AliEMCALCalibCoefs* AliEMCALCalibCoefs::GetCalibTableFromDb(const char *tn)
+calibCoef::calibCoef() : absId(-1), cc(-1), eCc(-1)
+{ }
+calibCoef::calibCoef(const Int_t id, const Double_t c, const Double_t ec) :
+absId(id), cc(c), eCc(ec)
+{}
+// ----------------------------------------------------------------------
+
+ClassImp(AliEMCALCalibCoefs)
+
+AliEMCALCalibCoefs::AliEMCALCalibCoefs() : TNamed("",""), fTable(0), fCurrentInd(0), fCalibMethod(0)
+{
+}
+
+AliEMCALCalibCoefs::AliEMCALCalibCoefs(const char* name, const Int_t nrow) : TNamed(name,"table of cell information") , fTable(0), fCurrentInd(0), fCalibMethod(0)
+{
+  fTable = new TObjArray(nrow);
+}
+
+void AliEMCALCalibCoefs::AddAt(calibCoef* r)
 {
-  // Initial cc with decalibration
-  char* dbString  = "local:///data/r22b/ALICE/PROD/CALIBRATION_May_2007/PI0/PDSF/10GEV/DECALIB/DeCalibDB"; 
+  (*fTable)[fCurrentInd] = new calibCoef(*r);
+  fCurrentInd++;
+}
+
+AliEMCALCalibCoefs::~AliEMCALCalibCoefs()
+{
+  if(fTable) {
+    fTable->Delete();
+    delete fTable;
+  }
+}
+
+calibCoef* AliEMCALCalibCoefs::GetTable(Int_t i) const
+{
+  return (calibCoef*)fTable->At(i);
+}
+
+
+// Get initial Calib Data from DB
+AliEMCALCalibCoefs* AliEMCALCalibCoefs::GetCalibTableFromDb(const char *tn, AliEMCALCalibData **calData)
+{ 
+  //
+  // See ~/macros/ALICE/sim.C  for choice of CDB
+  // Get calib. table which was used than calculated rec.points
+  // 
+  static char *calibType = "EMCAL/Calib/*";
+  static char *calibTypeData = "EMCAL/Calib/Data";
+  // Initial cc
+  calData[0] = 0;
+
+  AliCDBManager* man = AliCDBManager::Instance(); 
+  AliCDBStorage* specificStorage = man->GetSpecificStorage(calibType);
+
   AliEMCALCalibData* caldata = (AliEMCALCalibData*)
-    (AliCDBManager::Instance()->GetStorage(dbString)
-     ->Get("EMCAL/Calib/Data",gAlice->GetRunNumber())->GetObject());
+  specificStorage->Get(calibTypeData, gAlice->GetRunNumber())->GetObject();
   if(caldata == 0) return 0;
 
-  AliEMCALCalibCoefs *tab = new  AliEMCALCalibCoefs(tn);
+  AliEMCALGeometry *g = AliEMCALGeometry::GetInstance();
+
+  AliEMCALCalibCoefs *tab = new  AliEMCALCalibCoefs(tn,g->GetNCells());
   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);
@@ -56,13 +109,17 @@ AliEMCALCalibCoefs* AliEMCALCalibCoefs::GetCalibTableFromDb(const char *tn)
 
     calibCoef r;
     r.absId = id;
-    r.cc    = caldata->GetADCchannel(nSupMod, ietaCell, iphiCell);
+    r.cc    = caldata->GetADCchannel(nSupMod, ietaCell, iphiCell); // column(ietaCell) : row(iphiCell)
     r.eCc   = 0.0;
 
     tab->AddAt(&r);
   }
-  tab->Purge();
+  // tab->Purge();
   tab->SetCalibMethod(AliEMCALCalibCoefs::kPI0);
+  calData[0] = caldata; 
+
+  printf("\n <I> AliEMCALCalibCoefs::GetCalibTableFromDb \n name |%s| title |%s| | storage |%s|\n", 
+        caldata->GetName(), caldata->GetTitle(),  specificStorage->GetURI().Data());
 
   return tab;
 }
@@ -70,26 +127,78 @@ AliEMCALCalibCoefs* AliEMCALCalibCoefs::GetCalibTableFromDb(const char *tn)
 TH1F* AliEMCALCalibCoefs::GetHistOfCalibTableFromDb(const char *tn)
 {
   // First SM only
-  AliEMCALGeometry *g = AliEMCALGeometry::GetInstance("");
+  AliEMCALCalibData *calData[1]; // unused here
 
-  AliEMCALCalibCoefs* tab = GetCalibTableFromDb(tn);
+  AliEMCALCalibCoefs* tab = GetCalibTableFromDb(tn, calData);
   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++){
+  for(Int_t i=0; i<tab->GetSize(); i++){
     r = tab->GetTable(i);
     if(i>=1152) break;
-    h->Fill(r->cc*1000.); 
+    h->Fill(r->cc*1000.); // GEV ->MeV
   }
   delete tab;
   return h;
 }
 
+AliEMCALCalibData* AliEMCALCalibCoefs::GetCalibTableForDb(const AliEMCALCalibCoefs *tab, const char* dbLocation,
+const char* coment)
+{
+  printf("<I> AliEMCALCalibCoefs::GetCalibTableForDb started \n");
+  // See EMCAL/macros/CalibrationDB/AliEMCALSetCDB.C
+  // Define and save to CDB 
+  AliEMCALCalibData* caldata=0;
+  if(tab==0) return caldata;
+
+  Int_t firstRun   =  0;
+  Int_t lastRun    = 10;
+  Int_t beamPeriod =  1;
+  char* objFormat  = "";
+  caldata = new AliEMCALCalibData("EMCAL");
+  caldata->SetTitle(coment);
+
+  AliEMCALGeometry *g = AliEMCALGeometry::GetInstance();
+
+  for(int id=0; id<tab->GetSize(); id++){
+    Float_t ped=0.;
+
+    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 = tab->GetTable(id);
+    // ietaCell - column; iphiCell - row
+    caldata->SetADCchannel (nSupMod, ietaCell, iphiCell, r->cc);
+    caldata->SetADCpedestal(nSupMod, ietaCell, iphiCell, ped);
+  }
+  printf("<I> Fill AliEMCALCalibData table \n");
+  //Store calibration data into database
+  
+  AliCDBMetaData md;
+  md.SetComment(objFormat);
+  md.SetBeamPeriod(beamPeriod);
+  md.SetResponsible("Aleksei Pavlinov");
+
+  AliCDBManager* man = AliCDBManager::Instance();
+  if(man == 0) {
+    printf("<E> AliEMCALCalibCoefs::GetCalibTableForDb : define AliCDBManager, NO saving  !! \n"); 
+  } else {
+    printf("<I> AliCDBManager %p \n", man); 
+    AliCDBId id("EMCAL/Calib/Data",firstRun,lastRun); // create in EMCAL/Calib/Data DBFolder 
+    TString DBFolder(dbLocation);
+    AliCDBStorage* loc = man->GetStorage(DBFolder.Data());
+    loc->Put(caldata, id, &md);
+  }
+
+  return caldata;
+}
+
 calibCoef *AliEMCALCalibCoefs::GetRow(const int absId)
 {
   calibCoef *r=0;
-  for(int id=0; id<GetNRows(); id++){
+  for(int id=0; id<fTable->GetSize(); id++){
     r = GetTable(id);
     if(r->absId == absId) return r;
   }
@@ -98,13 +207,13 @@ calibCoef *AliEMCALCalibCoefs::GetRow(const int absId)
 
 void AliEMCALCalibCoefs::PrintTable()
 {
-  printf(" Table : %s : nrows %i \n", GetName(), int(GetNRows()));
-  for(int i=0; i<GetNRows(); i++) PrintTable(i);
+  printf(" Table : %s : nrows %i \n", GetName(), int(fTable->GetSize()));
+  for(int i=0; i<fTable->GetSize(); i++) PrintTable(i);
 }
 
 void AliEMCALCalibCoefs::PrintTable(const Int_t i)
 {
-  if(i>=GetNRows()) return;
+  if(i>=fTable->GetSize()) return;
   printf("row %i \n", i);
   PrintRec(GetTable(i));
 }
index 9b420b0..6bda4e8 100644 (file)
@@ -6,42 +6,76 @@
 /* $Id$ */
 
 //_________________________________________________________________________
-//    Table of Calibration coefficients  
+//  Table of Calibration coefficients  
+//  Should be extended.
 //                  
 //*-- Author: Aleksei Pavlinov (WSU, Detroit, USA) 
 
 // --- ROOT system ---
-#include <TTable.h>
+#include <TNamed.h>
+#include <TObjArray.h>
 
 // unit is GeV
-struct  calibCoef {
+class calibCoef : public TObject {
+ public:
+  virtual const char* GetName() const {return Form("CC%5.5i",absId);}
+  calibCoef();
+  calibCoef(const Int_t id, const Double_t c, const Double_t ec);
+  virtual ~calibCoef() {};
+
   Int_t    absId; // absolute id of cell 
   Double_t cc;    // Calib. coef
   Double_t eCc;   // Calib. coef. error
+  // 
+  ClassDef(calibCoef,1) // Cell calibration information 
 };
 
 class TH1F;
+class AliEMCALCalibData;
 
-class AliEMCALCalibCoefs : public TTable {
+class AliEMCALCalibCoefs : public TNamed {
  public:
   enum EEmcalCalibType {kMC, kEQUALIZATION, kMIP, kPI0}; // type of EMCAL calibrations 
 
+  AliEMCALCalibCoefs(); // default constractor
+  AliEMCALCalibCoefs(const char* name, const Int_t nrow);
+  virtual ~AliEMCALCalibCoefs();
+
+  AliEMCALCalibCoefs & operator = (const AliEMCALCalibCoefs  & /*rvalue*/) {
+    // assignement operator requested by coding convention but not needed
+    Fatal("operator =", "not implemented");
+    return *this;
+  };
+  void        AddAt(calibCoef* r);
+  calibCoef*  GetTable(Int_t i) const;
+  Int_t       GetSize()  const {return fTable->GetSize();}
+  Int_t       GetNRows() const {return fCurrentInd;}
+  void        Purge() {/* nothing */};
+
   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");
+  static AliEMCALCalibCoefs *GetCalibTableFromDb(const char *tn="CCIN", AliEMCALCalibData **calData=0);
+  //  const char* dbLocation="local:///data/r22b/ALICE/PROD/CALIBRATION_May_2007/PI0/PDSF/10GEV/DECALIB/DeCalibDB");
+  static TH1F *GetHistOfCalibTableFromDb(const char *tn="CCIN");
+  //const char* dbLocation="local:///data/r22b/ALICE/PROD/CALIBRATION_May_2007/PI0/PDSF/10GEV/DECALIB/DeCalibDB");
+  // Create DB calib table
+  static AliEMCALCalibData*  GetCalibTableForDb(const AliEMCALCalibCoefs *tab, 
+  const char* dbLocation="local://.", // current directory
+  const char* coment="pi0 calibration, MC, Jun 2007");
   // Menu
   void PrintTable();                 // *MENU*
   void PrintTable(const Int_t i);    // *MENU*
   void PrintRec(calibCoef *r);
 
+  //
  protected:
+  TObjArray *fTable;
+  Int_t fCurrentInd;
   Int_t fCalibMethod;  // method of calibration - EEmcalCalibType
 
-  ClassDefTable(AliEMCALCalibCoefs , calibCoef)
-  ClassDef(AliEMCALCalibCoefs,1) // Table of Calibration coefficients  
+  ClassDef(AliEMCALCalibCoefs,2) // Table of Calibration coefficients  
 };
 
 #endif // ALIEMCALCalibCoefs_H
index ef5cc37..3faef11 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Id:*/
+/*
+$Log$
+*/ 
 
 //_________________________________________________________________________
-// Top EMCAL folder which will keep all information about EMCAL itself,
-// super Modules (SM), modules, towers, set of hists and so on.
+// Cell folder which will keep all information about cell(tower) itself
+//  Initial version was created with TDataSet staf
+//  TObjectSet -> TFolder; Sep 6, 2007
 //
 //*-- Author: Aleksei Pavlinov (WSU, Detroit, USA) 
 
@@ -27,6 +30,7 @@
 #include "AliEMCALFolder.h"
 #include "AliEMCALSuperModule.h"
 #include "AliEMCALCalibData.h"
+#include "AliEMCALRecPointsQaESDSelector.h"
 
 #include "AliEMCALCalibCoefs.h"
 
@@ -36,7 +40,6 @@
 #include <TH1.h>
 #include <TF1.h>
 #include <TNtuple.h>
-#include <TObjectSet.h>
 
 typedef  AliEMCALHistoUtilities u;
 
@@ -47,18 +50,19 @@ Double_t MPI0         = 0.13498; // mass of pi0
 Double_t MPI02        = MPI0*MPI0; // mass**2
 
 AliEMCALCell::AliEMCALCell() : 
-TObjectSet(), 
+TFolder(), 
+ fParent(0),fLh(0),
 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)), 
+  TFolder(Form("Cell%4.4i",absId),title), 
+ fParent(0),fLh(0),
 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);
@@ -70,6 +74,7 @@ AliEMCALCell::~AliEMCALCell()
 {
   // dtor
 }
+//-------------------------------------------------------------------------------------
 
 void AliEMCALCell::SetCCfromDB(AliEMCALCalibData *ccDb)
 {
@@ -84,8 +89,14 @@ void AliEMCALCell::SetCCfromDB(AliEMCALCalibData *ccDb)
 
 void AliEMCALCell::SetCCfromCCTable(AliEMCALCalibCoefs *t)
 {
-  if(t == 0) return;
-  this->AddObject((TObject*)BookHists(), kTRUE);
+  if(t == 0) {
+    //    Dump();
+    return;
+  }
+  if(fLh == 0) {
+    fLh = BookHists();
+    Add(fLh);
+  }
 
   calibCoef *r = t->GetTable(fAbsId);
   if(r && r->absId == fAbsId) {
@@ -115,11 +126,11 @@ 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);
+  printf("<I> AliEMCALCell::FitHist : |%s| is started : opt %s\n", h->GetName(), opt);
   TString tit(h->GetTitle());
 
   TF1 *GausPol2 = 0, *g=0, *bg=0;
-  if(h->GetListOfFunctions()->GetSize() == 0) {
+  if(h->GetListOfFunctions()->GetSize() == 0 || 1) {
     g = u::Gausi(name, 0.0, 0.4, h); // gaus estimation
 
     g->SetParLimits(0, h->Integral()/20., h->Integral());
@@ -154,6 +165,7 @@ void AliEMCALCell::FitHist(TH1* h, const char* name, const char* opt)
   } else {
     GausPol2 = (TF1*)h->GetListOfFunctions()->At(0);
     optFit = "IME+";
+    printf("<I> Function is defined alredy : %s optFit %s \n", GausPol2->GetTitle(), optFit.Data());
   }
   //  optFit = "IME+";
   h->Fit(GausPol2, optFit.Data(),"", 0.01, 0.28);
@@ -162,12 +174,12 @@ void AliEMCALCell::FitHist(TH1* h, const char* name, const char* opt)
     gStyle->SetOptFit(111);
     u::DrawHist(h,2);
   }
-  printf("<I> AliEMCALCell::FitHist : h %p |%s| is ended \n\n", h, h->GetName());
+  printf("<I> AliEMCALCell::FitHist : |%s| is ended \n\n", h->GetName());
 }
 
 void AliEMCALCell::FitEffMassHist(const char* opt)
 {
-  AliEMCALFolder* EMCAL = (AliEMCALFolder*)(GetParent()->GetParent()->GetParent()); 
+  AliEMCALFolder* EMCAL = AliEMCALRecPointsQaESDSelector::GetEmcalFolder();
   Int_t it = EMCAL->GetIterationNumber();
 
   TH1* h = (TH1*)GetHists()->At(0);
@@ -178,7 +190,7 @@ void AliEMCALCell::FitEffMassHist(const char* opt)
   if(fFun) {
     Double_t mpi = fFun->GetParameter(1), mpi2 = mpi*mpi;
     Double_t ccTmp = fCcIn * MPI02 / mpi2;
-    if(it<=6) {
+    if(it<=1) { // Jul 16, 2007
       fCcOut = ccTmp;
     } else {
       fCcOut = (ccTmp + fCcIn)/2.;
@@ -188,14 +200,15 @@ void AliEMCALCell::FitEffMassHist(const char* opt)
   printf(" %s | fCcIn %6.5f -> % 6.5f <- fCcOut \n", GetTitle(), fCcIn , fCcOut);
 }
 
-void AliEMCALCell::Print()
+void AliEMCALCell::PrintInfo()
 {
+  printf(" %s %s \n", GetName(), GetTitle());
+  if(fLh == 0 ) return;
   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();
+  // if(f) f->Dump();
 }
 
 TList* AliEMCALCell::BookHists()
@@ -203,7 +216,7 @@ TList* AliEMCALCell::BookHists()
   gROOT->cd();
   TH1::AddDirectory(1);
 
-  AliEMCALFolder* EMCAL = (AliEMCALFolder*)(GetParent()->GetParent()->GetParent()); 
+  AliEMCALFolder* EMCAL = AliEMCALRecPointsQaESDSelector::GetEmcalFolder();
   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);
index a37efd9..8ead9fc 100644 (file)
@@ -7,10 +7,12 @@
 
 //_________________________________________________________________________
 //  EMCAL cell - keep everyrhing for calibration task     
+//  Initial version was created with TDataSet staf
+//  TObjectSet -> TFolder; Sep 6, 2007
 //                  
 //*-- Author: Aleksei Pavlinov (WSU, Detroit, USA) 
 
-#include <TObjectSet.h>
+#include <TFolder.h>
 
 class TList;
 class TH1;
@@ -20,7 +22,7 @@ class TNtuple;
 class AliEMCALCalibData;
 class AliEMCALCalibCoefs;
 
-class AliEMCALCell : public TObjectSet {
+class AliEMCALCell : public TFolder {
 
  public:
   
@@ -32,7 +34,9 @@ class AliEMCALCell : public TObjectSet {
   void SetCCfromDB(AliEMCALCalibData *ccDb);  // obsolete
   void SetCCfromCCTable(AliEMCALCalibCoefs *t);
 
-  TList*   GetHists() {return (TList*)fObj;}
+  TList*   GetHists() {return fLh;}
+  TObject* GetParent() {return fParent;}
+  void     SetParent(TObject *parent) {fParent=parent;}
   Int_t    GetAbsId()  const {return fAbsId;}
   Int_t    GetSupMod() const {return fSupMod;}
   Int_t    GetModule() const {return fModule;}
@@ -47,8 +51,11 @@ class AliEMCALCell : public TObjectSet {
   static void FitHist(TH1* h, const char* name="",const char* opt="");
   // Menu
   void FitEffMassHist(const char* opt=""); //*MENU*
-  void Print();                            //*MENU*
+  void PrintInfo();                            //*MENU*
  protected:
+  TObject* fParent; // parent
+  TList* fLh;
+  //
   Int_t fAbsId;   // abs cell id 
   Int_t fSupMod;  // super module number
   Int_t fModule;  // module number inside SM
@@ -58,7 +65,7 @@ class AliEMCALCell : public TObjectSet {
   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)
+  Double_t fCcOut; // output cc  in GeV (from fit now)
 
   TF1*   fFun;     //! fitting function - gaus + pol2
   //
index 91f1421..b0f95f1 100644 (file)
 #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>
+#include <TLegendEntry.h>
+#include <TLine.h>
 
   const TString AliEMCALFolder::fgkBaseFolderName("EMCAL");
   const TString AliEMCALFolder::fgkCCFirstName("CCFIRST");
   const TString AliEMCALFolder::fgkCCinName("CCIN");
   const TString AliEMCALFolder::fgkCCoutName("CCOUT");
+  const TString AliEMCALFolder::fgkDirOfRootFiles("$HOME/ALICE/SHISHKEBAB/RF/CALIB/JUL16/");
 
 typedef  AliEMCALHistoUtilities u;
 
@@ -73,27 +74,25 @@ ClassImp(AliEMCALFolder)
   TList *lObj = 0;
 
 AliEMCALFolder::AliEMCALFolder() : 
-  TObjectSet(), 
+  TFolder(), 
   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),
+  TFolder(name,title),
   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)),
+  TFolder(Form("%s_%2.2i", AliEMCALFolder::fgkBaseFolderName.Data(),it),title),
   fCounter(it), fGeometry(0), fNumOfCell(0), fLhists(0), fLofCells(0),fPi0SelPar(0),fCalibData(0),
   fCellNtuple(0)
 {
-  SetTitle(title);
   Init(putToBrowser);
 }
 
@@ -106,22 +105,25 @@ void AliEMCALFolder::Init(Bool_t putToBrowser)
 {
   lObj = new TList;
   lObj->SetName("Objects"); // name is good ?
-  this->AddObject((TObject*)lObj, kTRUE);
+  this->Add((TObject*)lObj);
+  //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());
+  // Jun 13, 2007;
+  // Jul 13 - See ~/macros/ALICE/sim.C  for choice of CDB
+  AliEMCALCalibData *calData[1];
+  // Get from defined arrea;
+  // First table is table which used in rec.points finder.
+  Add(AliEMCALCalibCoefs::GetCalibTableFromDb(fgkCCFirstName.Data(),calData));
+  fCalibData = calData[0];
   lObj->Add(fCalibData);
-  // Jun 13, 2007
-  Add(AliEMCALCalibCoefs::GetCalibTableFromDb(fgkCCFirstName.Data()));
-  if(GetIterationNumber()==1) {
-    Add(AliEMCALCalibCoefs::GetCalibTableFromDb()); // first iteration only
+  if(GetIterationNumber()<=1) {
+    Add(AliEMCALCalibCoefs::GetCalibTableFromDb(fgkCCinName.Data(), calData));
   }
+
   // Selection Parameter
   fPi0SelPar = AliEMCALPi0SelectionParam::Set1();
   this->Add(fPi0SelPar);
@@ -141,10 +143,9 @@ void AliEMCALFolder::Init(Bool_t putToBrowser)
 
 AliEMCALSuperModule* AliEMCALFolder::GetSuperModule(const Int_t nm)
 {
-  TDataSet *set = 0; 
   AliEMCALSuperModule *sm = 0;
 
-  set = FindByName(Form("SM%2.2i",nm));
+  TObject *set = FindObject(Form("SM%2.2i",nm));
   if(set) sm = (AliEMCALSuperModule*)set;
 
   return sm;
@@ -191,16 +192,22 @@ void AliEMCALFolder::FillPi0Candidate(const Double_t mgg, AliESDCaloCluster* 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();
+  const Short_t* absId   = cl->GetDigitIndex()->GetArray();
+
+  AliEMCALCalibCoefs *tFirst = GetCCFirst();
+  calibCoef *rFirst=0;
 
-  int    indMax = 0;
-  double emax   = cl->GetDigitAmplitude()->At(0);
+  int    indMax = 0, id=0;
+  id = Int_t(absId[0]);
+  rFirst = tFirst->GetTable(id);
+  double emax   = cl->GetTrueDigitEnergy(indMax, rFirst->cc);
   if(nDigits > 1) {
     for(int i=1; i<nDigits; i++) {
-      if(emax < cl->GetDigitAmplitude()->At(i)) {
+      id = Int_t(absId[i]);
+      rFirst = tFirst->GetTable(id);
+      if(emax < cl->GetTrueDigitEnergy(i, rFirst->cc)) {
         indMax = i;
-        emax   = cl->GetDigitAmplitude()->At(i);
+        emax   = cl->GetTrueDigitEnergy(i, rFirst->cc);
       }
     }
   }
@@ -218,7 +225,7 @@ void AliEMCALFolder::FillPi0Candidate(const Double_t mgg, Int_t absIdMax, Int_t
   //
   static Int_t nSupModMax, nModuleMax, nIphiMax, nIetaMax, iphiCellMax, ietaCellMax;
   static AliEMCALCell* cell;
-  static TDataSet *set; 
+  static TFolder *set; 
   static AliEMCALSuperModule *sm;
 
   fGeometry->GetCellIndex(absIdMax, nSupModMax, nModuleMax, nIphiMax, nIetaMax);
@@ -233,20 +240,24 @@ void AliEMCALFolder::FillPi0Candidate(const Double_t mgg, Int_t absIdMax, Int_t
     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));
+    set  = dynamic_cast<TFolder*>(FindObject(Form("SM%2.2i",nSupModMax)));
     if(set==0) {
       sm = new  AliEMCALSuperModule(nSupModMax);
       Add(sm);
+      sm->SetParent(this);
       sm->Init();
     } else {
-      sm = (AliEMCALSuperModule*)set;
+      sm = dynamic_cast<AliEMCALSuperModule*>(set);
+    }
+    if(sm) {
+      sm->AddCellToEtaRow(cell, ietaCellMax);
+      cell->SetParent(sm);
     }
-    sm->AddCellToEtaRow(cell, ietaCellMax);
     // 
     cell->SetCCfromCCTable(GetCCIn());
   } else {
     cell = GetCell(absIdMax);
-    set  = FindByName(Form("SM%2.2i",nm));
+    set  = dynamic_cast<TFolder*>(FindObject(Form("SM%2.2i",nm)));
     if(set) sm = (AliEMCALSuperModule*)set;
   }
   if(sm == 0) assert(0);
@@ -285,9 +296,9 @@ void AliEMCALFolder::FitAllSMs()
 
 AliEMCALCalibCoefs* AliEMCALFolder::GetCCTable(const char* name)
 {
-  TDataSet *ds = FindByName(name);
-  if(ds) return (AliEMCALCalibCoefs*)ds;
-  else   return 0;
+  TObject *obj = FindObject(name);
+  if(obj) return (AliEMCALCalibCoefs*)obj;
+  else    return 0;
 }
 
 Int_t AliEMCALFolder::GetSMNumber(AliESDCaloCluster* cl)
@@ -303,55 +314,75 @@ Int_t AliEMCALFolder::GetSMNumber(AliESDCaloCluster* cl)
 
 // Recalibration staf - Jun 18,2007
 AliEMCALRecPoint* AliEMCALFolder::GetRecPoint(AliESDCaloCluster *cl, AliEMCALCalibCoefs *tOld,AliEMCALCalibCoefs *tNew,
-TList *l)
+TList *l, Double_t deff, Double_t w0, Double_t phiSlope)
 {
   //
   // Static function;
   // Get recalibrated rec.point from ESD cluster
+  // If tNew == 0 -> get ideal calibration with  ADCCHANNELEC
   //
-  // 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; 
+  //static Float_t ECAW0 = 4.5; // hard coded now - see AliEMCALClusterizerv1::InitParameters()
+  static Double_t ECAW0 = 5.5; // Beter case for simulation 
+  Int_t ampDigi=0, indMax=-1;
+  Double_t eDigiNew=0.0, eDigiMax=0.0; 
+  static TArrayD ed;
 
+  if(w0 > 0.5) ECAW0 = w0;
   AliEMCALRecPoint *rp=0;
   calibCoef *rOld=0, *rNew=0;
-  if(cl && tOld && tNew){
+  //   printf(" AliEMCALFolder::GetRecPoint() : RECALIBRATION : w0 %f ECAW0 %f \n", w0, ECAW0);
+  if(cl && tOld){
     // cl->PrintClusterInfo(1);
     const Int_t ndg   = cl->GetNumberOfDigits();
-    //const Int_t nprim = cl->GetNumberOfPrimaries();
+    //    const Int_t nprim = cl->GetNumberOfPrimaries();
 
-    const Short_t* prim    = cl->GetLabels()->GetArray();
+    const Short_t  prim    = cl->GetLabel();
     const Short_t* dgAbsId = cl->GetDigitIndex()->GetArray();
-    //    const UShort_t* dgAmp   = cl->GetDigitAmplitude()(); // This is energy - bad definition
+    //    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;
+    ed.Set(ndg); // resize array
     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);
+      absId   = Int_t(dgAbsId[i]);
+      rOld    = tOld->GetTable(absId);
+      ampDigi = cl->GetTrueDigitAmplitude(i, rOld->cc); // True amplitude
 
-      new(digits[i]) AliEMCALDigit(Int_t(prim[i]),0,absId, ampDigi, 0.0, i, 0.0); 
+      new(digits[i]) AliEMCALDigit(Int_t(prim),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
+
+      if(tNew) rNew = tNew->GetTable(absId);
+      if(rNew) {
+        rNew  = tNew->GetTable(absId);
+        eDigiNew = Double_t(ampDigi) * rNew->cc;     // Recalibrate with new cc
+      } else {
+        eDigiNew = Double_t(ampDigi) * ADCCHANNELEC; // Ideal calibration
+      }
+      //eDigiNew = Double_t(cl->GetTrueDigitEnergy(i)); // Copy from ESD for checking
       rp->AddDigit(*dg, Float_t(eDigiNew));
+      ed[i] = eDigiNew;
+      if(eDigiMax<eDigiNew) {
+        eDigiMax = eDigiNew;
+        indMax   = i;
+      }
       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
+    if(indMax>=0) rp->SetIndMaxDigit(indMax);
+    if(deff > 0.0) { // for fit
+      rp->EvalLocalPositionFit(deff, ECAW0, phiSlope, &digits); // I need just position
+    } else { // get w0 and deff from parametrisation - Sep 4, 2007 
+      rp->EvalLocalPosition2(&digits, ed);
+    }
     digits.Delete();
   }
   //rp->Print("print");
@@ -362,41 +393,54 @@ void  AliEMCALFolder::Save(const char *fn, const char *opt)
 { 
   //
   // Jun 5, 2007; See TFileIter and StFMC.cxx
-  //
-  TString FN(fn);
+  // Jul 16 - added fgkDirOfRootFiles
+  // Sep 7, 2007 - should be changed without TFileIter
+  /*
+  TString FN = fgkDirOfRootFiles;
+  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
+  UInt_t runNumber = 0; // 0 now, - may statistics on selector
   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
+  // Jul 16 - added fgkDirOfRootFiles
   //
+  printf("<I> AliEMCALFolder::Read(%s,%s)\n",fn,opt); 
   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);
+  TString FN = fgkDirOfRootFiles;
+  FN += fn;
+  if(FN.Contains(".root")==0) FN += ".root";
 
-  TKey *key = (TKey*)l->At(0);
-  EMCAL = (AliEMCALFolder*)key->ReadObj();
-  f.Close();
-  if(EMCAL) EMCAL->InitAfterRead();
+  TFile f(FN.Data(),opt);
+  if(f.IsOpen()) {
+    TList *l = f.GetListOfKeys();
+    printf("<I> The total number of the objects: %d \n File %s\n", l->GetSize(), FN.Data());
 
+    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");
+{ 
+  lObj    = (TList*)FindObject("Objects");
+  if(lObj) {
+    fLhists = (TList*)lObj->FindObject("HistsOfEmcal");
+  }
 }
 
 void AliEMCALFolder::DrawQA(const int nsm)
@@ -477,11 +521,12 @@ void AliEMCALFolder::CreateCellNtuple()
   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));
+    TFolder *setEta = dynamic_cast<TFolder*>(sm->FindObject(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);
+      TList* l = (TList*)setEta->GetListOfFolders();
+      for(int phi=0; phi<l->GetSize(); phi++) { // cycle on cells (phi directions)
+        cell = (AliEMCALCell*)l->At(phi);
         if(cell) {
           cell->FillCellNtuple(fCellNtuple);
           //printf(" fill cell %s : %s \n", cell->GetName(), cell->GetTitle());
@@ -491,15 +536,15 @@ void AliEMCALFolder::CreateCellNtuple()
   }
   fCellNtuple->Print();
   lObj->Add(fCellNtuple); 
-  // --
-  CreateAndFillAdditionalHists();
 }
 
 void AliEMCALFolder::CreateAndFillAdditionalHists()
 {
+  CreateCellNtuple();
   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.));
+  fLhists->Add(new TH1F("04_CCoutOnEdge2", "cc out on edge of calorimeter(2) (in MeV)", 70, 12., 19.));
   // Fill
   Float_t* args;
   for(Int_t i=0; i<(Int_t)fCellNtuple->GetEntries(); i++){
@@ -508,8 +553,9 @@ void AliEMCALFolder::CreateAndFillAdditionalHists()
     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);
+    if     ((phi==0||phi==23) || (eta==0||eta==47)) u::FillH1(fLhists, 2, cc);
+    else if((phi==1||phi==22) || (eta==1||eta==46)) u::FillH1(fLhists, 4, cc); // next to edge
+    else                                            u::FillH1(fLhists, 3, cc);
   }
   // Drawing
   Int_t wh=530, ww=750;
@@ -525,9 +571,9 @@ void AliEMCALFolder::CreateAndFillAdditionalHists()
   h1->SetTitle("CC distribution after #pi^{0} calibration");
   h1->SetXTitle("  MeV  ");
   h1->SetYTitle("  N  ");
-  TLatex *lat1 = u::lat(Form("rel.width = %4.2f%%", 
+  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)", 
+  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) {
@@ -539,9 +585,35 @@ void AliEMCALFolder::CreateAndFillAdditionalHists()
   u::DrawHist(hccFirst,2,1,"same",3);
 
   
+  // Ideal calibration - Jul 18, 2007
+  Double_t ADCCHANNELEC = 0.0153, ccIdeal =  ADCCHANNELEC*1.e+3;
+  Double_t ym = h1->GetMaximum();
+  TLine *l = new TLine(ccIdeal,-ym*0.05, ccIdeal,ym*1.05);
+  l->SetLineColor(kBlue);
+  l->Draw(); 
+
   TLegend *leg = new TLegend(0.1,0.6, 0.45,0.85);
   leg->AddEntry(hccFirst, "Initial cc ", "L");
   leg->AddEntry(h1, "Final cc", "L");
+  TLegendEntry *le=0;
+  if(l) {
+    le  = leg->AddEntry(l, "Ideal calibration", "L");
+    le->SetTextColor(l->GetLineColor());
+  }
+
+  TH1 *hCcEdge  = (TH1*)fLhists->At(2);
+  TH1 *hCcEdge2 = (TH1*)fLhists->At(4);
+  if(1) {
+    u::DrawHist(hCcEdge,2,kGreen,"same",1);
+    le = leg->AddEntry(hCcEdge , "Edge cell", "L");
+    le->SetTextColor(hCcEdge->GetLineColor());
+
+    u::DrawHist(hCcEdge2,2, 28,"same",1); // 28 - like brown
+    le = leg->AddEntry(hCcEdge2 , "Edge cell (2)", "L");
+    le->SetTextColor(hCcEdge2->GetLineColor());
+    u::DrawHist(h1,2,1,"same");
+  }
+
   leg->Draw();
 
   c->Update();
@@ -553,6 +625,7 @@ void AliEMCALFolder::TestSMStruct()
   for(int m=0; m<12; m++) {
     AliEMCALSuperModule *sm = new AliEMCALSuperModule(m);
     Add(sm);
+    sm->SetParent(this);
   }
 }
 
index 523e209..7fae4bf 100644 (file)
@@ -3,16 +3,18 @@
 /* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice     */
 
-/* $log:$ */
+/* $Log$ */
 
 //_________________________________________________________________________
-//  Top EMCAL folder - keep everyrhing for calibration task     
+//  Top EMCAL folder - keep everyrhing for calibration task
+//  Initial version was created with TDataSet staf
+//  TObjectSet -> TFolder; Sep 5, 2007
 //                  
 //*-- Author: Aleksei Pavlinov (WSU, Detroit, USA) 
 
 // --- ROOT system ---
 
-#include <TObjectSet.h>
+#include <TFolder.h>
 #include <TString.h>
 
 class AliEMCALGeometry;
@@ -28,7 +30,7 @@ class AliEMCALRecPoint;
 class TList;
 class TNtuple;
 
-class AliEMCALFolder : public TObjectSet {
+class AliEMCALFolder : public TFolder {
 
  public:
   
@@ -53,14 +55,15 @@ class AliEMCALFolder : public TObjectSet {
   void FillPi0Candidate(const Double_t mgg, Int_t absIdMax, Int_t nm);
   // Define CC
   void FitAllSMs();    // SM0 now
-  // Service routine
+  // 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); 
+  static AliEMCALRecPoint *GetRecPoint(AliESDCaloCluster *cl,AliEMCALCalibCoefs *tOld,AliEMCALCalibCoefs *tNew, 
+  TList *l=0, Double_t deff=-1., Double_t w0=-1., Double_t phiSlope=0.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");
@@ -71,7 +74,7 @@ class AliEMCALFolder : public TObjectSet {
 
  protected:
   TList* BookHists();
-  Int_t fCounter; // Coonter of iteration 
+  Int_t fCounter; // Counter of iteration 
  //
   AliEMCALGeometry *fGeometry; //
   //
@@ -89,8 +92,9 @@ class AliEMCALFolder : public TObjectSet {
   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 
+  static const TString fgkDirOfRootFiles;  // name of directory for saving EMCAL folder
 
-  void TestSMStruct();  
+  void TestSMStruct();  // *MENU*
 
   ClassDef(AliEMCALFolder,2) // EMCAL folder
     
index 03a4e48..e6072c1 100644 (file)
@@ -35,6 +35,8 @@
 #include <TH1.h>
 #include <TH2.h>
 #include <TF1.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
 #include <TLatex.h>
 #include <TChain.h>
 #include <TList.h>
@@ -47,6 +49,7 @@
 
 #include "AliESDCaloCluster.h"
 #include "AliEMCALRecPoint.h"
+#include "AliRunLoader.h"
 
 using namespace std;
 
@@ -85,7 +88,7 @@ void AliEMCALHistoUtilities::FillH1(TList *l, Int_t ind, Double_t x, Double_t w)
   //fill 1d histogram
   static TH1* hid=0;
   if(l == 0) return;
-  if(ind < l->GetSize()){
+  if(ind>=0 && ind < l->GetSize()){
     hid = (TH1*)l->At(ind);
     hid->Fill(x,w);
   }
@@ -96,7 +99,7 @@ void AliEMCALHistoUtilities::FillH2(TList *l, Int_t ind, Double_t x, Double_t y,
   //fill 2d histogram
   static TH2* hid=0;
   if(l == 0) return;
-  if(ind < l->GetSize()){
+  if(ind>=0 && ind < l->GetSize()){
     hid = (TH2*)l->At(ind);
     hid->Fill(x,y,w);
   }
@@ -160,7 +163,17 @@ void AliEMCALHistoUtilities::AddToNameAndTitleToList(TList *l, const char *name,
   }
 }
 
-TLatex *AliEMCALHistoUtilities::lat(const char *text, Float_t x,Float_t y, Int_t align, Float_t tsize, short tcolor)
+void AliEMCALHistoUtilities::ResetListOfHists(TList *l)
+{
+  if(l==0) return;
+
+  for(int i=0; i<l->GetSize(); i++) {
+    TH1F* h = (TH1F*)l->At(i);
+    h->Reset(); 
+  }
+}
+
+TLatex *AliEMCALHistoUtilities::Lat(const char *text, Float_t x,Float_t y, Int_t align, Float_t tsize, short tcolor)
 { 
   double y1=y;
   TLatex *latex = new TLatex;
@@ -173,6 +186,138 @@ TLatex *AliEMCALHistoUtilities::lat(const char *text, Float_t x,Float_t y, Int_t
   return latex;
 }
 
+TGraph *AliEMCALHistoUtilities::DrawGraph(Int_t n, Double_t *x, Double_t *y, Int_t markerColor, 
+Int_t markerStyle, const char* opt, const char* tit, const char* xTit,const char* yTit, Int_t ifun,  
+const char *optFit, const char *fun)
+{
+  /* Drawing options 
+  chopt='L' :  A simple polyline between every points is drawn
+  chopt='F' :  A fill area is drawn ('CF' draw a smooth fill area)
+  chopt='A' :  Axis are drawn around the graph
+  chopt='C' :  A smooth Curve is drawn
+  chopt='*' :  A Star is plotted at each point
+  chopt='P' :  Idem with the current marker
+  chopt='B' :  A Bar chart is drawn at each point
+  chopt='1' :  ylow=rwymin
+  chopt='X+' : The X-axis is drawn on the top side of the plot.
+  chopt='Y+' : The Y-axis is drawn on the right side of the plot.
+
+    Fitting options
+   The list of fit options is given in parameter option.
+      option = "W"  Set all errors to 1
+             = "U" Use a User specified fitting algorithm (via SetFCN)
+             = "Q" Quiet mode (minimum printing)
+             = "V" Verbose mode (default is between Q and V)
+             = "B" Use this option when you want to fix one or more parameters
+                   and the fitting function is like "gaus","expo","poln","landau".
+             = "R" Use the Range specified in the function range
+             = "N" Do not store the graphics function, do not draw
+             = "0" Do not plot the result of the fit. By default the fitted function
+                   is drawn unless the option"N" above is specified.
+             = "+" Add this new fitted function to the list of fitted functions
+                   (by default, any previous function is deleted)
+             = "C" In case of linear fitting, not calculate the chisquare
+                    (saves time)
+             = "F" If fitting a polN, switch to minuit fitter
+             = "ROB" In case of linear fitting, compute the LTS regression
+                     coefficients (robust(resistant) regression), using
+                     the default fraction of good points
+               "ROB=0.x" - compute the LTS regression coefficients, using
+                           0.x as a fraction of good points
+
+  */
+  printf("AliEMCALHistoUtilities::drawGraph started \n");
+  printf("Drawing opt |%s| : Fitting opt |%s|\n", opt, optFit);
+
+    TGraph *gr = new TGraph(n, x, y);
+    gr->SetMarkerColor(markerColor);
+    gr->SetLineColor(markerColor);
+    gr->SetMarkerStyle(markerStyle);
+    gr->SetTitle(tit);
+    gr->Draw(opt);
+
+    TString ctmp(opt);
+    if(ctmp.Contains("A")) {
+       gr->GetHistogram()->SetXTitle(xTit);
+       gr->GetHistogram()->SetYTitle(yTit);
+    }
+    ctmp = optFit; 
+    if(ifun>=0) {
+      TString sf("pol"); sf += ifun;
+      gr->Fit(sf.Data(),optFit);
+      printf("\n ** Fit by Polynomial of degree %i : %s **\n", ifun, sf.Data());
+    } else if(!ctmp.Contains("no",TString::kIgnoreCase)){
+      printf("\n ** ifun %i : %s **\n", ifun, fun);
+      gr->Fit(fun, optFit);
+    }
+
+    return gr;
+}
+
+TGraphErrors *AliEMCALHistoUtilities::DrawGraphErrors(const Int_t n,Double_t *x,Double_t *y,Double_t *ex, 
+Double_t *ey, Int_t markerColor,  Int_t markerStyle, const char* opt, const char* tit, 
+const char* xTit,char* yTit, Int_t ifun, const char *optFit, const char *fun)
+{
+  printf("AliEMCALHistoUtilities::drawGraphErrors started \n");
+  printf("Drawing opt |%s| : ifun %i: Fitting opt |%s|, fun |%s|\n", 
+        opt, ifun, optFit, fun);
+
+  TGraphErrors *gr = new TGraphErrors(n, x,y,ex,ey);
+  gr->SetMarkerColor(markerColor);
+  gr->SetLineColor(markerColor);
+  gr->SetMarkerStyle(markerStyle);
+  if(tit&&strlen(tit)>0) gr->SetTitle(tit);
+
+  TString ctmp(opt);
+  if(ctmp.Contains("A")) {
+     gr->GetHistogram()->SetXTitle(xTit);
+     gr->GetHistogram()->SetYTitle(yTit);
+  }
+  if(ifun>0) {
+    if(ifun != 999) {
+      TString sf("pol"); sf += ifun;
+      gr->Fit(sf.Data(),optFit);
+      printf("\n ** Fit by Polynomial of degree %i : %s **\n", ifun, sf.Data());
+    } else {
+      gr->Fit(fun, optFit);
+      printf("\n ** Fit by %s **\n", fun);
+    }
+  } else {
+    if(strlen(optFit)) {
+      printf("\n ** ifun %i : %s **\n", ifun, fun);
+      gr->Fit(fun, optFit);
+    }
+  }
+
+  gr->Draw(opt);
+
+  return gr;
+}
+
+TF1* AliEMCALHistoUtilities::GetResolutionFunction(const char *opt, TString &latexName)
+{
+  TString OPT(opt);
+  OPT.ToUpper();
+  TF1 *fres=0;
+  if      (OPT.Contains("FRES1")) {
+    fres = new TF1("fres","[0]+[1]/sqrt(x)", 0.0, 101.);
+    latexName = "#frac{#sigma_{E}}{E} = A+#frac{B}{#sqrt{E}}";
+  } else if(OPT.Contains("FRES2")) { 
+    fres = new TF1("fres","sqrt([0]*[0]+[1]*[1]/x)", 0.0, 101.);
+    latexName = "#sqrt{A^{2}+#frac{B^{2}}{E}}";
+  }
+  if(fres) {
+    fres->SetParName(0,"A");
+    fres->SetParName(1,"B");
+
+    fres->SetParameter(0, 2.0);
+    fres->SetParameter(1, 6.6);
+    fres->SetLineWidth(2);
+    fres->SetLineColor(kRed);
+  }
+  return fres;
+}
+
 void AliEMCALHistoUtilities::InitChain(TChain *chain, const char* nameListOfFiles, Int_t nFileMax)
 {
   // Read name of files from text file with nameListOfFiles and added to the chain
@@ -199,6 +344,52 @@ void AliEMCALHistoUtilities::InitChain(TChain *chain, const char* nameListOfFile
   //  chain->Lookup();
 }
 
+AliRunLoader* AliEMCALHistoUtilities::InitKinematics(const Int_t nev, const char* galiceName)
+{
+  static AliRunLoader *RL = 0;
+  if(RL == 0 || nev%1000==0) {
+    if(RL)  {
+      RL->UnloadgAlice();
+      delete RL;
+    }
+    RL = AliRunLoader::Open(galiceName,AliConfig::GetDefaultEventFolderName(),"read");
+    RL->LoadgAlice(); // obligatory
+  }
+  if(RL) {
+    RL->GetEvent(nev%1000);
+    RL->LoadKinematics();
+    /* Get what you need after that
+      RL->LoadHits();
+      AliStack *stack=RL->Stack();
+      AliESDCaloCluster * clus = esd->GetCaloCluster(n);
+      Int_t label = clus->GetLabel(); // what is this 
+      TParticle *primary=stack->Particle(label); 
+    */
+  }
+  return RL;
+}
+
+Double_t AliEMCALHistoUtilities::GetMomentum(const char* nameListOfFiles)
+{
+  // Get momentum value from string  - like /....100GEV/.. 
+  Double_t p=-1; // negative if undefined 
+  TString st(nameListOfFiles);
+  if(st.Length()>=5) {
+    st.ToUpper();
+    Ssiz_t ind1 = st.Index("GEV"), ind2=0;
+    if(ind1>0) {
+      for(Int_t i=2; i<=10; i++) {
+        ind2 = st.Index("/",ind1-i);
+        if(ind2>0 && ind2<ind1) break;
+      }
+      TString mom  = st(ind2+1, ind1-ind2-1);
+      if(mom.IsFloat()) p = mom.Atof();
+      printf(" dir |%s| : mom |%s| : p %f : ind2,1 %i->%i\n", st.Data(), mom.Data(), p, ind2, ind1);
+    }
+  }
+  return p;
+}
+
 int AliEMCALHistoUtilities::ParseString(const TString &topt, TObjArray &Opt)
 { 
   // Moved from AliEMCALGeometry
@@ -217,6 +408,7 @@ int AliEMCALHistoUtilities::ParseString(const TString &topt, TObjArray &Opt)
   }
   return Opt.GetEntries();
 }
+
 // Analysis utilites
 Bool_t AliEMCALHistoUtilities::GetLorentzVectorFromESDCluster(TLorentzVector &v, const AliESDCaloCluster* cl)
 {
@@ -352,4 +544,34 @@ Double_t AliEMCALHistoUtilities::GiPol2(Double_t *x, Double_t *par)
   return f;
 }
 
+// Calibration stuff
+Double_t AliEMCALHistoUtilities::GetCorrectionCoefficientForGamma_1(const Double_t eRec)
+{
+  // Correction to rec.energy - Jul 15, 2007
+  // E(gamma) = corr * E(eRec);
+  /* corr = p0*(eRec-7.5)+p1*(eRec-7.5)*(eRec-7.5)+p2; 0.0<eRec<10.0
+   1  p0           6.07157e-05   1.15179e-04  -0.00000e+00   1.20997e-03
+   2  p1           1.50019e-04   3.13566e-05  -0.00000e+00   1.22531e-02
+   3  p2           9.99019e-01   4.08251e-04  -0.00000e+00   1.40606e-03
+
+     corr = p3 + p4*eRec + p5*eRec*eRec
+   1  p3           9.97135e-01   5.31970e-04   1.37962e-09   1.68120e-08
+   2  p4           3.15740e-04   2.53371e-05   1.11475e-11   1.74192e-04
+   3  p5          -1.35383e-06   2.19495e-07  -5.82864e-13   4.52277e-02
+  */
+  static Double_t p0=6.07157e-05, p1=1.50019e-04, p2=9.99019e-01;
+  static Double_t p3=9.97135e-01, p4=3.15740e-04, p5=-1.35383e-06;
+  static Double_t corr=1.0;
+  if(eRec>=0.0 && eRec <=10.0) {
+    corr = p0*(eRec-7.5) + p1*(eRec-7.5)*(eRec-7.5) + p2;
+  } else if(eRec>10.0 && eRec <=105.0) {
+    corr = p3 + p4*eRec + p5*eRec*eRec;
+  }
+  //printf(" eRec %f | corr %f \n", eRec, corr);
+  return corr;
+}
 
+Double_t AliEMCALHistoUtilities::GetCorrectedEnergyForGamma_1(const Double_t eRec)
+{
+  return GetCorrectionCoefficientForGamma_1(eRec) * eRec;
+}
index 1a22ba5..35ed776 100644 (file)
@@ -1,6 +1,6 @@
-#ifndef AliEMCALHistoUtilities_H
-#define AliEMCALHistoUtilities_H
-/* Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
+#ifndef ALIEMCALHISTOUTILITIES_H
+#define ALIEMCALHISTOUTILITIES_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice     */
 
 /* $Id$ */
 #include <TNamed.h>
 
 class TList;
+class TString;
 class TH1;
+class TGraph;
+class TGraphErrors;
 class TF1;
 class TLatex;
 class TChain;
@@ -24,10 +27,10 @@ class TLorentzVector;
 
 class AliESDCaloCluster;
 class AliEMCALRecPoint;
+class AliRunLoader;
 
 class AliEMCALHistoUtilities: public TNamed {
-  public:
-  AliEMCALHistoUtilities(const char *name="emcalUtilitiesRoutines",
+  public:  AliEMCALHistoUtilities(const char *name="emcalUtilitiesRoutines",
   const char *tit="EMCAL utility routines");
   AliEMCALHistoUtilities(const  AliEMCALHistoUtilities &) : TNamed("", ""){
     Fatal("cpy ctor", "not implemented") ; }
@@ -41,10 +44,21 @@ class AliEMCALHistoUtilities: public TNamed {
   const char* opt="RECREATE");
   static void AddToNameAndTitle(TH1   *h=0, const char *name=0, const char *title=0);
   static void AddToNameAndTitleToList(TList *l=0, const char *name=0, const char *title=0);
-  static TLatex *lat(const char *text="", Float_t x=0.0,Float_t y=0.0, Int_t align=12, Float_t tsize=0.05, short tcolor = 1); 
+  static void ResetListOfHists(TList *l);
+
+  static TLatex *Lat(const char *text="", Float_t x=0.0,Float_t y=0.0, Int_t align=12, Float_t tsize=0.05, short tcolor = 1); 
+  static TGraph *DrawGraph(Int_t n=4, Double_t *x=0, Double_t *y=0, Int_t markerColor=4,  
+  Int_t markerStyle=4, const char* opt="", const char* tit="", const char* xTit="  jet E_{t}  ",
+  const char* yTit="", Int_t ifun=0, const char *optFit="W+", const char *fun="");
+  static TGraphErrors *DrawGraphErrors(const Int_t n=4,Double_t *x=0,Double_t *y=0,Double_t *ex=0, 
+  Double_t *ey=0, Int_t markerColor=4,Int_t markerStyle=4, const char* opt="", 
+  const char* tit="", const char* xTit="  jet E_{t}  ",
+  char* yTit="", Int_t ifun=0, const char *optFit="W+", const char *fun="");
   // TChain
   static void InitChain(TChain *chain=0, const char* nameListOfFiles=0, Int_t nFileMax=0); 
-  // 
+  static AliRunLoader* InitKinematics(const Int_t nev=0, const char* galiceName="galice.root");
+  //
+  static Double_t GetMomentum(const char* nameListOfFiles); 
   static int ParseString(const TString &topt, TObjArray &Opt); 
   // Analysis utilites
   static Bool_t GetLorentzVectorFromESDCluster(TLorentzVector &v, const AliESDCaloCluster *cl);
@@ -60,6 +74,10 @@ class AliEMCALHistoUtilities: public TNamed {
   //
   static Double_t Gi(Double_t *x, Double_t *par);
   static Double_t GiPol2(Double_t *x, Double_t *par);
+  // Calibration stuff
+  static Double_t GetCorrectionCoefficientForGamma_1(const Double_t eRec);
+  static Double_t GetCorrectedEnergyForGamma_1(const Double_t eRec);
+  static TF1* GetResolutionFunction(const char *opt, TString &latexName);
 
   AliEMCALHistoUtilities & operator = (const AliEMCALHistoUtilities &) {
     Fatal("operator =", "not implemented") ; return *this ; }
@@ -67,4 +85,4 @@ class AliEMCALHistoUtilities: public TNamed {
   ClassDef(AliEMCALHistoUtilities,1) // EMCAL utility routines
 };
 
-#endif // AliEMCALHistoUtilities_H
+#endif // ALIEMCALHISTOUTILITIES_H
index 84003f4..8d5f865 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Log$ */
+/* $Log$co: warning: `/* $Log' is obsolescent; use ` * $Log'.
+
+ * Revision 1.1  2007/08/08 15:58:01  hristov
+ * New calibration classes. They depend on TTable, so libTable.so is added to the list of Root libraries. (Aleksei)
+ * */
 
 //_________________________________________________________________________
 //    Set of parameters for pi0 selection 
 
 #include "AliEMCALPi0SelectionParam.h"
 
-TableClassImpl(AliEMCALPi0SelectionParam,pi0SelectionParam)
+ClassImp(pi0SelectionParam)
+pi0SelectionParam::pi0SelectionParam() : 
+eOfRpMin(0.3), eOfRpMax(30.), massGGMin(0.03), massGGMax(0.28), momPi0Min(1.8), momPi0Max(12.)
+{ }
+//_________________________________________________________________________
+
+ClassImp(AliEMCALPi0SelectionParam)
+AliEMCALPi0SelectionParam::AliEMCALPi0SelectionParam() : TNamed("",""), fTable(0), fCurrentInd(0)
+{
+}
+
+AliEMCALPi0SelectionParam::AliEMCALPi0SelectionParam(const char* name, const Int_t nrow) : TNamed(name,"table of cell information") , fTable(0), fCurrentInd(0)
+{
+  fTable = new TObjArray(nrow);
+}
+
+void AliEMCALPi0SelectionParam::AddAt(pi0SelectionParam* r)
+{
+  (*fTable)[fCurrentInd] = new pi0SelectionParam(*r);
+  fCurrentInd++;
+}
+
+AliEMCALPi0SelectionParam::~AliEMCALPi0SelectionParam()
+{
+  if(fTable) {
+    fTable->Delete();
+    delete fTable;
+  }
+}
+
+pi0SelectionParam* AliEMCALPi0SelectionParam::GetTable(Int_t i) const
+{
+  return (pi0SelectionParam*)fTable->At(i);
+}
 
 void AliEMCALPi0SelectionParam::PrintTable()
 {
index 787c31a..121baf1 100644 (file)
@@ -6,26 +6,47 @@
 /* $Id$ */
 
 //_________________________________________________________________________
-//    Set of parameters for pi0 selection 
+//  Set of parameters for pi0 selection 
 //                  
 //*-- Author: Aleksei Pavlinov (WSU, Detroit, USA) 
 
 // --- ROOT system ---
-#include <TTable.h>
+#include <TNamed.h>
+#include <TObjArray.h>
 
 // unit is GeV
-struct  pi0SelectionParam {
+class  pi0SelectionParam :  public TObject{
+ public:
+  pi0SelectionParam();
+  virtual ~pi0SelectionParam() {};
+  virtual const char* GetName() const {return "Pi0Par";}
+
   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
+
+  ClassDef(pi0SelectionParam,1);
 };
 
-class AliEMCALPi0SelectionParam : public TTable {
+class AliEMCALPi0SelectionParam : public TNamed {
  public:
-  ClassDefTable(AliEMCALPi0SelectionParam , pi0SelectionParam)
+  AliEMCALPi0SelectionParam(); // default constractor
+  AliEMCALPi0SelectionParam(const char* name, const Int_t nrow);
+  virtual ~AliEMCALPi0SelectionParam();
+
+  AliEMCALPi0SelectionParam & operator = (const AliEMCALPi0SelectionParam  & /*rvalue*/) {
+    // assignement operator requested by coding convention but not needed
+    Fatal("operator =", "not implemented");
+    return *this;
+  };
+  // 
+  void AddAt(pi0SelectionParam* r);
+  pi0SelectionParam* GetTable(Int_t i) const;
+  Int_t       GetSize()  const {return fTable->GetSize();}
+  Int_t       GetNRows() const {return fCurrentInd;}
 
  // Menu
   void PrintTable();                 // *MENU*
@@ -34,8 +55,12 @@ class AliEMCALPi0SelectionParam : public TTable {
 
   // Set of parameter(s)
   static AliEMCALPi0SelectionParam* Set1();
+  //
+ protected:
+  TObjArray *fTable;
+  Int_t fCurrentInd;
 
-  ClassDef(AliEMCALPi0SelectionParam,1) // Set of Parameters For Pi0 Selection     
+  ClassDef(AliEMCALPi0SelectionParam, 2) // Set of Parameters For Pi0 Selection     
 };
 
 #endif // ALIEMCALPI0SELECTIONPARAM_H
index 05d3b41..7881e07 100644 (file)
@@ -1,5 +1,5 @@
 /**************************************************************************
- * Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
+ * 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.              *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Log$ */
-
+/* $Log$co: warning: `/* $Log' is obsolescent; use ` * $Log'.
+ * Revision 1.1  2007/08/08 15:58:01  hristov
+ * New calibration classes. (Aleksei)
+ * */
 //_________________________________________________________________________
 // This is selector for making a set of histogramms from ESD for rec.points
 //*--  Authors: Aleksei Pavlinov (WSU) 
-//*--  Feb 17, 2007 - May 23, 2007
+//*--  Feb 17, 2007 - Sep 11, 2007
 //*--  Pi0 calibration
+//*--  Recalibration, linearity, shower profile 
 //_________________________________________________________________________
 
 #include "AliEMCALRecPointsQaESDSelector.h"
 #include "AliEMCALFolder.h"
 #include "AliEMCALSuperModule.h"
 #include "AliEMCALCell.h"
+#include "AliEMCALCellInfo.h"
 #include "AliEMCALPi0SelectionParam.h"
 #include "AliEMCALHistoUtilities.h"
 #include "AliEMCALCalibCoefs.h"
 #include "AliLog.h"
 #include "AliESD.h"
 #include "AliEMCALRecPoint.h"
-
+#include "AliRunLoader.h"
+#include "AliStack.h"
+#include "jetfinder/AliEMCALJetMicroDst.h" // Sgpdge
+#include "AliEMCALDigit.h"
+
+#include <iostream>
+#include <iomanip>
+#include <fstream>
 #include <assert.h>
 #include <map>
 #include <string>
 #include <TCanvas.h>
 #include <TVector3.h>
 #include <TLorentzVector.h>
+#include <TParticle.h>
 #include <TChain.h>
 #include <TFile.h>
 #include <TH1F.h>
+#include <TH1D.h>
+#include <TH2F.h>
 #include <TF1.h>
 #include <TMath.h>
-#include <TArrayS.h>
 #include <TBrowser.h>
-#include <TDataSet.h>
+#include <TList.h>
+#include <TCanvas.h>
+#include <TLine.h>
+#include <TObjString.h>
+#include <TArrayI.h>
+#include <TGraphErrors.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TNtuple.h>
 
 using namespace std;
 
 typedef  AliEMCALHistoUtilities u;
 
+Double_t AliEMCALRecPointsQaESDSelector::fDistEff = 3.;
+Double_t AliEMCALRecPointsQaESDSelector::fW0      = 4.5;
+Double_t AliEMCALRecPointsQaESDSelector::fSlopePhiShift = 0.01; // ?? just guess
+
+AliEMCALFolder* AliEMCALRecPointsQaESDSelector::fEMCAL = 0;
+AliEMCALFolder* AliEMCALRecPointsQaESDSelector::fEMCALOld = 0;
+
+AliEMCALGeometry *geo=0;
+
+char *anaOpt[]={
+  "CORR1",   // GetCorrectedEnergyForGamma_1(Double_t eRec);
+  "RECALIB",
+  "IDEAL",
+  "PI0",
+  "GAMMA",
+  "KINE",   // reading kine file
+  "PROF",   // Shower profile: phi direction now
+  "FIT"     // define parameters : deff, w0 and phislope
+};
+enum keyOpt{
+  kCORR1,
+  kRECALIB,
+  kIDEAL,
+  kPI0,
+  kGAMMA,
+  kKINE,
+  kPROF,
+  kFIT,
+ //-- 
+  kEND
+};
+        
+int  nAnaOpt = sizeof(anaOpt) / sizeof(char*); 
+// Enumeration variables
+// enum {kUndefined=-1, kLed=0, kBeam=1};
+// enum {kPatchResponse=0, kClusterResponse=1};
+// enum {kSkip=0,  kSaveToMemory=1};
+
 ClassImp(AliEMCALRecPointsQaESDSelector)
 
 Int_t  nmaxCell = 4*12*24*11;
 
 AliEMCALRecPointsQaESDSelector::AliEMCALRecPointsQaESDSelector() :
   AliSelector(),
+  fPmom(0),
   fChain(0),
   fLofHistsPC(0),
   fLofHistsRP(0),
+  fLKineVsRP(0),
+  fLShowerProfile(0),
+  fCellsInfo(0),
   fEmcalPool(0),
-  fEMCAL(0),
-  fEMCALOld(0)
+  fRunOpts(0),
+  fArrOpts(0),
+  fKeyOpts(0)
 {
   //
   // Constructor. Initialization of pointers
@@ -92,13 +156,6 @@ void AliEMCALRecPointsQaESDSelector::Begin(TTree* tree)
   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.
@@ -110,19 +167,92 @@ void AliEMCALRecPointsQaESDSelector::SlaveBegin(TTree* tree)
   AliSelector::SlaveBegin(tree);
 }
 
+void AliEMCALRecPointsQaESDSelector::Init(TTree* tree)
+{
+  // Init function
+
+  AliSelector::Init(tree);
+}
+Bool_t AliEMCALRecPointsQaESDSelector::Notify()
+{
+  // The Notify() function is called when a new file is opened. This
+  // can be either for a new TTree in a TChain or when when a new TTree
+  // is started when using PROOF. Typically here the branch pointers
+  // will be retrieved. It is normaly not necessary to make changes
+  // to the generated code, but the routine can be extended by the
+  // user if needed.
+
+  return AliSelector::Notify();
+}
+
 void AliEMCALRecPointsQaESDSelector::InitStructure(Int_t it)
 {
-  AliEMCALGeometry::GetInstance("SHISH_TRD1_CURRENT_2X2"); // initialize geometry just once
-  //AliEMCALGeometry::GetInstance(""); // initialize geometry just once
+  geo = AliEMCALGeometry::GetInstance("SHISH_TRD1_CURRENT_2X2"); // initialize geometry just once
+  fCellsInfo = AliEMCALCellInfo::GetTableForGeometry(geo);
+
+  if(fRunOpts.Length()>0) CheckRunOpts();
 
-  fLofHistsPC = DefineHistsOfRP("PseudoCl");
-  fLofHistsRP = DefineHistsOfRP("RP");
+  Int_t key = 0;
+  if(GetKeyOptsValue(kGAMMA)) key = 1;
+  fLofHistsRP = DefineHistsOfRP("RP", fPmom, key);
+  fLofHistsPC = DefineHistsOfRP("PseudoCl", fPmom);
 
-  fEmcalPool = new TObjectSet("PoolOfEMCAL");
-  if(it == 1) {
+  fEmcalPool = new TFolder("PoolOfEMCAL","");
+  if(it <= 1) {
     fEMCAL     = new AliEMCALFolder(it); // folder for first itteration   
     fEmcalPool->Add(fEMCAL);
   }
+  //if(it<=0) SetName("GammaSel"); // For convinience
+  //else      SetName("Pi0Sel");
+  if(GetKeyOptsValue(kKINE)) fLKineVsRP = DefineHistsOfKineVsRP("KineVsRP",fPmom, key);
+  if(GetKeyOptsValue(kPROF)) fLShowerProfile = DefineHistsForShowerProfile("ProfY", fPmom);
+}
+
+void AliEMCALRecPointsQaESDSelector::CheckRunOpts()
+{
+  fRunOpts.ToUpper();
+  int nopt = u::ParseString(fRunOpts, fArrOpts);
+  printf("<I> AliEMCALRecPointsQaESDSelector::CheckRunOpts() analyze %i(%i) options : nAnaOpt %i\n", 
+         nopt, fArrOpts.GetEntries(), nAnaOpt);  
+  if(nopt <=0) return;
+
+  fKeyOpts = new TArrayI(nAnaOpt);
+  for(Int_t i=0; i<nAnaOpt; i++) (*fKeyOpts)[i] = 0;
+
+  for(Int_t i=0; i<fArrOpts.GetEntries(); i++ ) {
+    TObjString *o = (TObjString*)fArrOpts.At(i); 
+
+    TString runOpt = o->String();
+    Int_t indj=-1;
+
+    for(Int_t j=0; j<nAnaOpt; j++) {
+      TString opt = anaOpt[j];
+      if(runOpt.Contains(opt,TString::kIgnoreCase)) {
+       indj = j;
+        break;
+      }
+    }
+    if(indj<0) {
+      printf("<E> option |%s| unavailable **\n", 
+      runOpt.Data());
+      assert(0);
+    } else {
+      (*fKeyOpts)[indj] = 1;
+      printf("<I> option |%s| is valid : number %i : |%s| \n", 
+            runOpt.Data(), indj, anaOpt[indj]);
+    }
+  }
+}
+
+Int_t AliEMCALRecPointsQaESDSelector::GetKeyOptsValue(Int_t key)
+{
+  static Int_t val=0;
+  val = 0;
+  if(fKeyOpts && key>=0 && key<fKeyOpts->GetSize()) {
+    val = fKeyOpts->At(key);
+  }
+  // printf(" key %i : val %i : opt %s \n", key, val, anaOpt[key]);
+  return val;
 }
 
 Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
@@ -145,8 +275,13 @@ Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
   //  Assuming that fTree is the pointer to the TChain being processed,
   //  use fTree->GetTree()->GetEntry(entry).
 
-  if (AliSelector::Process(entry) == kFALSE)
+  if (AliSelector::Process(entry) == kFALSE) {
+    AliDebug(AliLog::kError, Form("Entry %i is not available ",entry));
     return kFALSE;
+  }
+
+  // esd->ReadFromTree(tree);
+  // AliESDEvent * esd = new AliESDEvent();
 
   // Check prerequisites
   if (!fESD)
@@ -154,6 +289,7 @@ Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
     AliDebug(AliLog::kError, "ESD branch not available");
     return kFALSE;
   }
+  
   static pi0SelectionParam* rPar = GetEmcalFolder()->GetPi0SelectionParRow(0);
 
   static Int_t nEmcalClusters, indOfFirstEmcalRP, nEmcalRP,nEmcalPseudoClusters;
@@ -161,6 +297,9 @@ Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
   indOfFirstEmcalRP = fESD->GetFirstEMCALCluster();
   u::FillH1(fLofHistsRP, 1, double(indOfFirstEmcalRP));
 
+  static AliRunLoader* RL = 0;
+  static Int_t nev = 0; // Temporary - 0nly for reading one file now !!
+  
   static AliESDCaloCluster *cl = 0; 
   nEmcalRP = nEmcalPseudoClusters = 0;
   TList *l=0;
@@ -171,6 +310,7 @@ Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
 
   static TLorentzVector v, vcl;
   int nrp = 0; // # of RP for gg analysis
+  static Double_t erec=0., ecorr=0.0;
   for(int i=indOfFirstEmcalRP; i<indOfFirstEmcalRP+nEmcalClusters; i++) {
     cl = fESD->GetCaloCluster(i);
     if(cl->GetClusterType() == AliESDCaloCluster::kPseudoCluster) {
@@ -178,10 +318,26 @@ Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
       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(fEMCAL->GetIterationNumber()>1||GetKeyOptsValue(kIDEAL)||GetKeyOptsValue(kRECALIB)||GetKeyOptsValue(kFIT)) {
+        AliEMCALRecPoint *rp=0;
+        if(GetKeyOptsValue(kFIT) == kFALSE) fDistEff = -1.; // No fitting ; Sep 4, 2007
+        if(GetKeyOptsValue(kIDEAL)) {
+          rp = AliEMCALFolder::GetRecPoint(cl, fEMCAL->GetCCFirst(), 0, fLofHistsRP, fDistEff, fW0, fSlopePhiShift);
+       } else {
+          rp = AliEMCALFolder::GetRecPoint(cl, fEMCAL->GetCCFirst(), fEMCAL->GetCCIn(), fLofHistsRP, fDistEff, fW0, fSlopePhiShift);
+       }
+        if(GetKeyOptsValue(kPROF)) {
+          FillHistsForShowerProfile( GetListShowerProfile(), rp, GetCellsInfo());
+       }
+       //if(rp->GetPointEnergy()>=rPar->eOfRpMin && u::GetLorentzVectorFromRecPoint(v, rp)) {
         if(u::GetLorentzVectorFromRecPoint(v, rp)) { // comparing with RP
+          if(GetKeyOptsValue(kCORR1)) {
+           erec  = v.Rho();
+            ecorr = u::GetCorrectedEnergyForGamma_1(erec);
+            v.SetRho(ecorr);
+            v.SetE(ecorr); // This is gamma
+           //printf("<1> erec %f | ecorr %f \n", erec, ecorr);
+         }
           new(lvM1[nrp]) TLorentzVector(v);
           indLv[nrp] = i;
           nrp++;
@@ -190,13 +346,25 @@ Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
          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()));
+
+         u::FillH1(fLofHistsRP, 16, v.P());
+         u::FillH2(fLofHistsRP, 17, vcl.P(), vcl.P()-v.P());
           l = 0; // no filling
+          if(GetKeyOptsValue(kIDEAL) || GetKeyOptsValue(kRECALIB)) l = fLofHistsRP;
         }
-        delete rp;
+        if(rp) 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 ! 
+          if(GetKeyOptsValue(kCORR1)) {
+           erec  = v.Rho();
+            ecorr = u::GetCorrectedEnergyForGamma_1(erec);
+            v.SetRho(ecorr);
+            v.SetE(ecorr); // This is gamma now
+            // printf("<2> erec %f | ecorr %f \n", erec, ecorr);
+            // printf("    v.Rho()  %f \n", v.Rho());
+         }
           new(lvM1[nrp]) TLorentzVector(v);
           indLv[nrp] = i;
           nrp++;
@@ -235,6 +403,10 @@ Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
   u::FillH1(fLofHistsPC, 0, double(nEmcalPseudoClusters));
 
   static TLorentzVector *lv1=0, *lv2=0, lgg;
+  for(int i1=0; i1<nrp; i1++){
+    lv1 = (TLorentzVector*)lvM1.At(i1);
+    u::FillH1(fLofHistsRP, 18, lv1->P());
+  }
   static Double_t mgg, pgg;
   mgg = pgg = 0.;
   nrp = lvM1.GetEntriesFast(); 
@@ -251,21 +423,40 @@ Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
  
        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());
+            if(fEMCAL && fEMCAL->GetIterationNumber()>=1) {
+              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());
+           }
          }
        }
       }
     }
   }
 
+  //  static Int_t fileNumber = 0;
+  static TString curFileName;
+  //  if(GetKeyOptsValue(kKINE) && nev<fChain->GetEntries()) {
+  if(GetKeyOptsValue(kKINE)) {
+  // Get galice.root file in current directory
+    if(nev%1000==0) {
+      printf(" current file |%s|\n", fChain->GetCurrentFile()->GetName());
+      curFileName = fChain->GetCurrentFile()->GetName();
+      curFileName.ReplaceAll("AliESDs.","galice."); 
+    }
+    RL = u::InitKinematics(nev, curFileName.Data());
+  // Compare kineamtics vs EMCal clusters
+    FillHistsOfKineVsRP(fLKineVsRP, RL, lvM1);
+  }
+
   lvM1.Delete();
 
   if(nEmcalClusters != (nEmcalRP+nEmcalPseudoClusters))
     Info("Process","nEmcalClusters %i : RP %i + PC %i ",nEmcalClusters, nEmcalRP, nEmcalPseudoClusters); 
 
+  nev++;
+
   return kTRUE;
 }
 
@@ -294,24 +485,18 @@ void AliEMCALRecPointsQaESDSelector::Terminate()
   // 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)
+TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfRP(const char *name,Double_t p,Int_t keyOpt)
 {
-  double ADCchannelEC = 0.00305; // ~3mev per adc count
+  printf("<I> DefineHistsOfRP :%s : p %f : keyOpt %i \n", name, p, keyOpt);
+  Double_t ADCchannelEC = 0.0153; // ~15mev per adc count
+  Double_t xma = p*1.4, xmi=0.0, step=0.0, xmic=xmi, xmac = xma;
+  if(xma<0) xma = 20.;
+  Int_t nmax=1000, scale=4, nmaxc = nmax;
 
   gROOT->cd();
   TH1::AddDirectory(1);
@@ -320,32 +505,310 @@ TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfRP(const char *name)
 
   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;
+  if(keyOpt==1 && p>0.1) {
+    if   (p>=100.)  scale=12;
+    else if(p>=50.) scale=10;
+    else if(p>=30.) scale=8;  
+    xma = p + scale*(0.15*TMath::Sqrt(p));
+    xmi = p - scale*(0.15*TMath::Sqrt(p));
+    step = (xma-xmi)/nmax;
+  }
+
+  if(step < 0.0153) {
+    nmax = int((xma-xmi) / ADCchannelEC)+1;
+    xma =  xmi + ADCchannelEC*nmax;
+  } 
+  new TH1F("04_EnergyOf", "energy of ", nmax, xmi, xma);
+  nmaxc = nmax; xmic=xmi; xmac = xma;
+
+  nmax = 10000;
+  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("05_DigitEnergyIn", "digit energy in ", nmaxc, xmic, xmac);
   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.);
+  Double_t pmax=11., dpmax=0.20;
+  if(p>0.1) {
+    pmax=p*1.2;
+    dpmax *= p;
+  }
+  new TH1F("11_MomClESD-RpRecalib", "difference of momentum cl(ESD) - rp(Recalib)", 100, -.01, dpmax);
   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);
+  //  AliEMCALGeometry* g = AliEMCALGeometry::GetInstance();
+  new TH1F("15_AbsIdRecalib", "abs Id of digits after recalibration", geo->GetNCells(),-0.5,Double_t(geo->GetNCells())-0.5);
+  new TH1F("16_EnergyOfRecalibRp_", "energy of recalibrated rec.points", nmaxc, xmic, xmac); // Jul 12, 2007
+  new TH2F("17_ShiftRecalib_", "E(clESD) - E(recalib)", 110,0.0, pmax, 50,0.0,dpmax); // Jul 13, 2007
   
-  TList *l = u::MoveHistsToList(Form("ListOfHists%s",name), kFALSE);
-  u::AddToNameAndTitleToList(l, name, name);
+  // Corrected staff
+  new TH1F("18_EnergyOfCorrectedLV", "energy of corrected LV", nmaxc, xmic, xmac);
+
+  TString st = Form("ListOfHists%s_P=%5.1f",name,p);
+  st.ReplaceAll(" ","");
+  TList *l = u::MoveHistsToList(st.Data(), kFALSE);
+  st = Form("%s_P=%5.1f",name,p);
+  st.ReplaceAll(" ","");
+  u::AddToNameAndTitleToList(l, st.Data(), st.Data());
+
+  return l;
+}
+
+TList* AliEMCALRecPointsQaESDSelector::DefineHistsOfKineVsRP(const char *name,  Double_t p, Int_t keyOpt)
+{
+  printf("<I>  DefineHistsOfKineVsRP :%s : p %f : keyOpt %i \n", name, p, keyOpt);
+
+  gROOT->cd();
+  TH1::AddDirectory(1);
+  new TH1F("00_hVx",Form("Vx of primary vertex"), 100, -5., +5.);   // 00
+  new TH1F("01_hVy",Form("Vy of primary vertex"), 100, -5., +5.);   // 01
+  new TH1F("02_hVz",Form("Vz of primary vertex"), 100, -50., +50.); // 02
+
+  //  Double_t ADCchannelEC = 0.0153; // ~15mev per adc count
+  Double_t xma = p*1.4, xmi=0.0, sig=0.15*TMath::Sqrt(p);
+  //  Double_t step=0.0, xmic=xmi, xmac = xma;
+  if(xma<0) xma = 20.;
+  //Int_t nmax=1000;
+  // scale=4, nmaxc = nmax;
+
+  new TH1F("03_hGidPrimar", "Geant Id of primary particle ", 101, 0.5, 100.5);   // 03
+  new TH1F("04_hPmomPrim","momentum of primary particle", 100, xmi, xma);        // 04
+  new TH1F("05_hEtaPrim","#eta of primary particle ", 200, -1., 1.);             // 05
+  new TH1F("06_hPhiPrim","#phi of primary particle ", 63*2, 0.0,TMath::Pi()*2.); // 06
+  new TH1F("07_hThetaPrim","#theta of primary particle", 63, 0.0,TMath::Pi());   // 07 
+  new TH1F("08_hPhiPrimInDegree","#phi of primary particle in degree", 120, 75.0, 195.0);  // 08
+  new TH1F("09_hThetaPrimInDegree","#theta of primary particle in degree", 90, 45., 135.); // 09
+
+  // Particle vs cluster
+  new TH1F("10_hE(P)-E(RP)"," E(p) - E(RP) ", 100, -10.*sig, 10.*sig);   // 10
+
+  new TH1F("11_hPvsRpAngleInDegree","angle between P and RP (in degree)", 110, 0.0, 1.0); // 11
+  double dphi=0.5; // for fitting
+  new TH1F("12_hPvsRpDPhiInDegree","dif(#phi) between P and RP (in degree)", 100, -dphi, dphi); // 12
+  new TH1F("13_hPvsRpDThetaInDegree","dif(#theta) between P and RP (in degree)", 100, -0.5, 0.5); // 13
+
+  new TH1F("14_hPvsRpDPhiInDegree","dif(#phi) between P and RP (in degree) right part of SM", 100, -dphi, +dphi); // 14
+  new TH1F("15_hPvsRpDPhiInDegree","dif(#phi) between P and RP (in degree) left part of SM", 100, -dphi, dphi); // 15
+
+  new TH1F("16_hPvsRpDThetaInDegree","dif(#theta) between P and RP (even index)", 100, -0.5, 0.5); // 16
+  new TH1F("17_hPvsRpDThetaInDegree","dif(#theta) between P and RP (odd  index)", 100, -0.5, 0.5); // 17
+
+  TString st = Form("ListOfHists%s_P=%5.1f",name,p);
+  st.ReplaceAll(" ","");
+  TList *l = u::MoveHistsToList(st.Data(), kFALSE);
+  st = Form("%s_P=%5.1f",name,p);
+  st.ReplaceAll(" ","");
+  u::AddToNameAndTitleToList(l, st.Data(), st.Data());
+  return l;
+}
+
+TList *AliEMCALRecPointsQaESDSelector::DefineHistsForShowerProfile(const char *name, Double_t p)
+{
+  // Aug 1, 2007 - shower profile business as in testbeam analysis
+  //               Phi direction here is Y direction in testbeam 
+  printf("<I> DefineHistsForShowerProfile: %s : p %f \n", 
+  name, p);
+
+  gROOT->cd();
+  TH1::AddDirectory(1);
+
+  // Cell unit
+  new TH1F("00_hPHiMean"," mean in #phi direction ", 24*10, 0.0, 24.);
+  new TH1F("01_hPHiLog"," mean with log in #phi direction ", 24*10, 0.0, 24.);
+  new TH2F("02_XmeanVsXlog", "xmean vs xlog", 240,0.0, 24.,240,0.0, 24.);
+  new TH2F("03_XmeanVsXlog", "xmean vs xlog (system of cell with max energy, first half on #eta)", 
+  50,-0.5,+0.5, 100,-1.,+1.);
+  new TH2F("04_XmeanVsXlog", "xmean vs xlog (system of cell with max energy, second half on #eta)", 
+  50,-0.5,+0.5, 100,-1.,+1.);
+
+  new TH2F("05_PhiShowerProfile", " #phi shower profile - cumulative distribution",
+  6*25,-3.,3., 201, -0.0025, 1.0025);   // 00; x direction in cell unit
+
+  TString st = Form("L%s_P=%5.1f",name,p);
+  st.ReplaceAll(" ","");
+  TList *l = u::MoveHistsToList(st.Data(), kFALSE);
+  st = Form("%s_P=%5.1f",name,p);
+  st.ReplaceAll(" ","");
+  u::AddToNameAndTitleToList(l, st.Data(), st.Data());
 
   return l;
 }
 
+void AliEMCALRecPointsQaESDSelector::FillHistsOfKineVsRP(TList *l, AliRunLoader* RL,  TClonesArray &lvM)
+{
+  //
+  // lvM - array of TLorentzVector's which was cretaef from AliESDCaloCluster's
+  //
+
+  if(l==0 || RL==0) return;
+
+  // TNtuple for qucik analysis
+  static TNtuple *nt=0;
+  if(nt==0) {
+    gROOT->cd();
+    Int_t bsize = int(1.e+7);
+    nt = new TNtuple("angle","angle ntuple for quick analysis",
+    "dPhi:dTheta:phiP:thetaP", bsize);
+    gROOT->GetListOfBrowsables()->Add(nt);
+  }
+  static AliStack* st=0;
+  static TParticle *p=0;
+  static Int_t gid=0, ic=0, pdg=0, i=0;
+  gid = ic = pdg = 0;
+
+  st = RL->Stack();
+  if(st == 0) return;
+  // first primary particle
+  p = st->Particle(0);
+  if(p == 0) return;
+
+  u::FillH1(l, ic++, p->Vx());
+  u::FillH1(l, ic++, p->Vy());
+  u::FillH1(l, ic++, p->Vz());
+
+  pdg = p->GetPdgCode();
+  //gid = gMC->IdFromPDG(pdg); // gMc should be defined
+  AliEMCALJetMicroDst::Sgpdge(pdg,gid); // temporary
+
+  u::FillH1(l, ic++, Double_t(gid));
+  u::FillH1(l, ic++, p->P());
+  u::FillH1(l, ic++, p->Eta());
+  u::FillH1(l, ic++, TVector2::Phi_0_2pi(p->Phi()) );
+  u::FillH1(l, ic++, p->Theta());
+  // In degree
+  u::FillH1(l, ic++, TVector2::Phi_0_2pi(p->Phi())*TMath::RadToDeg());
+  u::FillH1(l, ic++, p->Theta()*TMath::RadToDeg()); // 09
+
+  if(lvM.GetEntriesFast() == 0) return;
+  //
+  // Initial kinematics vs Calo clusters
+  //
+  static TLorentzVector *lv = 0, lvp;
+  lvp.SetPxPyPzE(p->Px(),p->Py(),p->Pz(),p->Energy()); // Particle 
+
+  Double_t angle = 0.0, angleMin = 180.0, eDiff=0.0, dPhi=0., dTheta=0., phiP=0.0, thetaP=0.0;
+  Int_t indMin = 0;
+  phiP   = lvp.Phi() * TMath::RadToDeg();
+  if (phiP<81. || phiP>99.)  return;     // cut phi boundaries
+  thetaP = lvp.Theta() * TMath::RadToDeg();
+  if (thetaP<56. || thetaP>89.)  return; // cut theta boundaries
+
+  for(i=0; i<lvM.GetEntriesFast(); i++) {
+    lv = (TLorentzVector*)lvM.At(i);
+    angle = lv->Angle(lvp.Vect())*TMath::RadToDeg();
+    if(angleMin > angle) {
+      angleMin = angle;
+      indMin   = i;
+    }
+  }
+  lv = (TLorentzVector*)lvM.At(indMin);
+  eDiff   = lvp.E() - lv->E();
+  u::FillH1(l, ic++, eDiff); 
+
+  dPhi    = TVector2::Phi_mpi_pi(lvp.Phi()-lv->Phi())*TMath::RadToDeg();
+  dTheta  = TVector2::Phi_mpi_pi(lvp.Theta()-lv->Theta())*TMath::RadToDeg();
+  u::FillH1(l, ic++, angleMin);
+  u::FillH1(l, ic++, dPhi); 
+  u::FillH1(l, ic++, dTheta); 
+
+  if       (phiP>=81. && phiP<=90.) {
+    u::FillH1(l, 14, dPhi); 
+  } else if(phiP> 90. && phiP<=99.) {
+    u::FillH1(l, 15, dPhi); 
+  }
+  // Unfinished - have to get absid of digit with max energy 
+  u::FillH1(l, 16, dTheta); 
+  u::FillH1(l, 17, dTheta); 
+  if(nt) {
+    nt->Fill(dPhi,dTheta,phiP,thetaP);
+    /*
+     tv__tree = (TTree *) gROOT->FindObject("angle");
+     tv__tree->Draw("dPhi:phiP","abs(dPhi)<0.5","", 9182, 0);
+     tv__tree->Draw("dTheta:thetaP","abs(dTheta)<0.5","", 9182, 0);
+     */
+  }
+}
+
+void AliEMCALRecPointsQaESDSelector::FillHistsForShowerProfile
+(TList *l, AliEMCALRecPoint *rp, AliEMCALCellInfo * t)
+{
+  // Aug 1, 2007
+  if(l==0 || rp==0 || t==0) return;
+  // --
+  static Double_t xmean=0., xlog=0.;
+  static cellInfo rMax;
+  static Int_t phiSize;
+
+  if(rp->GetPointEnergy() < 1.0) return;
+  //if(rp->GetPointEnergy() > 8.0) return; // discard merged clusters for pi0 case
+
+  EvalLocalPhiPosition(1.,  rp, t, xmean, phiSize, rMax);
+  if(phiSize == 0) return; // just one row in cell directions
+
+  EvalLocalPhiPosition(5.5, rp, t, xlog, phiSize, rMax);
+  if(rMax.iPhi>1.5&&rMax.iPhi<21.5 && rMax.iEta>1.5&&rMax.iEta<46.5) { 
+    u::FillH1(l, 0, xmean); 
+    u::FillH1(l, 1, xlog); 
+    u::FillH2(l, 2, xlog, xmean);
+    // Select two central modules
+    if((rMax.iPhim==5 || rMax.iPhim==6)){
+    // Transition to system of cell with max energy
+      xmean -= (double(rMax.iPhi)+0.5);  
+      xlog  -= (double(rMax.iPhi)+0.5);  
+      if(rMax.iEtam>=2 && rMax.iEtam<=12){ // approximatively first half on eta 
+        u::FillH2(l, 3, xlog, xmean);
+      } else {// approximatively second half on eta 
+        u::FillH2(l, 4, xlog, xmean);
+      }
+    }
+  }
+}
+
+void   AliEMCALRecPointsQaESDSelector::EvalLocalPhiPosition(const Double_t wlog, const AliEMCALRecPoint *rp, const AliEMCALCellInfo* t, Double_t &xcog, Int_t &phiSize, cellInfo &rMax)
+{
+  // wlog = 1 - usual center of gravity; >1 - with logarithmic weight.
+  // digits - array of digits
+  // t      - 
+  //==
+  // xcog - position of center of gravity in cell unit 
+  if(wlog<1.0 || rp==0 || t==0) return;
+  xcog = 0;
+
+  static Double_t wtot=0., w=0., edigi=0., e=0., edigiMax=0.;
+  static Int_t absid = 0, phiMin=0, phiMax=0;
+  static cellInfo* r=0;
+
+  e = rp->GetPointEnergy();
+  wtot = 0.0;
+  edigiMax=0.;
+
+  phiMin = 23; phiMax = 0;
+  for(Int_t iDigit=0; iDigit<rp->GetMultiplicity(); iDigit++) {
+    absid = rp->GetAbsId()[iDigit];
+    edigi = rp->GetEnergiesList()[iDigit];
+    if(wlog > 1.0)  w = TMath::Max( 0., wlog + TMath::Log(edigi/e));
+    else            w = edigi; // just energy
+
+    r     = t->GetTable(absid);
+    xcog += w*(Double_t(r->iPhi) + 0.5);
+    wtot += w;
+    if(edigi > edigiMax) {
+      edigiMax = edigi;
+      rMax = (*r);
+    }
+    if(phiMin > r->iPhi) phiMin = r->iPhi; 
+    if(phiMax < r->iPhi) phiMax = r->iPhi; 
+  }
+  xcog /= wtot;
+  phiSize = phiMax - phiMin;
+  //  printf("phiSize %i \n", phiSize); 
+}
 /* unused now
 TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfTowers(const char *name)
 {
@@ -374,14 +837,27 @@ void AliEMCALRecPointsQaESDSelector::FitEffMassHist()
 void AliEMCALRecPointsQaESDSelector::PrintInfo()
 {
   // Service routine
+  printf("\n %i Entrie(s) | Option(s) |%s| \n", GetOptsArray().GetEntries(), fRunOpts.Data());  
+  for(int i=0; i<nAnaOpt; i++) {
+    if(GetKeyOptsValue(i)) printf(" %i |%s| \n", i, anaOpt[i]);
+  }
+
   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()));
   }
+  printf(" fDistEff %f fW0 %f fSlopePhiShift %f \n", fDistEff, fW0, fSlopePhiShift);
+}
+
+void  AliEMCALRecPointsQaESDSelector::SetMomentum(Double_t p) 
+{ // Jul 9, 2007
+  fPmom = p;
+  //  SetName(Form("%s_p_%f5.1", GetName(), p));
 }
 
+
 AliEMCALFolder*  AliEMCALRecPointsQaESDSelector::CreateEmcalFolder(const Int_t it)
 {
   AliEMCALFolder* newFolder = new AliEMCALFolder(it); // folder for iteration #it   
@@ -401,7 +877,7 @@ AliEMCALFolder*  AliEMCALRecPointsQaESDSelector::CreateEmcalFolder(const Int_t i
 AliEMCALFolder* AliEMCALRecPointsQaESDSelector::GetEmcalOldFolder(const Int_t nsm)
 {
   AliEMCALFolder* folder=0;
-  if(fEmcalPool) folder =  (AliEMCALFolder*)fEmcalPool->FindByName(Form("EMCAL_%2.2i",nsm));
+  if(fEmcalPool) folder =  (AliEMCALFolder*)fEmcalPool->FindObject(Form("EMCAL_%2.2i",nsm));
   return folder;
 }
 
@@ -421,10 +897,15 @@ void AliEMCALRecPointsQaESDSelector::SetEmcalOldFolder(AliEMCALFolder* folder)
 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(geo)         b->Add(geo);
+  if(fCellsInfo)  b->Add(fCellsInfo);
+  //
+  if(fLofHistsPC) b->Add(fLofHistsPC);
+  if(fLofHistsRP) b->Add(fLofHistsRP);
+  if(fLKineVsRP)  b->Add(fLKineVsRP);
+  if(fLShowerProfile) b->Add(fLShowerProfile);
   //  if(u) b->Add(u);
 }
 
@@ -434,10 +915,45 @@ Bool_t AliEMCALRecPointsQaESDSelector::IsFolder() const
   return kFALSE;
 }
 
+void AliEMCALRecPointsQaESDSelector::Save(Int_t ver, const char *optIO)
+{ // Aug 3, 2007
+  TString dir("/home/pavlinov/ALICE/SHISHKEBAB/RF/CALIB/"); // Root directory for saving
+  TString nf=dir;
+  if(GetKeyOptsValue(kPROF)) {
+    nf += "PROF/PROFILE_";
+    nf += ver;
+    nf += ".root";
+    TFile f(nf.Data(), optIO);
+    if(f.IsOpen()) {
+      this->Write();
+      f.Close();
+      printf("<I> Save selectort to file |%s| : optIO %s \n",nf.Data(), optIO);
+    } else {
+      printf("<W> File %s exits ! Increase version : ver %i !\n", nf.Data(), ver);
+    }
+  } else {
+    printf("No PROF option -> no saving now !\n");
+  }
+}
+
+AliEMCALRecPointsQaESDSelector* AliEMCALRecPointsQaESDSelector::ReadSelector(const char* nf)
+{
+  AliEMCALRecPointsQaESDSelector* selector=0;
+
+  TH1::AddDirectory(0);
+  TFile f(nf,"READ");
+  if(f.IsOpen()) {
+    TObject *o = f.Get("AliEMCALRecPointsQaESDSelector");
+    if(o) selector = dynamic_cast<AliEMCALRecPointsQaESDSelector *>(o);
+  }
+  printf("<I> read selector %p : file |%s| \n", selector, nf);
+  return selector;
+}
+
 void AliEMCALRecPointsQaESDSelector::ReadAllEmcalFolders()
 {
   if(fEmcalPool==0) {
-    fEmcalPool = new TObjectSet("PoolOfEMCAL");
+    fEmcalPool = new TFolder("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);
@@ -463,8 +979,9 @@ void AliEMCALRecPointsQaESDSelector::PictVsIterNumber(const Int_t ind, const Int
    
   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));
+  TList* col = (TList*)fEmcalPool->GetListOfFolders();
+  for(Int_t i=0; i<col->GetSize(); i++) { // cycle on EMCAL folders
+    AliEMCALFolder* folder = static_cast<AliEMCALFolder*>(col->At(i));
     if(folder==0) continue;
     it = folder->GetIterationNumber();
 
@@ -503,3 +1020,610 @@ void AliEMCALRecPointsQaESDSelector::PictVsIterNumber(const Int_t ind, const Int
   hout->SetMinimum(0.0);
 
 }
+
+TH1F* AliEMCALRecPointsQaESDSelector::FitHistOfRecPointEnergy(const char *opt)
+{
+  TH1::AddDirectory(0);
+
+  TString OPT(opt);
+  OPT.ToUpper();
+
+  Int_t ind = 4, ind2 = 16;
+  if(GetKeyOptsValue(kIDEAL)) {
+    ind  = 16; // Jul 12, 2007
+    ind2 =  4;
+  } else if(GetKeyOptsValue(kCORR1)) {
+    ind  = 18;
+    ind2 =  4;
+  }
+  TH1F *hold = (TH1F*)fLofHistsRP->At(ind), *h=0;
+  if(hold == 0) return 0;
+  if(hold->GetEntries() <10.) return hold;
+
+  if(OPT.Contains("CLONE")) {
+    TString newName(Form("C_%s",hold->GetName()));
+    h = (TH1F*)hold->Clone(newName.Data());
+    printf(" Clone hist %s -> |%s|%s| \n",hold->GetName(),h->GetName(),h->GetTitle()); 
+  } else {
+    h = hold;
+  }
+  TH1F* h2 = (TH1F*)fLofHistsRP->At(ind2);
+
+  Double_t xmax = h->GetXaxis()->GetXmax(), xmin = 0.4*xmax;;
+  TF1 *g = u::Gausi("fRecPointE", xmin, xmax, h);
+  g->SetLineColor(kRed);
+  gStyle->SetOptFit(111);
+
+  h->Fit(g,"N","", xmin, xmax);
+  printf(" (1) xmin %f : xmax %f \n", xmin, xmax);
+
+  xmin = g->GetParameter(1) - 4.*g->GetParameter(2);
+  xmax = g->GetParameter(1) + 4.*g->GetParameter(2);
+  h->Fit(g,"Q+","", xmin, xmax);
+  u::DrawHist(h2,1, 1, "same",2);
+  printf(" (2) xmin %f : xmax %f \n", xmin, xmax);
+
+  return h;
+}
+
+TCanvas *AliEMCALRecPointsQaESDSelector::Linearity(TList *l, int ifun)
+{ //Jul 10, 2007
+  if(l==0) {
+    printf("<E> AliEMCALRecPointsQaESDSelector::Linearity :TList is zero ! Bye ! \n");
+    return 0;
+  }
+  Double_t p[9]={0.5, 1., 3., 5., 10., 20., 30., 50., 100.}, ep[9];
+  // see macro rHits.C ->defineSampleFraction
+  TH1F* hErecOverEin = new TH1F("hErecOverEin","Ratio E_{rec}/E_{#gamma} vs E_{#gamma}", 101,0.5,101.5);
+  // Fill hist
+  TArrayD invRat(9), eInvRat(9), erec(9),derec(9), residual(9);
+  TArrayD res(9), eres(9); // resolution
+  Int_t np=9;
+  for(Int_t var=0; var<9; var++){
+    TH1F *h = (TH1F*)l->At(var);
+    TF1  *f = (TF1*)h->GetListOfFunctions()->At(0);
+    Double_t  mean = f->GetParameter(1);
+    Double_t emean = f->GetParError(1);
+
+    Int_t bin = hErecOverEin->FindBin(p[var]);
+    hErecOverEin->SetBinContent(bin, mean/p[var]);
+    hErecOverEin->SetBinError(bin,  emean/p[var]); 
+    //
+    invRat[var]  = p[var]/mean;
+    eInvRat[var] = emean/p[var]*invRat[var];
+    erec[var]    = mean;
+    derec[var]   = emean;
+    // Resolution in %
+    res[var]  = 100.*f->GetParameter(2)/p[var];
+    eres[var] = 100.*f->GetParError(2)/p[var];
+    ep[var] = 0.0;
+  }
+
+  TCanvas *c = new TCanvas("Linearity","Linearity", 20,20, 700, 500);
+  gStyle->SetOptStat(0);
+  if(0) { // E(rec) / E(gamma)
+    u::DrawHist(hErecOverEin,2);
+    hErecOverEin->SetMinimum(0.984);
+    hErecOverEin->SetMaximum(1.001);
+  }
+  Int_t markerColor=1;
+  char *fun="", *optFit="";
+  TF1 *f = 0;
+  if(0) {
+    if(ifun==-5) {
+      fun = "fun5"; optFit="R+";
+      f = new TF1(fun,"([0]*(x-7.5)*(x-7.5))*exp([1]+[2]*x+[3]*x*x)+1.", 0.0, 100.);
+      f->FixParameter(0, 1.95380e-05);
+      //      f->SetParameter(0, 1.95380e-05);
+      //      f->SetParameter(1, 1.0);
+    } else if(ifun==-4) {
+      fun = "fun4"; optFit="R+";
+      f = new TF1(fun,"([0]*(x-7.5)*(x-7.5))*([1]+[2]*abs(x)+[3]*(x-[4])*(x-[4]))+1.", 0.0, 100.);
+      f->FixParameter(0, 1.95380e-05);
+      f->SetParameter(4, 50.);
+      //      f = new TF1(fun,"exp(([0]+[1]*x)*(x-[2])*(x-[3]))", 0.0, 100.);
+      //f = new TF1(fun,"([0]+[1]*x)*(x-[2])*(x-[3])", 0.0, 100.);
+      //f->SetParameter(0, 1.);
+      //f->SetParameter(1, 1.);
+
+      //      f->SetParameter(2, 5.);
+      //f->SetParLimits(2, 3.,8.);
+      //f->SetParameter(3, 10.);
+      //f->SetParLimits(3, 9.,15.);
+      //f->SetParameter(2, 2.e-4);
+    } else if(ifun==-3) {
+      fun = "fun3"; optFit="R+";
+      f = new TF1(fun,"[0]+[1]*x+[2]*x*x+[3]*x*x*x", 0.0, 11.);
+    } else if(ifun==-2) {
+      fun = "fun2"; optFit="R+";
+
+      /*
+      f = new TF1(fun,"[0]*(x-7.5)+[1]*(x-7.5)*(x-7.5)+[2]", 0.0, 10.1); 
+      f->SetParameter(0, 5.98727e-04);
+      f->SetParameter(1, 2.12045e-04);
+      f->SetParameter(2, 1.);
+      */
+      f = new TF1(fun,"[0]+[1]*x+[2]*x*x", 9.0, 100.1); 
+    } else if(ifun==-1) {
+      fun = "fun0"; optFit="R+";
+      f = new TF1(fun,"[0]", 0.0, 100.1); 
+    } else if(ifun>=2) {
+      fun = Form("pol%i",ifun); 
+      optFit="+";
+    }
+    TGraphErrors *gr = u::DrawGraphErrors(np, erec.GetArray(),invRat.GetArray(),derec.GetArray(),eInvRat.GetArray(),
+   //    markerColor,21+markerColor,"AP", " Ratio E_{#gamma}/E_{Rec} ", "E_{Rec}  ","   E_{#gamma}/E_{Rec}",
+ markerColor,21+markerColor,"AP", " Ratio E_{#gamma}/E_{Rec}^{corr} ", "E_{Rec}^{corr}  ","   E_{#gamma}/E_{Rec}^{corr}",
+                                         ifun, optFit, fun);
+    gr->GetHistogram()->SetAxisRange(0.0,100.);
+    //    double xmi=0.999, xma=1.017;
+    double xmi=0.995, xma=1.005;
+    gr->GetHistogram()->SetMaximum(xma);
+    gr->GetHistogram()->SetMinimum(xmi);
+    gr->GetHistogram()->SetTitleOffset(1.4,"y");
+    if(ifun==0) {
+      f = new TF1("fres", "AliEMCALHistoUtilities::EnergyCorrectionForGamma_1(x)", 0., 101.); 
+      f->Draw("same");
+    }
+  }
+  TLine *line = new TLine(0.0,1.0, 100.,1.0);
+  line->Draw();
+
+  if(0) {
+    c->Clear();
+    for(int i=0; i<9; i++) {
+      residual[i] =  100.*(invRat[i] - u::GetCorrectionCoefficientForGamma_1(erec[i])); // in percent
+      printf(" erec %f : residual %5.3f \n", erec[i], residual[i]);
+    }
+    markerColor = 2;
+    TGraphErrors *gr2=u::DrawGraphErrors(np, erec.GetArray(),residual.GetArray(),derec.GetArray(),eInvRat.GetArray(),
+                 markerColor,21+markerColor,"AP"," residual in %, rec.point level", "E_{Rec}  ","  residual in %",
+    -1, "", 0);
+    gr2->GetHistogram()->SetAxisRange(0.0,100.);
+    gr2->GetHistogram()->SetMaximum(0.2);
+    gr2->GetHistogram()->SetMinimum(-0.1);
+    line = new TLine(0.0,0.0, 101.,0.0);
+    line->Draw();
+    TLatex *lat = u::Lat("linearity better 0.2% after correction",20., 0.15, 12, 0.06, 1); 
+    if(lat);
+  }
+  if(1) {
+    TString latexName;
+    c->Clear();
+    gStyle->SetOptFit(111);
+    markerColor = 1;
+    ifun        = -11;
+    f = u::GetResolutionFunction("FRES2", latexName);  
+    f->SetNpx(10000);
+    f->SetLineColor(kBlack);
+
+    TGraphErrors *gr3=u::DrawGraphErrors(np, p,res.GetArray(), ep, eres.GetArray(),
+    markerColor,21+markerColor,"AP"," resolution in %, rec.point level","E_{#gamma} "," resolution in %",
+                                        ifun, "+", f->GetName());
+    gr3->GetHistogram()->SetAxisRange(0.0,101.);
+    gr3->GetHistogram()->SetMaximum(14.);
+    gr3->GetHistogram()->SetMinimum(0.0);
+    gr3->SetMarkerSize(1.5);
+
+    TLatex *lat = u::Lat(latexName.Data(),82., 11., 12, 0.06, 1);
+    if(lat);
+    // Exp. data
+    TF1 *fexp = new TF1(*f);
+    fexp->SetName("fexp");
+    fexp->SetParameter(0, 2.168);
+    fexp->SetParameter(1, 8.818);
+    fexp->SetLineWidth(1);
+    fexp->SetLineColor(kRed);
+    fexp->SetLineStyle(1);
+    fexp->Draw("same");
+
+    TLegend *leg = new TLegend(0.21,0.36, 0.68,0.82);
+    leg->AddEntry(gr3,  "MC data", "P");
+    leg->AddEntry(f,    "fit of MC data", "L");
+    TLegendEntry *ent3 = leg->AddEntry(fexp, "#splitline{fit of exp.data}{FNAL, Nov 2005}", "L");
+    ent3->SetTextColor(fexp->GetLineColor());
+    leg->Draw();
+  }
+  c->Update();
+
+  return c;
+}
+
+TCanvas *AliEMCALRecPointsQaESDSelector::DrawKineVsRP(TList *l)
+{ //Jul 25, 2007
+  if(l==0) {
+    printf("<W> AliEMCALRecPointsQaESDSelector::DrawKineVsRP : TList is zero ! \n");
+    return 0;
+  }
+  TCanvas *c = new TCanvas("KineVsRP","KineVsRP", 20,20, 700, 500);
+  gStyle->SetOptStat(1110);
+  c->Divide(2,2);
+
+  c->cd(1);
+  TH1F* h1 = (TH1F*)l->At(10);
+  u::DrawHist(h1,2);
+
+  c->cd(2);
+  TH1F* h2 = (TH1F*)l->At(11);
+  u::DrawHist(h2,2);
+
+  c->cd(3);
+  TH1F* h3 = (TH1F*)l->At(12);
+  u::DrawHist(h3,2);
+  u::DrawHist((TH1F*)l->At(14), 1,kRed,"same");
+  u::DrawHist((TH1F*)l->At(15), 1,kGreen,"same");
+
+  c->cd(4);
+  TH1F* h4 = (TH1F*)l->At(13);
+  u::DrawHist(h4, 2);
+  u::DrawHist((TH1F*)l->At(16), 1,kRed,"same");
+  u::DrawHist((TH1F*)l->At(17), 1,kGreen,"same");
+
+  /*
+  TH1F* h2 = (TH1F*)l->At(11);
+  u::DrawHist(h2,2);
+  */
+
+  c->Update();
+  return c;
+}
+
+TCanvas* AliEMCALRecPointsQaESDSelector::DrawMeanVsLog(TH2F *h2)
+{
+  // h - input histogramm : mean vds log coordinates 
+  if(h2==0) return 0;
+
+  TCanvas *c = new TCanvas("KineVsRP","KineVsRP", 20,20, 700, 500);
+
+  TH1F *hid1 = new TH1F("h1","h1 title ", h2->GetNbinsX(), 
+  h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());  
+
+  gROOT->cd();
+  TH1::AddDirectory(1);
+  TString newName;
+  for(int ix=1; ix<=h2->GetNbinsX();ix++) {
+    newName = "hd_"; newName += ix;
+    TH1D *hd = h2->ProjectionY(newName.Data(),ix,ix);
+    if(hd->Integral()>=4.) {
+      //      lhd->Add(hd);
+      hid1->SetBinContent(ix, hd->GetMean());
+      hid1->SetBinError(ix, hd->GetRMS());
+    }
+  }
+  TF1 *f1 = new TF1("fcg", "0.5*TMath::SinH(x/[0])/TMath::SinH(0.5/[0])", -0.5, 0.5);
+  f1->SetParameter(0,0.13);
+  f1->SetLineColor(kRed);
+  f1->SetLineWidth(1);
+
+  gStyle->SetOptStat(0);
+  gStyle->SetOptFit(111);
+  hid1->Fit(f1,"R+");
+
+  u::DrawHist(hid1,2);
+
+//u::DrawHist(h2,2);
+
+  c->Update();
+  return c;
+}
+
+TCanvas* AliEMCALRecPointsQaESDSelector::DrawPhiEtaAnglesDistribution(const char* gn)
+{ 
+  // Aug 6, 2007
+  // Proper geometry should be defined already 
+  // dTheta distibution has two peaks which coresponding different cells inside module ! 
+
+  TCanvas *c = new TCanvas("Geometry","Geometry", 20,20, 700, 500);
+  c->Divide(2,2);
+
+  if(geo==0) geo = AliEMCALGeometry::GetInstance(gn);
+
+  gROOT->cd();
+  TH1::AddDirectory(1);
+  TH1F *hDtheta  = new TH1F("hDtheta","#Delta#theta in one SM", 60, -2.0, +1.0); // in degree
+  TH2F *hDtheta2 = new TH2F("hDtheta2","#Delta#theta vs iEta of cell", 48, -0.5, 47.5, 60,-2.0,+1.0);
+
+  TH1F *hDphi  = new TH1F("hDphi","#Delta#ph in one SM", 2000, -10.0, +10.0); // in degree
+
+  TVector3 vg3;
+  AliEMCALCellInfo* t = GetCellsInfo();
+  Double_t thetaCell=0., thetaModule=0.;
+  Double_t phiCell=0., dphi=0.;
+
+  for(int absid=0; absid<12*24*4; absid++){
+    cellInfo *r = t->GetTable(absid);
+    geo->GetGlobal(absid, vg3);
+
+    thetaCell   = vg3.Theta()*TMath::RadToDeg();
+    thetaModule = 90. - 1.5*r->iEtam;
+    hDtheta->Fill(thetaCell - thetaModule);
+    hDtheta2->Fill(double(r->iEta), thetaCell - thetaModule);
+
+    phiCell   = vg3.Phi()*TMath::RadToDeg();
+    dphi      = phiCell - 90.;
+    hDphi->Fill(dphi);
+  }
+
+  c->cd(1);
+  u::DrawHist(hDtheta,2);
+  c->cd(2);
+  u::DrawHist(hDtheta2,2);
+  c->cd(3);
+  u::DrawHist(hDphi,2);
+
+  c->Update();
+  return c;
+}
+
+TCanvas* AliEMCALRecPointsQaESDSelector::DrawDeffVsEnergy()
+{
+  // Aug 31, 2007 - obsolete
+  Double_t p[]={0.5, 1.0, 3., 5., 10., 20.,30.,50.,100.};
+  Double_t ep[]={0., 0., 0., 0., 0., 0., 0., 0., 0.};
+  // 2 pars
+  Double_t  deff[]={9.07515, 10.0835, 11.2032, 11.4229, 12.3578, 13.0332, 13.3281, 13.7910, 14.3319};
+  Double_t edeff[]={2.69645e-03, 3.52180e-03, 4.19236e-02, 1.21201e-01, 1.58886e-01, 3.96680e-01, 3.29985e-02, 1.17113e-01, 5.22763e-01};
+  //
+  Double_t w0[]={3.44679e+00, 3.82714e+00, 4.15035e+0, 4.36650e+00, 4.51511e+00, 4.65590, 4.63289e+00, 4.66568, 4.68125};
+  Double_t ew0[]={7.58982e-01, 1.26420e-02, 2.36129e-10, 1.21201e-01, 2.12999e-01, 7.95650e-02, 6.15307e-03, 1.88803e-01, 5.18022e-05};
+  int np=sizeof(p)/sizeof(Double_t);
+  printf("<I> AliEMCALRecPointsQaESDSelector::DrawDeffVsEnergy() | np %i \n", np);
+
+  TCanvas *c = new TCanvas("Deff","Deff", 20,20, 700, 500);
+  c->Divide(2,1);
+
+  Int_t markerColor=1;
+  c->cd(1);
+  TF1 *fdeff= new TF1("fdeff","[0]+[1]*log(x)",0.4, 100.4);
+  if(fdeff);
+  TGraphErrors *gr = u::DrawGraphErrors(np, p,deff, ep, edeff,
+    markerColor,21+markerColor,"AP"," D_{eff} vs E_{#gamma} ","E_{#gamma}         "," D_{eff} in cm ",
+                                       -1, "", 0);
+  //                                   -1, "+", fdeff->GetName());
+  gr->GetHistogram()->SetMaximum(15.);
+  gPad->SetLogx(1);
+
+  c->cd(2);
+  TGraphErrors *gr2 = u::DrawGraphErrors(np, p,w0, ep, ew0,
+    markerColor,21+markerColor,"AP"," w_{0} vs E_{#gamma} ","E_{#gamma}         "," w_{0}    ",
+    -1, "", 0);
+
+  gr2->GetHistogram()->SetMaximum(5.);
+  gPad->SetLogx(1); 
+
+  c->Update();
+  return c;
+}
+
+TCanvas* AliEMCALRecPointsQaESDSelector::DrawDeffVsEnergy2(const char *opt)
+{
+  // Aug 28, 2008 - read pars and pars errors from files
+  Double_t p[]={0.5, 1.0, 3., 5., 10., 20.,30.,50.,100.};
+  Double_t ep[]={0., 0., 0., 0., 0., 0., 0., 0., 0.};
+  // 2 pars
+  Double_t deff[9], edeff[9], w0[9], ew0[9]; // max size now 
+  TString OPT(opt);
+
+  int np = sizeof(p)/sizeof(Double_t);
+  printf("<I> AliEMCALRecPointsQaESDSelector::DrawDeffVsEnergy2() | np %i \n", np);
+  ReadParsDeffAndW0("/data/r22b/ALICE/CALIB/FIT/", deff, edeff, w0, ew0, 1);
+
+  TCanvas *c = new TCanvas("Deff","Deff", 20,20, 700, 500);
+  c->Divide(2,1);
+
+  TF1 *fdeff = 0, *fw0 = 0;
+  TString optFit(""), funName("");
+  if(OPT.Contains("fit1")) {
+    fdeff= new TF1("fdeff","[0]+[1]*log(x)",0.1, 101.);
+    fdeff->SetLineColor(kRed);
+    fdeff->SetLineWidth(1);
+
+    /* good description - "[0]/(1.+exp([1]*x))"
+   1  p0           4.82208e+00   9.93617e-04  -0.00000e+00   7.16648e-06
+   2  p1          -7.79655e-01   4.07500e-03  -0.00000e+00   2.01009e-03
+      better description - "[0]/(1.+exp([1]*(x-[2])))" (like the Woods-Saxon potential)
+   1  p0           4.83713e+00   1.00437e-03  -6.98984e-10   4.90371e-07
+   2  p1          -2.77970e-01   3.35587e-03  -1.79239e-09   2.67312e-07
+   3  p2           4.41116e+00   8.72191e-02   1.82791e-04   1.55643e-05
+     */
+    fw0= new TF1("fw0","[0]/(1.+exp([1]*(x+[2])))",0.1, 101.);
+    fw0->SetLineColor(kRed);
+    fw0->SetLineWidth(1);
+    fw0->SetParameter(0, 4.8);
+    fw0->SetParameter(1, -2.77970e-01);
+    fw0->SetParameter(2,  4.41116);
+    optFit = "+";
+  }
+  if(fdeff) funName =  fdeff->GetName();
+
+  Int_t markerColor=1;
+  c->cd(1);
+  gStyle->SetOptFit(111);
+  TGraphErrors *gr = u::DrawGraphErrors(np, p,deff, ep, edeff,
+    markerColor,21+markerColor,"AP"," D_{eff} vs E_{#gamma} ","E_{#gamma}         "," D_{eff} in cm ",
+                                       -1, optFit.Data(), funName.Data());
+  gr->GetHistogram()->SetMaximum(15.);
+  gPad->SetLogx(1);
+  TLegend *leg1 = new TLegend(0.12,0.76, 0.70,0.90);
+  TLegendEntry *le1 = leg1->AddEntry(fdeff, Form("%s",fdeff->GetTitle()), "lp");
+  //TLegendEntry *le1 = leg1->AddEntry(fdeff, Form("%4.2f+%4.2f*log(E_{#gamma})",
+  //fdeff->GetParameter(0),fdeff->GetParameter(1)), "lp");
+  le1->SetTextColor(fdeff->GetLineColor());
+  leg1->Draw();
+
+  c->cd(2);
+  gStyle->SetOptFit(111);
+  if(fw0) funName =  fw0->GetName();
+  TGraphErrors *gr2 = u::DrawGraphErrors(np, p,w0, ep, ew0,
+    markerColor,21+markerColor,"AP"," w_{0} vs E_{#gamma} ","E_{#gamma}         "," w_{0}    ",
+                                       -1, optFit.Data(), funName.Data());
+
+  gr2->GetHistogram()->SetMaximum(5.);
+
+  TLegend *leg2 = new TLegend(0.17,0.6, 0.99,0.72);
+  TLegendEntry *le2 = leg2->AddEntry(fw0, Form("%s",fw0->GetTitle()), "lp");
+  //TLegendEntry *le2 = leg2->AddEntry(fw0, Form("#frac{%4.2f}{1.+exp(%4.2f*(x+%4.2f)}",
+  //fw0->GetParameter(0),fw0->GetParameter(1),fw0->GetParameter(2)), "lp");
+  le2->SetTextColor(fw0->GetLineColor());
+  leg2->Draw();
+  //gPad->SetLogx(1); 
+
+  c->Update();
+  return c;
+}
+
+void AliEMCALRecPointsQaESDSelector::ReadParsDeffAndW0
+(const char *dirName, double *deff, double *edeff, double *w0, double *ew0, const Int_t pri)
+{
+  int strategy = 0, itmp=0;
+  char line[100];
+  for(int var=11; var<=19; var++){
+    int ind = var -11;
+    ifstream fin;
+    TString fname = Form("%s/fitVar%iStrategy%i.txt", dirName, var, strategy);
+    //printf(" open file %s \n", fname.Data());
+    fin.open(fname.Data());
+   
+    for(int i=1; i<=2; i++) { // skip to lines
+      fin.getline(line,100).eof();
+      if(pri>=2) printf("%s \n", line);
+    }
+
+    fin >> itmp >> deff[ind] >> edeff[ind];
+    fin >> itmp >> w0[ind]   >> ew0[ind];
+    if(pri>=1) 
+    printf(" %i def %f+/-%f : def %f+/-%f\n", ind, deff[ind],edeff[ind],w0[ind],ew0[ind]); 
+    fin.getline(line,100).eof(); // skip last line
+    fin.close();
+  }
+}
+
+TCanvas* AliEMCALRecPointsQaESDSelector::DrawSpaceResolution()
+{
+  // Sep 4, 2007;
+  // Space resolution vs optimised space resolution
+  Double_t p[]={0.5, 1.0, 3., 5., 10., 20.,30.,50.,100.};
+  Double_t ep[]={0., 0., 0., 0., 0., 0., 0., 0., 0.};
+  Double_t  spAng[]={0.1877, 0.1351, 0.08144, 0.06331, 0.05384, 0.04876, 0.04706, 0.04656, 0.04726};
+  Double_t rmsSpAng[]={0.1033, 0.07685, 0.05235, 0.03992, 0.04012, 0.04257, 0.03472, 0.02814, 0.02784};
+  Double_t  spAngOpt[]={0.1733, 0.1311, 0.07961, 0.06401, 0.05347, 0.04618, 0.04288, 0.04, 0.03802};
+  Double_t rmsSpAngOpt[]={0.09476, 0.07472, 0.05011, 0.04242, 0.04075, 0.04304, 0.03545, 0.02744, 0.02593};
+
+  int np=sizeof(p)/sizeof(Double_t);
+  printf("<I> AliEMCALRecPointsQaESDSelector::DrawSpaceResolution() | np %i \n", np);
+
+  Double_t* eSpAng    = new Double_t[np];
+  Double_t* eSpAngOpt = new Double_t[np];
+  Double_t C=TMath::Sqrt(8000.), C2 = 1000.*TMath::DegToRad();
+  for(int i=0; i<np; i++){
+    spAng[i]       *= C2;
+    rmsSpAng[i]    *= C2;
+    spAngOpt[i]    *= C2;
+    rmsSpAngOpt[i] *= C2;
+
+    eSpAng[i]    = spAng[i]/C;
+    eSpAngOpt[i] = spAngOpt[i]/C;
+  }
+
+  TCanvas *c = new TCanvas("Deff","Deff", 20,20, 700, 500);
+  //c->Divide(2,1);
+
+  c->cd(1);
+  gStyle->SetOptFit(111);
+
+  TGraphErrors *gr = u::DrawGraphErrors(np, p,spAng, ep, eSpAng,
+  kBlack, 21,"AP","Angle resolution (mrad) vs E_{#gamma} ","E_{#gamma}       "," anlgle resolution (mrad) ",
+                                       -1, "", 0);
+  gr->GetHistogram()->SetMaximum(4.);
+  gr->GetHistogram()->SetMinimum(0.6);
+  gPad->SetLogx(1);
+  gPad->SetLogy(1);
+
+  TF1 *fang = new TF1("fang","[0]+[1]/sqrt(x)",0.1, 101.);
+  fang->SetLineColor(kRed);
+  fang->SetLineWidth(1);
+
+  TGraphErrors *gr2 = u::DrawGraphErrors(np, p,spAngOpt, ep, eSpAngOpt,
+  kRed,22,"P"," space vs E_{#gamma} ","E_{#gamma}         "," anlgle resolution (mrad) ",
+                                        -1, "+", fang->GetName());
+  TLegend *leg1 = new TLegend(0.40,0.57, 0.92,0.87);
+  TLegendEntry *le1 = leg1->AddEntry(gr, "Initial angle resolution", "p");
+  le1->SetTextColor(gr->GetLineColor());
+
+  TLegendEntry *le2 = leg1->AddEntry(gr2, "Optimized angle resolution", "p");
+  le2->SetTextColor(gr2->GetLineColor());
+
+  TLegendEntry *le3 = leg1->AddEntry(fang, 
+  Form("%3.2f+#frac{%4.2f}{#sqrt{E}}",fang->GetParameter(0),fang->GetParameter(1)), "lp");
+  le3->SetTextColor(fang->GetLineColor());
+
+  leg1->Draw();
+
+  c->Update();
+
+  return c;
+}
+
+void AliEMCALRecPointsQaESDSelector::ResetAllListOfHists()
+{
+  u::ResetListOfHists(fLofHistsPC);
+  u::ResetListOfHists(fLofHistsRP);
+  u::ResetListOfHists(fLKineVsRP);
+  u::ResetListOfHists(fLShowerProfile);
+}
+
+void AliEMCALRecPointsQaESDSelector::ReloadChain(Long64_t entry)
+{ // unused now
+  if(fChain) {
+    fChain->LoadTree(entry);
+  }
+}
+
+void AliEMCALRecPointsQaESDSelector::GetInitialParsForFit(const Int_t var, Double_t &deff, Double_t &w0, Double_t &phislope, const int phiCase)
+{
+  //  int phiCase=0; // 0 or 1
+  phislope = 0.001;
+  if(var==11)         { // stay here
+    deff = 9.07515;
+    w0   = 3.44679;
+  }  else if(var==12) {
+    deff = 1.00835e+01;
+    w0   = 3.82714e+00;
+  }  else if(var==13) {
+    deff = 1.12032e+01;
+    w0   = 4.15035e+00;
+  }  else if(var==14) {
+    deff = 1.14229e+01;
+    w0   = 4.36650;
+  }  else if(var==15) {
+    deff = 1.21361e+01; 
+    w0   = 4.72875;
+  }  else if(var==16) {
+    deff = 1.28096e+01;
+    w0   = 4.81753e+00;
+  }  else if(var==17) {
+    deff = 1.33281e+01; 
+    w0   = 4.63289e+00;
+  }  else if(var==18) {
+    deff = 1.37910e+01; 
+    w0   = 4.66568;
+  }  else if(var==19) {// Aug 15, 2007
+    switch (phiCase) {
+    case 0: // phislope was defined than deff and w0 were fixed
+      deff = 1.43319e+01;
+      w0   = 4.69279;
+      phislope = 2.41419e-04;
+      break;
+    case 1: // all parameters free
+      deff = 1.46327e+01;
+      w0   = 4.68125;
+      phislope = 8.95559e-04;
+      break;
+    }
+  } else { 
+    printf("<E> var % i -> no definition ! \n", var);
+    assert(0);
+  }
+  printf("<I> AliEMCALRecPointsQaESDSelector::GetInitialParForFit()\n deff %f w0 %f phislope %f \n", 
+        deff,w0, phislope);
+}
index 9837810..a447ee5 100644 (file)
 
 #include "AliSelector.h"
 
+#include <TObjArray.h>
+
 class AliEMCALFolder;
+class AliRunLoader;
+class AliEMCALRecPoint;
+class AliEMCALCellInfo;
+class cellInfo;
 
 class TList;
+class TCanvas;
+class TH1F;
+class TH2F;
 class TBrowser;
 class TChain;
-class TObjectSet;
+class TFolder;
+class TArrayI;
 
 class AliEMCALRecPointsQaESDSelector :  public AliSelector {
   public:
     AliEMCALRecPointsQaESDSelector();
     virtual ~AliEMCALRecPointsQaESDSelector();
 
-    virtual void    Begin(TTree* tree);
-    virtual void    Init(TTree* tree);
-    virtual void    SlaveBegin(TTree *tree);
+    virtual void    Begin(TTree*);
+    virtual void    SlaveBegin(TTree* tree);
+    virtual void    Init(TTree *tree);
+    virtual Bool_t  Notify();
     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");
+  static TList *DefineHistsOfRP(const char *name="RP",  Double_t p=110.0, Int_t keyOpt=0);
+  static TList *DefineHistsOfKineVsRP(const char *name="KineVsRP",  Double_t p=110.0, Int_t keyOpt=0);
+  static TList *DefineHistsForShowerProfile(const char *name="ProfY", Double_t p=1.);
+  static void   FillHistsOfKineVsRP(TList *l, AliRunLoader* RL, TClonesArray &lvM);
+  static void   FillHistsForShowerProfile(TList *l, AliEMCALRecPoint *rp, AliEMCALCellInfo* t);
+  static void   EvalLocalPhiPosition(const Double_t wlog, const AliEMCALRecPoint *rp, const AliEMCALCellInfo* t, Double_t &xcog, Int_t &phiSize, cellInfo &rMax); 
+  //
+  TList *GetListKineVsRP()      {return fLKineVsRP;}
+  TList *GetListShowerProfile() {return fLShowerProfile;}
+  //  static TList *DefineHistsOfTowers(const char *name="towers");
   // 
-  void FitEffMassHist(); //*MENU*  
+  void FitEffMassHist(); // *MENU*  
   void PrintInfo();      // *MENU*
   //
   void    SetChain(TChain *chain)  {fChain = chain;}
   TChain* GetChain()               {return fChain;}
+  void    SetMomentum(Double_t p);
+  Double_t GetMomentum() {return fPmom;}
 
   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);
-  // 
+  AliEMCALCellInfo* GetCellsInfo() {return fCellsInfo;}
+  //
+  void      SetStringOfRunOpts(const char *st) {fRunOpts = st;}
+  TObjArray GetOptsArray() const {return fArrOpts;}
+  Int_t     GetKeyOptsValue(Int_t key);  // *MENU*
+  void      CheckRunOpts();
+  //
   virtual void Browse(TBrowser* b);
   virtual Bool_t  IsFolder() const;
+  //
+  void    Save(Int_t ver=0, const char *optIO="NEW");   // *MENU*
+  static  AliEMCALRecPointsQaESDSelector* ReadSelector(const char* nf = "/home/pavlinov/ALICE/SHISHKEBAB/RF/CALIB/PROF/PROFILE_0.root");
+  
 
   //
   //// Pictures staf - Jun 26, 2007
   //
   void ReadAllEmcalFolders();
   void PictVsIterNumber(const Int_t ind=0, const Int_t nsm=0);
-
+  // Gamma
+  TH1F* FitHistOfRecPointEnergy(const char *opt="CLONE");
+  static TCanvas *Linearity(TList *l, Int_t ifun=3);
+  static TCanvas *DrawKineVsRP(TList *l);
+  // Profile
+  static TCanvas *DrawMeanVsLog(TH2F *h2);
+  // Geometry staff
+  TCanvas *DrawPhiEtaAnglesDistribution(const char *gn="SHISH_TRD1_CURRENT_2X2"); // *MENU*
+  // Geometry constants 
+  TCanvas *DrawDeffVsEnergy();  // *MENU*
+  TCanvas *DrawDeffVsEnergy2(const char *opt="fit1"); // *MENU*
+  void     ReadParsDeffAndW0(const char *dirName="/data/r22b/ALICE/CALIB/FIT/",
+          double *deff=0, double *edeff=0, double *w0=0, double *ew0=0, const Int_t pri=0);
+  TCanvas *DrawSpaceResolution();
  protected:
+  Double_t fPmom; // positive if defined
   //
   TChain* fChain; //! chain if ESD files
-  TList* fLofHistsPC; //! list of histograms of pseudo clusters 
-  TList* fLofHistsRP; //! list of histograms of rec.points 
+  TList* fLofHistsPC; // list of histograms of pseudo clusters 
+  TList* fLofHistsRP; // list of histograms of rec.points 
+  TList* fLKineVsRP;  // list of histograms kinematics vs rec.points 
+  TList* fLShowerProfile;  // list of histograms for shower profile business
+  //
+  AliEMCALCellInfo *fCellsInfo; //
+  TFolder*         fEmcalPool;  //
+  //
+  // Options - Jul 10, 2007
   //
-  TObjectSet*     fEmcalPool;
-  AliEMCALFolder* fEMCAL;    //! current  EMCAL object
-  AliEMCALFolder* fEMCALOld; //! previous EMCAL object
+  TString   fRunOpts;        // String of running options
+  TObjArray fArrOpts;        // Array of options 
+  // Options keys
+  TArrayI  *fKeyOpts;
+  // Static parameters
+ public:
+  static AliEMCALFolder* GetEmcalFolder() {return fEMCAL;}
+  static AliEMCALFolder* GetEmcalOldFolder() {return fEMCALOld;}
+  static AliEMCALFolder*  fEMCAL;      // current  EMCAL object
+  static AliEMCALFolder*  fEMCALOld;   // previous EMCAL object
+  //
+  static Double_t fDistEff;  // effective depth of electromagnetic shower
+  static Double_t fW0;       // 
+  static Double_t fSlopePhiShift; // phi shift of cluster = fSlopePhiShift * phi
+  static void SetFitParameters(Double_t deff, Double_t w0, Double_t slope) 
+  {
+    fDistEff = deff; fW0 = w0; fSlopePhiShift = slope;
+  }
+  static void GetFitParameters(Double_t &deff, Double_t &w0, Double_t &slope)
+  {
+    deff = fDistEff; w0 = fW0; slope = fSlopePhiShift;
+  } 
+  void ResetAllListOfHists();
+  void ReloadChain(Long64_t entry=0);
+  void GetInitialParsForFit(const Int_t var, Double_t &deff, Double_t &w0, Double_t &phislope, const int phiCase=0);
 
  private:
   AliEMCALRecPointsQaESDSelector(const AliEMCALRecPointsQaESDSelector&);
@@ -70,5 +142,4 @@ class AliEMCALRecPointsQaESDSelector :  public AliSelector {
 
   ClassDef(AliEMCALRecPointsQaESDSelector, 1);
 };
-
 #endif
index ba73b60..bf698e4 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Id$ */
+/* 
+$Log$ 
+*/
 
 //_________________________________________________________________________
 // Top EMCAL folder which will keep all information about EMCAL itself,
 // super Modules (SM), modules, towers, set of hists and so on.
+//  TObjectSet -> TFolder; Sep 6, 2007
 //
 //*-- Author: Aleksei Pavlinov (WSU, Detroit, USA) 
 
@@ -39,16 +42,15 @@ typedef  AliEMCALHistoUtilities u;
 
 ClassImp(AliEMCALSuperModule)
 
-AliEMCALSuperModule::AliEMCALSuperModule() : TObjectSet(),
-fSMNumber(0)
+AliEMCALSuperModule::AliEMCALSuperModule() : TFolder()
+, fParent(0),fLh(0),fSMNumber(0)
 {
 }
 
 AliEMCALSuperModule::AliEMCALSuperModule(const Int_t m, const char* title) : 
-TObjectSet(Form("SM%2.2i",m)),
-fSMNumber(m)
-{
-  SetTitle(title);
+TFolder(Form("SM%2.2i",m), title)
+, fParent(0),fLh(0),fSMNumber(m)
+{ 
 } 
 
 AliEMCALSuperModule::~AliEMCALSuperModule()
@@ -58,15 +60,18 @@ AliEMCALSuperModule::~AliEMCALSuperModule()
 
 void AliEMCALSuperModule::Init()
 {
-  if(GetHists()==0) this->AddObject(BookHists(), kTRUE);
+  if(GetHists()==0) {
+    fLh = BookHists();
+    Add(fLh);
+  }
 }
 
 void  AliEMCALSuperModule::AddCellToEtaRow(AliEMCALCell *cell, const Int_t etaRow)
 {
-  static TDataSet *set;
-  set = FindByName(Form("ETA%2.2i",etaRow)); 
+  static TFolder *set;
+  set = dynamic_cast<TFolder*>(FindObject(Form("ETA%2.2i",etaRow))); 
   if(set==0) {
-    set = new  TDataSet(Form("ETA%2.2i",etaRow));
+    set = new  TFolder(Form("ETA%2.2i",etaRow),"eta row");
     Add(set);
   }
   set->Add(cell);
@@ -76,11 +81,14 @@ void AliEMCALSuperModule::FitForAllCells()
 {
   Int_t ncells=0;
   for(int eta=0; eta<48; eta++) { // eta row
-    TDataSet *setEta = FindByName(Form("ETA%2.2i",eta));
+    TFolder *setEta = dynamic_cast<TFolder*>(FindObject(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);
+      TList* l = (TList*)setEta->GetListOfFolders();
+      printf(" eta %i : %s : cells %i \n", eta, setEta->GetName(), l->GetSize());
+      for(int phi=0; phi<l->GetSize(); phi++) { // cycle on cells (phi directions)
+        AliEMCALCell* cell = dynamic_cast<AliEMCALCell*>(l->At(phi));
+        if(cell == 0) continue; 
+
         cell->FitEffMassHist();
 
         u::FillH1(GetHists(), 1, cell->GetCcIn()*1.e+3);
@@ -124,9 +132,12 @@ void AliEMCALSuperModule::DrawCC(int iopt)
   TH1 *hCCIn  = (TH1*)GetHists()->At(1);
   TH1 *hCCOut = (TH1*)GetHists()->At(2);
 
+  if(hCCIn == 0)              return;
+  if(hCCIn->GetEntries()<10.) return;
+
   hCCIn->SetStats(kFALSE);
   hCCOut->SetStats(kFALSE);
-  hCCOut->SetTitle("CC in and out");
+  hCCOut->SetTitle("CC in and out; Iter 1; Jul 26; All Statistics");
   hCCOut->SetXTitle("cc in MeV");
   hCCOut->SetYTitle("  N  ");
 
@@ -137,7 +148,10 @@ void AliEMCALSuperModule::DrawCC(int iopt)
   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());
+
+  if(hCCOut->GetEntries()>10.) 
   leg->AddEntry(hCCOut, Form("output cc : %6.2f #pm %6.2f", hCCOut->GetMean(),hCCOut->GetRMS()), "L");
+
   leg->Draw();
 
   if(c) c->Update();
@@ -146,9 +160,13 @@ void AliEMCALSuperModule::DrawCC(int iopt)
 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();
+  TList* l = (TList*)GetListOfFolders();
+  for(int eta=0; eta<l->GetSize(); eta++) { // cycle on eta row
+    TFolder *setEta = dynamic_cast<TFolder*>(l->At(eta));
+    if(setEta==0) continue;
+
+    TList* le = (TList*)setEta->GetListOfFolders();
+    ncells += le->GetSize();
   }
   return ncells;
 }
index 0b7c66f..226f4f5 100644 (file)
@@ -8,16 +8,18 @@
 //_________________________________________________________________________
 //  Emcal Super Module     
 //                  
-//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA) 
+//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
+//  Super Module folder
+//  Initial version was created with TDataSet staf
+//  TObjectSet -> TFolder; Sep 6, 2007
 
-// --- ROOT system ---
 
-#include <TObjectSet.h>
+#include <TFolder.h>
 
 class TList;
 class AliEMCALCell;
 
-class AliEMCALSuperModule : public TObjectSet {
+class AliEMCALSuperModule : public TFolder {
 
  public:
   
@@ -28,7 +30,9 @@ class AliEMCALSuperModule : public TObjectSet {
 
   void Init();
   void   AddCellToEtaRow(AliEMCALCell *cell, const Int_t etaRow);
-  TList* GetHists() {return (TList*)fObj;}
+  TList*   GetHists()  {return fLh;}
+  TObject* GetParent() {return fParent;}
+  void     SetParent(TObject *parent) {fParent=parent;}
   // MENU
   void FitForAllCells(); //*MENU*  
   void FitEffMassHist(); //*MENU*  
@@ -40,9 +44,11 @@ class AliEMCALSuperModule : public TObjectSet {
  protected:
   TList* BookHists();
   //
-  Int_t  fSMNumber;
+  TObject* fParent; // parent
+  TList*   fLh;     // List of hists
+  Int_t    fSMNumber;
 
-  ClassDef(AliEMCALSuperModule,1) // EMCAL SuperModule
+  ClassDef(AliEMCALSuperModule,2) // EMCAL SuperModule
     
 };