New method SetNCarlos in order to set the number of CARLOS chips connected to the...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALFolder.cxx
CommitLineData
16d3c94d 1/**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18//_________________________________________________________________________
19// Top EMCAL folder which will keep all information about EMCAL itself,
20// super Modules (SM), modules, towers, set of hists and so on.
21//
22//*-- Author: Aleksei Pavlinov (WSU, Detroit, USA)
23
24#include "AliEMCALFolder.h"
25#include "AliEMCALHistoUtilities.h"
26#include "AliEMCALGeometry.h"
27#include "AliEMCALSuperModule.h"
28#include "AliEMCALCell.h"
29#include "AliESDCaloCluster.h"
30
31#include "AliRun.h"
32
33#include "AliEMCALCalibData.h"
34#include "AliCDBMetaData.h"
35#include "AliCDBId.h"
36#include "AliCDBEntry.h"
37#include "AliCDBManager.h"
38#include "AliCDBStorage.h"
39
40#include "AliEMCALCalibCoefs.h"
41#include "AliEMCALDigit.h"
42#include "AliEMCALRecPoint.h"
43
44#include "AliEMCALPi0SelectionParam.h"
45
46#include <cassert>
47
48#include <TROOT.h>
49#include <TStyle.h>
50#include <TList.h>
51#include <TH1.h>
52#include <TF1.h>
53#include <TFile.h>
54#include <TFileIter.h>
55#include <TCanvas.h>
56#include <TClonesArray.h>
57#include <TKey.h>
58#include <TNtuple.h>
59#include <TLegend.h>
60#include <TArrayS.h>
61
62 const TString AliEMCALFolder::fgkBaseFolderName("EMCAL");
63 const TString AliEMCALFolder::fgkCCFirstName("CCFIRST");
64 const TString AliEMCALFolder::fgkCCinName("CCIN");
65 const TString AliEMCALFolder::fgkCCoutName("CCOUT");
66
67typedef AliEMCALHistoUtilities u;
68
69ClassImp(AliEMCALFolder)
70
71//AliEMCALGeometry* AliEMCALFolder::fGeometry = 0;
72
73 TList *lObj = 0;
74
75AliEMCALFolder::AliEMCALFolder() :
76 TObjectSet(),
77 fCounter(0), fGeometry(0), fNumOfCell(0), fLhists(0), fLofCells(0),fPi0SelPar(0),fCalibData(0),
78 fCellNtuple(0)
79{
80}
81
82AliEMCALFolder::AliEMCALFolder(const char* name, const char* title, Bool_t putToBrowser) :
83 TObjectSet(name),
84 fCounter(-1), fGeometry(0), fNumOfCell(0), fLhists(0), fLofCells(0),fPi0SelPar(0),fCalibData(0),
85 fCellNtuple(0)
86{
87 SetTitle(title);
88 Init(putToBrowser);
89}
90
91AliEMCALFolder::AliEMCALFolder(const Int_t it, const char* title, Bool_t putToBrowser) :
92 TObjectSet(Form("%s_%2.2i", AliEMCALFolder::fgkBaseFolderName.Data(),it)),
93 fCounter(it), fGeometry(0), fNumOfCell(0), fLhists(0), fLofCells(0),fPi0SelPar(0),fCalibData(0),
94 fCellNtuple(0)
95{
96 SetTitle(title);
97 Init(putToBrowser);
98}
99
100AliEMCALFolder::~AliEMCALFolder()
101{
102 // dtor
103}
104
105void AliEMCALFolder::Init(Bool_t putToBrowser)
106{
107 lObj = new TList;
108 lObj->SetName("Objects"); // name is good ?
109 this->AddObject((TObject*)lObj, kTRUE);
110 // Get default geometry - "SHISH_77_TRD1_2X2_FINAL_110DEG"; May 29, 2007
111 fGeometry = AliEMCALGeometry::GetInstance(); // should be define before
112 lObj->Add(fGeometry);
113
114 // Initial cc with decalibration
115 char* dbString = "local:///data/r22b/ALICE/PROD/CALIBRATION_May_2007/PI0/PDSF/10GEV/DECALIB/DeCalibDB"; // hard coded now
116 fCalibData = (AliEMCALCalibData*)
117 (AliCDBManager::Instance()->GetStorage(dbString)
118 ->Get("EMCAL/Calib/Data",gAlice->GetRunNumber())->GetObject());
119 lObj->Add(fCalibData);
120 // Jun 13, 2007
121 Add(AliEMCALCalibCoefs::GetCalibTableFromDb(fgkCCFirstName.Data()));
122 if(GetIterationNumber()==1) {
123 Add(AliEMCALCalibCoefs::GetCalibTableFromDb()); // first iteration only
124 }
125 // Selection Parameter
126 fPi0SelPar = AliEMCALPi0SelectionParam::Set1();
127 this->Add(fPi0SelPar);
128 //
129 fLhists = BookHists();
130 lObj->Add(fLhists);
131
132 // dimension should be get from geometry - 4*12*24*11);
133 fNumOfCell = fGeometry->GetNCells();
134 fLofCells = new AliEMCALCell*[fNumOfCell];
135 for(int i=0; i<fNumOfCell; i++) fLofCells[i] = 0;
136
137 printf("<I> Create AliEMCALFolder : it %i : name %s\n ", fCounter, GetName());
138
139 if(putToBrowser) gROOT->GetListOfBrowsables()->Add(this); // for testing purpuse
140}
141
142AliEMCALSuperModule* AliEMCALFolder::GetSuperModule(const Int_t nm)
143{
144 TDataSet *set = 0;
145 AliEMCALSuperModule *sm = 0;
146
147 set = FindByName(Form("SM%2.2i",nm));
148 if(set) sm = (AliEMCALSuperModule*)set;
149
150 return sm;
151}
152
153
154AliEMCALCell* AliEMCALFolder::GetCell(const Int_t absId)
155{ // May 30, 2007
156 if(absId<0 || absId >= fNumOfCell) return 0;
157 else return fLofCells[absId];
158}
159
160void AliEMCALFolder::SetCell(AliEMCALCell *cell, const Int_t absId)
161{
162 if(absId>=0 && absId < fNumOfCell) {
163 fLofCells[absId] = cell;
164 }
165}
166
167pi0SelectionParam* AliEMCALFolder::GetPi0SelectionParRow(Int_t nrow)
168{
169 pi0SelectionParam* r=0;
170 if(fPi0SelPar) {
171 r = fPi0SelPar->GetTable(nrow);
172 }
173 return r;
174}
175
176void AliEMCALFolder::FillPi0Candidate(const Double_t mgg, AliESDCaloCluster* cl1, AliESDCaloCluster* cl2)
177{
178 static Int_t absIdMax, nm1, nm2;
179
180 u::FillH1(fLhists, 0, 1.); // number entries
181 u::FillH1(fLhists, 1, mgg);
182
183 nm1 = GetSMNumber(cl1);
184 nm2 = GetSMNumber(cl2);
185
186 if(nm1==-1 || nm2==-1) assert(0);
187
188 if(nm1 != nm2) return; // Both cluster should be in the same SM
189
190 AliESDCaloCluster* cl = cl1;
191 if(cl1->E() < cl2->E()) cl = cl2; // Get cluster with highest energy
192
193 const Int_t nDigits = cl->GetNumberOfDigits();
194 // const UShort_t* adc = cl->GetDigitAmplitude(); // for future
195 const Short_t* absId = cl->GetDigitIndex()->GetArray();
196
197 int indMax = 0;
198 double emax = cl->GetDigitAmplitude()->At(0);
199 if(nDigits > 1) {
200 for(int i=1; i<nDigits; i++) {
201 if(emax < cl->GetDigitAmplitude()->At(i)) {
202 indMax = i;
203 emax = cl->GetDigitAmplitude()->At(i);
204 }
205 }
206 }
207 if(emax/cl->E() > 0.5) { // more than 50% of cluster energy
208 u::FillH1(fLhists, 0, 2.); // number "good" entries
209 absIdMax = Int_t(absId[indMax]);
210 FillPi0Candidate(mgg, absIdMax, nm1);
211 }
212}
213
214void AliEMCALFolder::FillPi0Candidate(const Double_t mgg, Int_t absIdMax, Int_t nm)
215{
216 //
217 // Jun 08;
218 //
219 static Int_t nSupModMax, nModuleMax, nIphiMax, nIetaMax, iphiCellMax, ietaCellMax;
220 static AliEMCALCell* cell;
221 static TDataSet *set;
222 static AliEMCALSuperModule *sm;
223
224 fGeometry->GetCellIndex(absIdMax, nSupModMax, nModuleMax, nIphiMax, nIetaMax);
225 if(nm != nSupModMax) assert(0);
226
227 fGeometry->GetCellPhiEtaIndexInSModule(nSupModMax, nModuleMax, nIphiMax, nIetaMax, iphiCellMax, ietaCellMax);
228
229 cell = 0;
230 set = 0;
231 sm = 0;
232 if(GetCell(absIdMax)==0) {
233 cell = new AliEMCALCell(absIdMax, Form("sm%2.2i:phi%2.2i:eta%2.2i(%4.4i)",nSupModMax,iphiCellMax,ietaCellMax,absIdMax));
234 SetCell(cell, absIdMax);
235 // For browser
236 set = FindByName(Form("SM%2.2i",nSupModMax));
237 if(set==0) {
238 sm = new AliEMCALSuperModule(nSupModMax);
239 Add(sm);
240 sm->Init();
241 } else {
242 sm = (AliEMCALSuperModule*)set;
243 }
244 sm->AddCellToEtaRow(cell, ietaCellMax);
245 //
246 cell->SetCCfromCCTable(GetCCIn());
247 } else {
248 cell = GetCell(absIdMax);
249 set = FindByName(Form("SM%2.2i",nm));
250 if(set) sm = (AliEMCALSuperModule*)set;
251 }
252 if(sm == 0) assert(0);
253 if(nm != sm->GetSMNumber()) assert(0);
254
255 u::FillH1(sm->GetHists(), 0, mgg);
256 cell->FillEffMass(mgg);
257}
258
259void AliEMCALFolder::FitAllSMs()
260{ // Jun 14, 2007
261 // Only first SM now - should be changed in the future
262 AliEMCALSuperModule *sm0 = GetSuperModule(0);
263 sm0->FitForAllCells();
264 // Get input calibration table
265 AliEMCALCalibCoefs *ccIn = GetCCIn();
266 if(ccIn == 0) {
267 printf("<E> no input cc \n");
268 return;
269 }
270 // New calibration table
271 AliEMCALCalibCoefs *ccOut = new AliEMCALCalibCoefs(fgkCCoutName.Data(), ccIn->GetNRows());
272 calibCoef *rIn=0, rOut;
273 for(Int_t i=0; i<ccIn->GetNRows(); i++){
274 rIn = ccIn->GetTable(i);
275 rOut = (*rIn);
276 AliEMCALCell* cell = GetCell(rIn->absId);
277 if(cell && cell->GetSupMod() == 0) { // only first module now
278 rOut.cc = cell->GetCcOut();
279 }
280 ccOut->AddAt(&rOut);
281 }
282 ccOut->Purge();
283 Add(ccOut);
284}
285
286AliEMCALCalibCoefs* AliEMCALFolder::GetCCTable(const char* name)
287{
288 TDataSet *ds = FindByName(name);
289 if(ds) return (AliEMCALCalibCoefs*)ds;
290 else return 0;
291}
292
293Int_t AliEMCALFolder::GetSMNumber(AliESDCaloCluster* cl)
294{
295 static Int_t absId, nSupMod, nModule, nIphi, nIeta;
296 nSupMod = -1; // if negative something wrong
297 if(cl) {
298 absId = Int_t(cl->GetDigitIndex()->At(0));
299 fGeometry->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
300 }
301 return nSupMod;
302}
303
304// Recalibration staf - Jun 18,2007
305AliEMCALRecPoint* AliEMCALFolder::GetRecPoint(AliESDCaloCluster *cl, AliEMCALCalibCoefs *tOld,AliEMCALCalibCoefs *tNew,
306TList *l)
307{
308 //
309 // Static function;
310 // Get recalibrated rec.point from ESD cluster
311 //
312 // printf(" AliEMCALFolder::GetRecPoint() : RECALIBRATION \n");
313 static Double_t ADCCHANNELEC = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
314 static Float_t ECAW0 = 4.5; // hard coded now - see AliEMCALClusterizerv1::InitParameters()
315 Int_t ampDigi=0;
316 Double_t eDigiNew=0.0;
317
318 AliEMCALRecPoint *rp=0;
319 calibCoef *rOld=0, *rNew=0;
320 if(cl && tOld && tNew){
321 // cl->PrintClusterInfo(1);
322 const Int_t ndg = cl->GetNumberOfDigits();
323 //const Int_t nprim = cl->GetNumberOfPrimaries();
324
325 const Short_t* prim = cl->GetLabels()->GetArray();
326 const Short_t* dgAbsId = cl->GetDigitIndex()->GetArray();
327 // const UShort_t* dgAmp = cl->GetDigitAmplitude()(); // This is energy - bad definition
328
329 rp = new AliEMCALRecPoint(""); // opt=""
330 rp->SetClusterType(AliESDCaloCluster::kClusterv1);
331 AliEMCALDigit* dg=0;
332 TClonesArray digits("AliEMCALDigit", ndg);
333 Int_t absId = 0;
334 for(Int_t i=0; i<ndg; i++){
335 // Save just abs id and amplitude of digits which will be used for recalculation of
336 // cluster energy and position
337 absId = Int_t(dgAbsId[i]);
338 rOld = tOld->GetTable(absId);
339 rNew = tNew->GetTable(absId);
340 ampDigi = Int_t(cl->GetDigitAmplitude()->At(i) / rOld->cc);
341
342 new(digits[i]) AliEMCALDigit(Int_t(prim[i]),0,absId, ampDigi, 0.0, i, 0.0);
343 dg = (AliEMCALDigit*)digits[i];
344 eDigiNew = ampDigi * rNew->cc;
345 //eDigiNew = ampDigi * ADCCHANNELEC;
346 //eDigiNew = Double_t(cl->GetDigitAmplitude()->At(i)); // Copy from ESD for checking
347 rp->AddDigit(*dg, Float_t(eDigiNew));
348 u::FillH1(l, 14, eDigiNew);
349 u::FillH1(l, 15, Double_t(absId));
350 //printf("<I> digit %i amp %i rOld->cc %6.5f GeV rNew->cc %6.5f GeV\n", i, ampDigi, rOld->cc, rNew->cc);
351 }
352 //printf("<I> recalibration of digits was done ! \n");
353 // rp->EvalAll(ECAW0, &digits);
354 rp->EvalLocalPosition(ECAW0, &digits); // I need just position
355 digits.Delete();
356 }
357 //rp->Print("print");
358 return rp;
359}
360
361void AliEMCALFolder::Save(const char *fn, const char *opt)
362{
363 //
364 // Jun 5, 2007; See TFileIter and StFMC.cxx
365 //
366 TString FN(fn);
367 if(FN.Contains(".root")==0) FN += ".root";
368 TFileIter f(FN.Data(),opt,"EMCAL object");
369 UInt_t eventNum = 0; // just one object
370 UInt_t runNumber = 0; // 0 now, - may statistics on selectorr
371 f.NextEventPut(this, eventNum, runNumber);
372 printf(" Save %s to file %s\n", GetName(), FN.Data());
373}
374
375AliEMCALFolder* AliEMCALFolder::Read(const char *fn, const char *opt)
376{
377 //
378 // Jun 27, 2007
379 //
380 AliEMCALFolder* EMCAL = 0;
381 TH1::AddDirectory(0); // this is obligatory
382
383 TFile f(fn,opt);
384 TList *l = f.GetListOfKeys();
385 printf(" The total number of the objects: %d in file %s\n", l->GetSize(), fn);
386
387 TKey *key = (TKey*)l->At(0);
388 EMCAL = (AliEMCALFolder*)key->ReadObj();
389 f.Close();
390 if(EMCAL) EMCAL->InitAfterRead();
391
392 return EMCAL;
393}
394
395
396void AliEMCALFolder::InitAfterRead()
397{
398 lObj = (TList*)fObj;
399 fLhists = (TList*)lObj->FindObject("HistsOfEmcal");
400}
401
402void AliEMCALFolder::DrawQA(const int nsm)
403{
404 //
405 // Jun 25, 2007
406 //
407
408 AliEMCALSuperModule* sm = GetSuperModule(nsm);
409 if(sm==0) return;
410 TList *l = sm-> GetHists();
411 Int_t nx=2, ny=2, wh=530, ww=750;
412
413 TCanvas *c = new TCanvas(Form("QA_%i",GetIterationNumber()), Form("QA_%i",GetIterationNumber()),
414 10, 10, ww, wh);
415 c->Divide(nx,ny);
416
417 int ic=1;
418 c->cd(ic++);
419 TH1 *h1 = (TH1*)l->At(0);
420 sm->FitEffMassHist();
421 u::DrawHist(h1,2);
422 h1->SetAxisRange(0.03, 0.28);
423
424 c->cd(ic++);
425 sm->DrawCC(0);
426 TH1 *hccin = (TH1*)l->At(1);
427 hccin->SetAxisRange(14., 20.);
428
429 gStyle->SetOptStat(1111);
430 c->cd(ic++);
431 TH1 *hmass = (TH1*)l->At(3);
432 u::DrawHist(hmass, 2);
433 hmass->SetAxisRange(0.12, 0.16);
434
435 c->cd(ic++);
436 TH1 *hres = (TH1*)l->At(4);
437 u::DrawHist(hres, 2);
438 hres->SetAxisRange(0.05, 0.120);
439
440 if(ic<nx*ny) {
441 c->cd(ic++);
442 u::DrawHist((TH1*)l->At(5), 2);
443
444 c->cd(ic++);
445 u::DrawHist((TH1*)l->At(6), 2);
446 }
447 c->Update();
448}
449
450TList* AliEMCALFolder::BookHists()
451{
452 gROOT->cd();
453 TH1::AddDirectory(1);
454
455 new TH1F("00_HStat", "hist of common EMCAL statistics", 100, 0.5, 100.5);
456 new TH1F("01_EffMassAll", "effective mass of #gamma,#gamma(m_{#pi^{0}}=134.9766 MeV) - whole EMCAL", 250,0.0,0.5);
457
458 TList *l = AliEMCALHistoUtilities::MoveHistsToList("HistsOfEmcal", kFALSE);
459
460 TH1::AddDirectory(0);
461 return l;
462}
463
464void AliEMCALFolder::CreateCellNtuple()
465{
466 // Jun 28, 2007
467 if(fCellNtuple) { // Already exist
468 fCellNtuple->Print();
469 return;
470 }
471 // Create ntuple
472 Int_t bsize = int(1.e+5);
473 fCellNtuple = new TNtuple("cells","Cells Ntuple for quick analysis",
474 "fAbsId:fSupMod:fModule:fPhi:fEta:fPhiCell:fEtaCell:fCcIn:fCcOut", bsize);
475 AliEMCALCell *cell=0;
476 // Fill ntuple
477 AliEMCALSuperModule* sm = GetSuperModule(0);
478 if(sm) printf(" TNtuple was created ! sm0 %s \n", sm->GetName());
479 for(int eta=0; eta<48; eta++) { // eta row
480 TDataSet *setEta = sm->FindByName(Form("ETA%2.2i",eta));
481 if(setEta) {
482 printf(" ***** eta row %s ******\n", setEta->GetName());
483 for(int phi=0; phi<setEta->GetListSize(); phi++) { // cycle on cells (phi directions)
484 cell = (AliEMCALCell*)setEta->At(phi);
485 if(cell) {
486 cell->FillCellNtuple(fCellNtuple);
487 //printf(" fill cell %s : %s \n", cell->GetName(), cell->GetTitle());
488 }
489 }
490 }
491 }
492 fCellNtuple->Print();
493 lObj->Add(fCellNtuple);
494 // --
495 CreateAndFillAdditionalHists();
496}
497
498void AliEMCALFolder::CreateAndFillAdditionalHists()
499{
500 TH1::AddDirectory(0);
501 fLhists->Add(new TH1F("02_CCoutOnEdge", "cc out on edge of calorimeter (in MeV)", 70, 12., 19.));
502 fLhists->Add(new TH1F("03_CCoutInside", "cc out inside of calorimeter (in MeV)", 70, 12., 19.));
503 // Fill
504 Float_t* args;
505 for(Int_t i=0; i<(Int_t)fCellNtuple->GetEntries(); i++){
506 fCellNtuple->GetEvent(i);
507 args = fCellNtuple->GetArgs();
508 Int_t phi = (Int_t)args[5];
509 Int_t eta = (Int_t)args[6];
510 Double_t cc = (Double_t)args[8]*1000.;
511 if((phi==0||phi==23) || (eta==0||eta==47)) u::FillH1(fLhists, 2, cc);
512 else u::FillH1(fLhists, 3, cc);
513 }
514 // Drawing
515 Int_t wh=530, ww=750;
516 TCanvas *c = new TCanvas("c_edge","CEDGE", 10, 10, ww, wh);
517
518 gStyle->SetOptStat(1100);
519 gStyle->SetOptFit(111);
520 TH1 *h1 = (TH1*)fLhists->At(3);
521 TF1 *g = u::Gausi("ccInside", 14.7, 16.4, h1);
522 g->SetLineColor(kRed);
523 h1->Fit(g,"Q+","", 14.7, 16.4);
524 u::DrawHist(h1,2);
525 h1->SetTitle("CC distribution after #pi^{0} calibration");
526 h1->SetXTitle(" MeV ");
527 h1->SetYTitle(" N ");
528 TLatex *lat1 = u::lat(Form("rel.width = %4.2f%%",
529 100.*h1->GetRMS()/ h1->GetMean()), 16.5, 100., 12, 0.045);
530 TLatex *lat2 = u::lat(Form("rel.width = %4.2f%% (from fit)",
531 100.*g->GetParameter(2)/ g->GetParameter(1)), 16.5, 70., 12, 0.045);
532
533 if(0) {
534 TH1 *h2 = (TH1*)fLhists->At(2);
535 u::DrawHist(h2,2,1,"same",2);
536 }
537
538 TH1F *hccFirst = AliEMCALCalibCoefs::GetHistOfCalibTableFromDb("ccTmp");
539 u::DrawHist(hccFirst,2,1,"same",3);
540
541
542 TLegend *leg = new TLegend(0.1,0.6, 0.45,0.85);
543 leg->AddEntry(hccFirst, "Initial cc ", "L");
544 leg->AddEntry(h1, "Final cc", "L");
545 leg->Draw();
546
547 c->Update();
548}
549
550void AliEMCALFolder::TestSMStruct()
551{
552 // testing May 22, 2007
553 for(int m=0; m<12; m++) {
554 AliEMCALSuperModule *sm = new AliEMCALSuperModule(m);
555 Add(sm);
556 }
557}
558