The data members fEMCALClusterCluster and fPHOSCluster are removed from AliESDCaloClu...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecPointsQaESDSelector.cxx
CommitLineData
16d3c94d 1/**************************************************************************
0fc11500 2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
16d3c94d 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
8ada0ffe 16/* $Log$
17/* Revision 1.2 2007/09/11 19:38:15 pavlinov
18/* added pi0 calibration, linearity, shower profile
19/*co: warning: `/* $Log' is obsolescent; use ` * $Log'.
0fc11500 20 * Revision 1.1 2007/08/08 15:58:01 hristov
21 * New calibration classes. (Aleksei)
22 * */
16d3c94d 23//_________________________________________________________________________
24// This is selector for making a set of histogramms from ESD for rec.points
25//*-- Authors: Aleksei Pavlinov (WSU)
0fc11500 26//*-- Feb 17, 2007 - Sep 11, 2007
16d3c94d 27//*-- Pi0 calibration
0fc11500 28//*-- Recalibration, linearity, shower profile
16d3c94d 29//_________________________________________________________________________
30
31#include "AliEMCALRecPointsQaESDSelector.h"
32#include "AliEMCALFolder.h"
33#include "AliEMCALSuperModule.h"
34#include "AliEMCALCell.h"
0fc11500 35#include "AliEMCALCellInfo.h"
16d3c94d 36#include "AliEMCALPi0SelectionParam.h"
37#include "AliEMCALHistoUtilities.h"
38#include "AliEMCALCalibCoefs.h"
39#include "AliEMCALGeometry.h"
40#include "AliLog.h"
41#include "AliESD.h"
42#include "AliEMCALRecPoint.h"
0fc11500 43#include "AliRunLoader.h"
44#include "AliStack.h"
45#include "jetfinder/AliEMCALJetMicroDst.h" // Sgpdge
46#include "AliEMCALDigit.h"
47
48#include <iostream>
49#include <iomanip>
50#include <fstream>
16d3c94d 51#include <assert.h>
52#include <map>
53#include <string>
54
55#include <TROOT.h>
56#include <TStyle.h>
57#include <TSystem.h>
58#include <TCanvas.h>
59#include <TVector3.h>
60#include <TLorentzVector.h>
0fc11500 61#include <TParticle.h>
16d3c94d 62#include <TChain.h>
63#include <TFile.h>
64#include <TH1F.h>
0fc11500 65#include <TH1D.h>
66#include <TH2F.h>
16d3c94d 67#include <TF1.h>
68#include <TMath.h>
16d3c94d 69#include <TBrowser.h>
0fc11500 70#include <TList.h>
71#include <TCanvas.h>
72#include <TLine.h>
73#include <TObjString.h>
74#include <TArrayI.h>
75#include <TGraphErrors.h>
76#include <TLegend.h>
77#include <TLegendEntry.h>
78#include <TNtuple.h>
16d3c94d 79
80using namespace std;
81
82typedef AliEMCALHistoUtilities u;
83
0fc11500 84Double_t AliEMCALRecPointsQaESDSelector::fDistEff = 3.;
85Double_t AliEMCALRecPointsQaESDSelector::fW0 = 4.5;
86Double_t AliEMCALRecPointsQaESDSelector::fSlopePhiShift = 0.01; // ?? just guess
87
88AliEMCALFolder* AliEMCALRecPointsQaESDSelector::fEMCAL = 0;
89AliEMCALFolder* AliEMCALRecPointsQaESDSelector::fEMCALOld = 0;
90
91AliEMCALGeometry *geo=0;
92
93char *anaOpt[]={
94 "CORR1", // GetCorrectedEnergyForGamma_1(Double_t eRec);
95 "RECALIB",
96 "IDEAL",
97 "PI0",
98 "GAMMA",
99 "KINE", // reading kine file
100 "PROF", // Shower profile: phi direction now
101 "FIT" // define parameters : deff, w0 and phislope
102};
103enum keyOpt{
104 kCORR1,
105 kRECALIB,
106 kIDEAL,
107 kPI0,
108 kGAMMA,
109 kKINE,
110 kPROF,
111 kFIT,
112 //--
113 kEND
114};
115
116int nAnaOpt = sizeof(anaOpt) / sizeof(char*);
117// Enumeration variables
118// enum {kUndefined=-1, kLed=0, kBeam=1};
119// enum {kPatchResponse=0, kClusterResponse=1};
120// enum {kSkip=0, kSaveToMemory=1};
121
16d3c94d 122ClassImp(AliEMCALRecPointsQaESDSelector)
123
124Int_t nmaxCell = 4*12*24*11;
125
126AliEMCALRecPointsQaESDSelector::AliEMCALRecPointsQaESDSelector() :
127 AliSelector(),
0fc11500 128 fPmom(0),
16d3c94d 129 fChain(0),
130 fLofHistsPC(0),
131 fLofHistsRP(0),
0fc11500 132 fLKineVsRP(0),
133 fLShowerProfile(0),
134 fCellsInfo(0),
16d3c94d 135 fEmcalPool(0),
0fc11500 136 fRunOpts(0),
137 fArrOpts(0),
138 fKeyOpts(0)
16d3c94d 139{
140 //
141 // Constructor. Initialization of pointers
142 //
143}
144
145AliEMCALRecPointsQaESDSelector::~AliEMCALRecPointsQaESDSelector()
146{
147 //
148 // Destructor
149 //
150
151 // histograms are in the output list and deleted when the output
152 // list is deleted by the TSelector dtor
153}
154
155void AliEMCALRecPointsQaESDSelector::Begin(TTree* tree)
156{
157 // Begin function
158
159 AliSelector::Begin(tree);
160}
161
16d3c94d 162void AliEMCALRecPointsQaESDSelector::SlaveBegin(TTree* tree)
163{
164 // The SlaveBegin() function is called after the Begin() function.
165 // When running with PROOF SlaveBegin() is called on each slave server.
166 // The tree argument is deprecated (on PROOF 0 is passed).
167
168 printf("**** <I> Slave Begin \n **** ");
169
170 AliSelector::SlaveBegin(tree);
171}
172
0fc11500 173void AliEMCALRecPointsQaESDSelector::Init(TTree* tree)
174{
175 // Init function
176
177 AliSelector::Init(tree);
178}
179Bool_t AliEMCALRecPointsQaESDSelector::Notify()
180{
181 // The Notify() function is called when a new file is opened. This
182 // can be either for a new TTree in a TChain or when when a new TTree
183 // is started when using PROOF. Typically here the branch pointers
184 // will be retrieved. It is normaly not necessary to make changes
185 // to the generated code, but the routine can be extended by the
186 // user if needed.
187
188 return AliSelector::Notify();
189}
190
16d3c94d 191void AliEMCALRecPointsQaESDSelector::InitStructure(Int_t it)
192{
0fc11500 193 geo = AliEMCALGeometry::GetInstance("SHISH_TRD1_CURRENT_2X2"); // initialize geometry just once
194 fCellsInfo = AliEMCALCellInfo::GetTableForGeometry(geo);
195
196 if(fRunOpts.Length()>0) CheckRunOpts();
16d3c94d 197
0fc11500 198 Int_t key = 0;
199 if(GetKeyOptsValue(kGAMMA)) key = 1;
200 fLofHistsRP = DefineHistsOfRP("RP", fPmom, key);
201 fLofHistsPC = DefineHistsOfRP("PseudoCl", fPmom);
16d3c94d 202
0fc11500 203 fEmcalPool = new TFolder("PoolOfEMCAL","");
204 if(it <= 1) {
16d3c94d 205 fEMCAL = new AliEMCALFolder(it); // folder for first itteration
206 fEmcalPool->Add(fEMCAL);
207 }
0fc11500 208 //if(it<=0) SetName("GammaSel"); // For convinience
209 //else SetName("Pi0Sel");
210 if(GetKeyOptsValue(kKINE)) fLKineVsRP = DefineHistsOfKineVsRP("KineVsRP",fPmom, key);
211 if(GetKeyOptsValue(kPROF)) fLShowerProfile = DefineHistsForShowerProfile("ProfY", fPmom);
212}
213
214void AliEMCALRecPointsQaESDSelector::CheckRunOpts()
215{
216 fRunOpts.ToUpper();
217 int nopt = u::ParseString(fRunOpts, fArrOpts);
218 printf("<I> AliEMCALRecPointsQaESDSelector::CheckRunOpts() analyze %i(%i) options : nAnaOpt %i\n",
219 nopt, fArrOpts.GetEntries(), nAnaOpt);
220 if(nopt <=0) return;
221
222 fKeyOpts = new TArrayI(nAnaOpt);
223 for(Int_t i=0; i<nAnaOpt; i++) (*fKeyOpts)[i] = 0;
224
225 for(Int_t i=0; i<fArrOpts.GetEntries(); i++ ) {
226 TObjString *o = (TObjString*)fArrOpts.At(i);
227
228 TString runOpt = o->String();
229 Int_t indj=-1;
230
231 for(Int_t j=0; j<nAnaOpt; j++) {
232 TString opt = anaOpt[j];
233 if(runOpt.Contains(opt,TString::kIgnoreCase)) {
234 indj = j;
235 break;
236 }
237 }
238 if(indj<0) {
239 printf("<E> option |%s| unavailable **\n",
240 runOpt.Data());
241 assert(0);
242 } else {
243 (*fKeyOpts)[indj] = 1;
244 printf("<I> option |%s| is valid : number %i : |%s| \n",
245 runOpt.Data(), indj, anaOpt[indj]);
246 }
247 }
248}
249
250Int_t AliEMCALRecPointsQaESDSelector::GetKeyOptsValue(Int_t key)
251{
252 static Int_t val=0;
253 val = 0;
254 if(fKeyOpts && key>=0 && key<fKeyOpts->GetSize()) {
255 val = fKeyOpts->At(key);
256 }
257 // printf(" key %i : val %i : opt %s \n", key, val, anaOpt[key]);
258 return val;
16d3c94d 259}
260
261Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
262{
263 // The Process() function is called for each entry in the tree (or possibly
264 // keyed object in the case of PROOF) to be processed. The entry argument
265 // specifies which entry in the currently loaded tree is to be processed.
266 // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
267 // to read either all or the required parts of the data. When processing
268 // keyed objects with PROOF, the object is already loaded and is available
269 // via the fObject pointer.
270 //
271 // This function should contain the "body" of the analysis. It can contain
272 // simple or elaborate selection criteria, run algorithms on the data
273 // of the event and typically fill histograms.
274
275 // WARNING when a selector is used with a TChain, you must use
276 // the pointer to the current TTree to call GetEntry(entry).
277 // The entry is always the local entry number in the current tree.
278 // Assuming that fTree is the pointer to the TChain being processed,
279 // use fTree->GetTree()->GetEntry(entry).
280
0fc11500 281 if (AliSelector::Process(entry) == kFALSE) {
282 AliDebug(AliLog::kError, Form("Entry %i is not available ",entry));
16d3c94d 283 return kFALSE;
0fc11500 284 }
285
286 // esd->ReadFromTree(tree);
287 // AliESDEvent * esd = new AliESDEvent();
16d3c94d 288
289 // Check prerequisites
290 if (!fESD)
291 {
292 AliDebug(AliLog::kError, "ESD branch not available");
293 return kFALSE;
294 }
0fc11500 295
16d3c94d 296 static pi0SelectionParam* rPar = GetEmcalFolder()->GetPi0SelectionParRow(0);
297
298 static Int_t nEmcalClusters, indOfFirstEmcalRP, nEmcalRP,nEmcalPseudoClusters;
299 nEmcalClusters = fESD->GetNumberOfEMCALClusters();
300 indOfFirstEmcalRP = fESD->GetFirstEMCALCluster();
301 u::FillH1(fLofHistsRP, 1, double(indOfFirstEmcalRP));
302
0fc11500 303 static AliRunLoader* RL = 0;
304 static Int_t nev = 0; // Temporary - 0nly for reading one file now !!
305
16d3c94d 306 static AliESDCaloCluster *cl = 0;
307 nEmcalRP = nEmcalPseudoClusters = 0;
308 TList *l=0;
309 double eDigi=0;
310
311 TClonesArray lvM1("TLorentzVector", 100); // for convenience; 100 is enough now
312 TArrayI indLv(100); // index of RP
313
314 static TLorentzVector v, vcl;
315 int nrp = 0; // # of RP for gg analysis
0fc11500 316 static Double_t erec=0., ecorr=0.0;
16d3c94d 317 for(int i=indOfFirstEmcalRP; i<indOfFirstEmcalRP+nEmcalClusters; i++) {
318 cl = fESD->GetCaloCluster(i);
8ada0ffe 319 if(cl->GetClusterType() == AliESDCaloCluster::kEMCALPseudoCluster) {
16d3c94d 320 nEmcalPseudoClusters++;
321 l = fLofHistsPC;
8ada0ffe 322 } else if(cl->GetClusterType() == AliESDCaloCluster::kEMCALClusterv1){
16d3c94d 323 nEmcalRP++;
0fc11500 324 if(fEMCAL->GetIterationNumber()>1||GetKeyOptsValue(kIDEAL)||GetKeyOptsValue(kRECALIB)||GetKeyOptsValue(kFIT)) {
325 AliEMCALRecPoint *rp=0;
326 if(GetKeyOptsValue(kFIT) == kFALSE) fDistEff = -1.; // No fitting ; Sep 4, 2007
327 if(GetKeyOptsValue(kIDEAL)) {
328 rp = AliEMCALFolder::GetRecPoint(cl, fEMCAL->GetCCFirst(), 0, fLofHistsRP, fDistEff, fW0, fSlopePhiShift);
329 } else {
330 rp = AliEMCALFolder::GetRecPoint(cl, fEMCAL->GetCCFirst(), fEMCAL->GetCCIn(), fLofHistsRP, fDistEff, fW0, fSlopePhiShift);
331 }
332 if(GetKeyOptsValue(kPROF)) {
333 FillHistsForShowerProfile( GetListShowerProfile(), rp, GetCellsInfo());
334 }
335 //if(rp->GetPointEnergy()>=rPar->eOfRpMin && u::GetLorentzVectorFromRecPoint(v, rp)) {
16d3c94d 336 if(u::GetLorentzVectorFromRecPoint(v, rp)) { // comparing with RP
0fc11500 337 if(GetKeyOptsValue(kCORR1)) {
338 erec = v.Rho();
339 ecorr = u::GetCorrectedEnergyForGamma_1(erec);
340 v.SetRho(ecorr);
341 v.SetE(ecorr); // This is gamma
342 //printf("<1> erec %f | ecorr %f \n", erec, ecorr);
343 }
16d3c94d 344 new(lvM1[nrp]) TLorentzVector(v);
345 indLv[nrp] = i;
346 nrp++;
347 // Conroling of recalibration
348 u::GetLorentzVectorFromESDCluster(vcl, cl);
349 u::FillH1(fLofHistsRP, 11, vcl.P()-v.P());
350 u::FillH1(fLofHistsRP, 12, TMath::RadToDeg()*(vcl.Theta()-v.Theta()));
351 u::FillH1(fLofHistsRP, 13, TMath::RadToDeg()*(vcl.Phi()-v.Phi()));
0fc11500 352
353 u::FillH1(fLofHistsRP, 16, v.P());
354 u::FillH2(fLofHistsRP, 17, vcl.P(), vcl.P()-v.P());
16d3c94d 355 l = 0; // no filling
0fc11500 356 if(GetKeyOptsValue(kIDEAL) || GetKeyOptsValue(kRECALIB)) l = fLofHistsRP;
16d3c94d 357 }
0fc11500 358 if(rp) delete rp;
16d3c94d 359 } else { // first iteration
360 // if(cl->E()>=rPar->eOfRpMin && u::GetLorentzVectorFromESDCluster(v, cl)) {
361 if(u::GetLorentzVectorFromESDCluster(v, cl)) { // comparing with RP
362 // cut 0.4 GeV may be high !
0fc11500 363 if(GetKeyOptsValue(kCORR1)) {
364 erec = v.Rho();
365 ecorr = u::GetCorrectedEnergyForGamma_1(erec);
366 v.SetRho(ecorr);
367 v.SetE(ecorr); // This is gamma now
368 // printf("<2> erec %f | ecorr %f \n", erec, ecorr);
369 // printf(" v.Rho() %f \n", v.Rho());
370 }
16d3c94d 371 new(lvM1[nrp]) TLorentzVector(v);
372 indLv[nrp] = i;
373 nrp++;
374 l = fLofHistsRP;
375 }
376 }
377
378 } else {
379 printf(" wrong cluster type : %i\n", cl->GetClusterType());
380 assert(0);
381 }
382 u::FillH1(l, 2, double(cl->GetClusterType()));
383
384 u::FillH1(l, 3, double(cl->GetNumberOfDigits()));
385 u::FillH1(l, 4, double(cl->E()));
386 // Cycle on digits (towers)
387 Short_t *digiAmpl = cl->GetDigitAmplitude()->GetArray();
388 Short_t *digiTime = cl->GetDigitTime()->GetArray();
389 Short_t *digiAbsId = cl->GetDigitIndex()->GetArray();
390 for(int id=0; id<cl->GetNumberOfDigits(); id++) {
391 eDigi = double(digiAmpl[id]) / 500.; // See AliEMCALClusterizerv1
392 // if(eDigi <= 0.0) { // sometimes it is happen
8ada0ffe 393 //if(eDigi > 10.0 && cl->GetClusterType() == AliESDCaloCluster::kEMCALClusterv1) {
16d3c94d 394 // printf(" %i digiAmpl %i : %f \n", id, int(digiAmpl[id]), eDigi);
395 //}
396 u::FillH1(l, 5, eDigi);
397 u::FillH1(l, 6, double(digiTime[id]));
398 u::FillH1(l, 7, double(digiAbsId[id]));
399 if(int(digiAbsId[id]) >= nmaxCell) {
400 printf(" id %i : digiAbsId[id] %i (%i) : %s \n",
401 id, int(digiAbsId[id]), nmaxCell, l->GetName());
402 }
403 }
404 }
405 u::FillH1(fLofHistsRP, 0, double(nEmcalRP));
406 u::FillH1(fLofHistsPC, 0, double(nEmcalPseudoClusters));
407
408 static TLorentzVector *lv1=0, *lv2=0, lgg;
0fc11500 409 for(int i1=0; i1<nrp; i1++){
410 lv1 = (TLorentzVector*)lvM1.At(i1);
411 u::FillH1(fLofHistsRP, 18, lv1->P());
412 }
16d3c94d 413 static Double_t mgg, pgg;
414 mgg = pgg = 0.;
415 nrp = lvM1.GetEntriesFast();
416 if(nrp >= 2) {
417 // eff.mass analysis
418 for(int i1=0; i1<nrp-1; i1++){
419 lv1 = (TLorentzVector*)lvM1.At(i1);
420 for(int i2=i1+1; i2<nrp; i2++){
421 lv2 = (TLorentzVector*)lvM1.At(i2);
422 lgg = (*lv1) + (*lv2);
423 mgg = lgg.M(); // eff.mass
424 pgg = lgg.P(); // momentum
425 u::FillH1(fLofHistsRP, 8, mgg);
426
427 if((mgg>=rPar->massGGMin && mgg<=rPar->massGGMax)) {// pi0 candidates
428 if((pgg>=rPar->momPi0Min && pgg>=rPar->momPi0Min)) {
0fc11500 429 if(fEMCAL && fEMCAL->GetIterationNumber()>=1) {
430 fEMCAL->FillPi0Candidate(mgg,fESD->GetCaloCluster(indLv[i1]),fESD->GetCaloCluster(indLv[i2]));
431 u::FillH1(fLofHistsRP, 9, pgg);
432 u::FillH1(fLofHistsRP,10, lv1->P());
433 u::FillH1(fLofHistsRP,10, lv2->P());
434 }
16d3c94d 435 }
436 }
437 }
438 }
439 }
440
0fc11500 441 // static Int_t fileNumber = 0;
442 static TString curFileName;
443 // if(GetKeyOptsValue(kKINE) && nev<fChain->GetEntries()) {
444 if(GetKeyOptsValue(kKINE)) {
445 // Get galice.root file in current directory
446 if(nev%1000==0) {
447 printf(" current file |%s|\n", fChain->GetCurrentFile()->GetName());
448 curFileName = fChain->GetCurrentFile()->GetName();
449 curFileName.ReplaceAll("AliESDs.","galice.");
450 }
451 RL = u::InitKinematics(nev, curFileName.Data());
452 // Compare kineamtics vs EMCal clusters
453 FillHistsOfKineVsRP(fLKineVsRP, RL, lvM1);
454 }
455
16d3c94d 456 lvM1.Delete();
457
458 if(nEmcalClusters != (nEmcalRP+nEmcalPseudoClusters))
459 Info("Process","nEmcalClusters %i : RP %i + PC %i ",nEmcalClusters, nEmcalRP, nEmcalPseudoClusters);
460
0fc11500 461 nev++;
462
16d3c94d 463 return kTRUE;
464}
465
466void AliEMCALRecPointsQaESDSelector::SlaveTerminate()
467{
468 // The SlaveTerminate() function is called after all entries or objects
469 // have been processed. When running with PROOF SlaveTerminate() is called
470 // on each slave server.
471
472 AliSelector::SlaveTerminate();
473
474 // Add the histograms to the output on each slave server
475 if (!fOutput)
476 {
477 AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
478 return;
479 }
480
481 // fOutput->Add(fMultiplicity);
482}
483
484void AliEMCALRecPointsQaESDSelector::Terminate()
485{
486 // The Terminate() function is the last function to be called during
487 // a query. It always runs on the client, it can be used to present
488 // the results graphically or save the results to file.
489
490 AliSelector::Terminate();
16d3c94d 491
16d3c94d 492 PrintInfo();
493}
0fc11500 494
16d3c94d 495//
0fc11500 496TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfRP(const char *name,Double_t p,Int_t keyOpt)
16d3c94d 497{
0fc11500 498 printf("<I> DefineHistsOfRP :%s : p %f : keyOpt %i \n", name, p, keyOpt);
499 Double_t ADCchannelEC = 0.0153; // ~15mev per adc count
500 Double_t xma = p*1.4, xmi=0.0, step=0.0, xmic=xmi, xmac = xma;
501 if(xma<0) xma = 20.;
502 Int_t nmax=1000, scale=4, nmaxc = nmax;
16d3c94d 503
504 gROOT->cd();
505 TH1::AddDirectory(1);
506 new TH1F("00_EmcalMultiplicity", "multiplicity of EMCAL ", 201, -0.5, 200.5); // real and pseudo
507 new TH1F("01_IndexOfFirstEmcal", "index of first emcal rec.points ", 201, -0.5, 200.5);
508
509 new TH1F("02_NumberOf", "number of ", 6, -0.5, 5.5);
510 new TH1F("03_NumberOfDigitsIn", "number of digits(towers) in rec.points ", 101, -0.5, 100.5);
16d3c94d 511
0fc11500 512 if(keyOpt==1 && p>0.1) {
513 if (p>=100.) scale=12;
514 else if(p>=50.) scale=10;
515 else if(p>=30.) scale=8;
516 xma = p + scale*(0.15*TMath::Sqrt(p));
517 xmi = p - scale*(0.15*TMath::Sqrt(p));
518 step = (xma-xmi)/nmax;
519 }
520
521 if(step < 0.0153) {
522 nmax = int((xma-xmi) / ADCchannelEC)+1;
523 xma = xmi + ADCchannelEC*nmax;
524 }
525 new TH1F("04_EnergyOf", "energy of ", nmax, xmi, xma);
526 nmaxc = nmax; xmic=xmi; xmac = xma;
527
528 nmax = 10000;
529 xmi = ADCchannelEC/2.; xma = xmi + ADCchannelEC*nmax;
16d3c94d 530 // All energy(momentum) unit is GeV if don't notice
0fc11500 531 new TH1F("05_DigitEnergyIn", "digit energy in ", nmaxc, xmic, xmac);
16d3c94d 532 new TH1F("06_DigitTimeIn", "digit time in 10ps(0.01ns) ", 1000, 0.0, 3.e+3); // ns/100 = 10 ps
533 new TH1F("07_DigitAbsIdIn", "digit abs id in ", nmaxCell, -0.5, double(nmaxCell)-0.5);
534 new TH1F("08_EffMass", "effective mass of #gamma,#gamma(m_{#pi^{0}}=134.9766 MeV)", 100, 0.0, 0.5);
535 new TH1F("09_MomOfPi0Candidate", "momentum of #pi^{0} candidates (0.085 <mgg<0.185)", 600, 0.0, 30.0);
536 new TH1F("10_MomOfRpPi0Candidate", "momentum of RP for #pi^{0} candidates (0.085 <mgg<0.185)", 600, 0.0, 30.0);
537 // Recalibration staf
0fc11500 538 Double_t pmax=11., dpmax=0.20;
539 if(p>0.1) {
540 pmax=p*1.2;
541 dpmax *= p;
542 }
543 new TH1F("11_MomClESD-RpRecalib", "difference of momentum cl(ESD) - rp(Recalib)", 100, -.01, dpmax);
16d3c94d 544 new TH1F("12_ThetaClESD-RpRecalib", "difference of #theta cl(ESD) - rp(Recalib) in degree", 100, -0.05, +0.05);
545 new TH1F("13_PhiClESD-RpRecalib", "difference of #phi cl(ESD) - rp(Recalib) in degree ", 100, -0.05, +0.05);
546 // Digi
547 new TH1F("14_EDigiRecalib", "energy of digits after recalibration", 2000, 0.0, 20.);
0fc11500 548 // AliEMCALGeometry* g = AliEMCALGeometry::GetInstance();
549 new TH1F("15_AbsIdRecalib", "abs Id of digits after recalibration", geo->GetNCells(),-0.5,Double_t(geo->GetNCells())-0.5);
550 new TH1F("16_EnergyOfRecalibRp_", "energy of recalibrated rec.points", nmaxc, xmic, xmac); // Jul 12, 2007
551 new TH2F("17_ShiftRecalib_", "E(clESD) - E(recalib)", 110,0.0, pmax, 50,0.0,dpmax); // Jul 13, 2007
16d3c94d 552
0fc11500 553 // Corrected staff
554 new TH1F("18_EnergyOfCorrectedLV", "energy of corrected LV", nmaxc, xmic, xmac);
555
556 TString st = Form("ListOfHists%s_P=%5.1f",name,p);
557 st.ReplaceAll(" ","");
558 TList *l = u::MoveHistsToList(st.Data(), kFALSE);
559 st = Form("%s_P=%5.1f",name,p);
560 st.ReplaceAll(" ","");
561 u::AddToNameAndTitleToList(l, st.Data(), st.Data());
562
563 return l;
564}
565
566TList* AliEMCALRecPointsQaESDSelector::DefineHistsOfKineVsRP(const char *name, Double_t p, Int_t keyOpt)
567{
568 printf("<I> DefineHistsOfKineVsRP :%s : p %f : keyOpt %i \n", name, p, keyOpt);
569
570 gROOT->cd();
571 TH1::AddDirectory(1);
572 new TH1F("00_hVx",Form("Vx of primary vertex"), 100, -5., +5.); // 00
573 new TH1F("01_hVy",Form("Vy of primary vertex"), 100, -5., +5.); // 01
574 new TH1F("02_hVz",Form("Vz of primary vertex"), 100, -50., +50.); // 02
575
576 // Double_t ADCchannelEC = 0.0153; // ~15mev per adc count
577 Double_t xma = p*1.4, xmi=0.0, sig=0.15*TMath::Sqrt(p);
578 // Double_t step=0.0, xmic=xmi, xmac = xma;
579 if(xma<0) xma = 20.;
580 //Int_t nmax=1000;
581 // scale=4, nmaxc = nmax;
582
583 new TH1F("03_hGidPrimar", "Geant Id of primary particle ", 101, 0.5, 100.5); // 03
584 new TH1F("04_hPmomPrim","momentum of primary particle", 100, xmi, xma); // 04
585 new TH1F("05_hEtaPrim","#eta of primary particle ", 200, -1., 1.); // 05
586 new TH1F("06_hPhiPrim","#phi of primary particle ", 63*2, 0.0,TMath::Pi()*2.); // 06
587 new TH1F("07_hThetaPrim","#theta of primary particle", 63, 0.0,TMath::Pi()); // 07
588 new TH1F("08_hPhiPrimInDegree","#phi of primary particle in degree", 120, 75.0, 195.0); // 08
589 new TH1F("09_hThetaPrimInDegree","#theta of primary particle in degree", 90, 45., 135.); // 09
590
591 // Particle vs cluster
592 new TH1F("10_hE(P)-E(RP)"," E(p) - E(RP) ", 100, -10.*sig, 10.*sig); // 10
593
594 new TH1F("11_hPvsRpAngleInDegree","angle between P and RP (in degree)", 110, 0.0, 1.0); // 11
595 double dphi=0.5; // for fitting
596 new TH1F("12_hPvsRpDPhiInDegree","dif(#phi) between P and RP (in degree)", 100, -dphi, dphi); // 12
597 new TH1F("13_hPvsRpDThetaInDegree","dif(#theta) between P and RP (in degree)", 100, -0.5, 0.5); // 13
598
599 new TH1F("14_hPvsRpDPhiInDegree","dif(#phi) between P and RP (in degree) right part of SM", 100, -dphi, +dphi); // 14
600 new TH1F("15_hPvsRpDPhiInDegree","dif(#phi) between P and RP (in degree) left part of SM", 100, -dphi, dphi); // 15
601
602 new TH1F("16_hPvsRpDThetaInDegree","dif(#theta) between P and RP (even index)", 100, -0.5, 0.5); // 16
603 new TH1F("17_hPvsRpDThetaInDegree","dif(#theta) between P and RP (odd index)", 100, -0.5, 0.5); // 17
604
605 TString st = Form("ListOfHists%s_P=%5.1f",name,p);
606 st.ReplaceAll(" ","");
607 TList *l = u::MoveHistsToList(st.Data(), kFALSE);
608 st = Form("%s_P=%5.1f",name,p);
609 st.ReplaceAll(" ","");
610 u::AddToNameAndTitleToList(l, st.Data(), st.Data());
611 return l;
612}
613
614TList *AliEMCALRecPointsQaESDSelector::DefineHistsForShowerProfile(const char *name, Double_t p)
615{
616 // Aug 1, 2007 - shower profile business as in testbeam analysis
617 // Phi direction here is Y direction in testbeam
618 printf("<I> DefineHistsForShowerProfile: %s : p %f \n",
619 name, p);
620
621 gROOT->cd();
622 TH1::AddDirectory(1);
623
624 // Cell unit
625 new TH1F("00_hPHiMean"," mean in #phi direction ", 24*10, 0.0, 24.);
626 new TH1F("01_hPHiLog"," mean with log in #phi direction ", 24*10, 0.0, 24.);
627 new TH2F("02_XmeanVsXlog", "xmean vs xlog", 240,0.0, 24.,240,0.0, 24.);
628 new TH2F("03_XmeanVsXlog", "xmean vs xlog (system of cell with max energy, first half on #eta)",
629 50,-0.5,+0.5, 100,-1.,+1.);
630 new TH2F("04_XmeanVsXlog", "xmean vs xlog (system of cell with max energy, second half on #eta)",
631 50,-0.5,+0.5, 100,-1.,+1.);
632
633 new TH2F("05_PhiShowerProfile", " #phi shower profile - cumulative distribution",
634 6*25,-3.,3., 201, -0.0025, 1.0025); // 00; x direction in cell unit
635
636 TString st = Form("L%s_P=%5.1f",name,p);
637 st.ReplaceAll(" ","");
638 TList *l = u::MoveHistsToList(st.Data(), kFALSE);
639 st = Form("%s_P=%5.1f",name,p);
640 st.ReplaceAll(" ","");
641 u::AddToNameAndTitleToList(l, st.Data(), st.Data());
16d3c94d 642
643 return l;
644}
645
0fc11500 646void AliEMCALRecPointsQaESDSelector::FillHistsOfKineVsRP(TList *l, AliRunLoader* RL, TClonesArray &lvM)
647{
648 //
649 // lvM - array of TLorentzVector's which was cretaef from AliESDCaloCluster's
650 //
651
652 if(l==0 || RL==0) return;
653
654 // TNtuple for qucik analysis
655 static TNtuple *nt=0;
656 if(nt==0) {
657 gROOT->cd();
658 Int_t bsize = int(1.e+7);
659 nt = new TNtuple("angle","angle ntuple for quick analysis",
660 "dPhi:dTheta:phiP:thetaP", bsize);
661 gROOT->GetListOfBrowsables()->Add(nt);
662 }
663 static AliStack* st=0;
664 static TParticle *p=0;
665 static Int_t gid=0, ic=0, pdg=0, i=0;
666 gid = ic = pdg = 0;
667
668 st = RL->Stack();
669 if(st == 0) return;
670 // first primary particle
671 p = st->Particle(0);
672 if(p == 0) return;
673
674 u::FillH1(l, ic++, p->Vx());
675 u::FillH1(l, ic++, p->Vy());
676 u::FillH1(l, ic++, p->Vz());
677
678 pdg = p->GetPdgCode();
679 //gid = gMC->IdFromPDG(pdg); // gMc should be defined
680 AliEMCALJetMicroDst::Sgpdge(pdg,gid); // temporary
681
682 u::FillH1(l, ic++, Double_t(gid));
683 u::FillH1(l, ic++, p->P());
684 u::FillH1(l, ic++, p->Eta());
685 u::FillH1(l, ic++, TVector2::Phi_0_2pi(p->Phi()) );
686 u::FillH1(l, ic++, p->Theta());
687 // In degree
688 u::FillH1(l, ic++, TVector2::Phi_0_2pi(p->Phi())*TMath::RadToDeg());
689 u::FillH1(l, ic++, p->Theta()*TMath::RadToDeg()); // 09
690
691 if(lvM.GetEntriesFast() == 0) return;
692 //
693 // Initial kinematics vs Calo clusters
694 //
695 static TLorentzVector *lv = 0, lvp;
696 lvp.SetPxPyPzE(p->Px(),p->Py(),p->Pz(),p->Energy()); // Particle
697
698 Double_t angle = 0.0, angleMin = 180.0, eDiff=0.0, dPhi=0., dTheta=0., phiP=0.0, thetaP=0.0;
699 Int_t indMin = 0;
700 phiP = lvp.Phi() * TMath::RadToDeg();
701 if (phiP<81. || phiP>99.) return; // cut phi boundaries
702 thetaP = lvp.Theta() * TMath::RadToDeg();
703 if (thetaP<56. || thetaP>89.) return; // cut theta boundaries
704
705 for(i=0; i<lvM.GetEntriesFast(); i++) {
706 lv = (TLorentzVector*)lvM.At(i);
707 angle = lv->Angle(lvp.Vect())*TMath::RadToDeg();
708 if(angleMin > angle) {
709 angleMin = angle;
710 indMin = i;
711 }
712 }
713 lv = (TLorentzVector*)lvM.At(indMin);
714 eDiff = lvp.E() - lv->E();
715 u::FillH1(l, ic++, eDiff);
716
717 dPhi = TVector2::Phi_mpi_pi(lvp.Phi()-lv->Phi())*TMath::RadToDeg();
718 dTheta = TVector2::Phi_mpi_pi(lvp.Theta()-lv->Theta())*TMath::RadToDeg();
719 u::FillH1(l, ic++, angleMin);
720 u::FillH1(l, ic++, dPhi);
721 u::FillH1(l, ic++, dTheta);
722
723 if (phiP>=81. && phiP<=90.) {
724 u::FillH1(l, 14, dPhi);
725 } else if(phiP> 90. && phiP<=99.) {
726 u::FillH1(l, 15, dPhi);
727 }
728 // Unfinished - have to get absid of digit with max energy
729 u::FillH1(l, 16, dTheta);
730 u::FillH1(l, 17, dTheta);
731 if(nt) {
732 nt->Fill(dPhi,dTheta,phiP,thetaP);
733 /*
734 tv__tree = (TTree *) gROOT->FindObject("angle");
735 tv__tree->Draw("dPhi:phiP","abs(dPhi)<0.5","", 9182, 0);
736 tv__tree->Draw("dTheta:thetaP","abs(dTheta)<0.5","", 9182, 0);
737 */
738 }
739}
740
741void AliEMCALRecPointsQaESDSelector::FillHistsForShowerProfile
742(TList *l, AliEMCALRecPoint *rp, AliEMCALCellInfo * t)
743{
744 // Aug 1, 2007
745 if(l==0 || rp==0 || t==0) return;
746 // --
747 static Double_t xmean=0., xlog=0.;
748 static cellInfo rMax;
749 static Int_t phiSize;
750
751 if(rp->GetPointEnergy() < 1.0) return;
752 //if(rp->GetPointEnergy() > 8.0) return; // discard merged clusters for pi0 case
753
754 EvalLocalPhiPosition(1., rp, t, xmean, phiSize, rMax);
755 if(phiSize == 0) return; // just one row in cell directions
756
757 EvalLocalPhiPosition(5.5, rp, t, xlog, phiSize, rMax);
758 if(rMax.iPhi>1.5&&rMax.iPhi<21.5 && rMax.iEta>1.5&&rMax.iEta<46.5) {
759 u::FillH1(l, 0, xmean);
760 u::FillH1(l, 1, xlog);
761 u::FillH2(l, 2, xlog, xmean);
762 // Select two central modules
763 if((rMax.iPhim==5 || rMax.iPhim==6)){
764 // Transition to system of cell with max energy
765 xmean -= (double(rMax.iPhi)+0.5);
766 xlog -= (double(rMax.iPhi)+0.5);
767 if(rMax.iEtam>=2 && rMax.iEtam<=12){ // approximatively first half on eta
768 u::FillH2(l, 3, xlog, xmean);
769 } else {// approximatively second half on eta
770 u::FillH2(l, 4, xlog, xmean);
771 }
772 }
773 }
774}
775
776void AliEMCALRecPointsQaESDSelector::EvalLocalPhiPosition(const Double_t wlog, const AliEMCALRecPoint *rp, const AliEMCALCellInfo* t, Double_t &xcog, Int_t &phiSize, cellInfo &rMax)
777{
778 // wlog = 1 - usual center of gravity; >1 - with logarithmic weight.
779 // digits - array of digits
780 // t -
781 //==
782 // xcog - position of center of gravity in cell unit
783 if(wlog<1.0 || rp==0 || t==0) return;
784 xcog = 0;
785
786 static Double_t wtot=0., w=0., edigi=0., e=0., edigiMax=0.;
787 static Int_t absid = 0, phiMin=0, phiMax=0;
788 static cellInfo* r=0;
789
790 e = rp->GetPointEnergy();
791 wtot = 0.0;
792 edigiMax=0.;
793
794 phiMin = 23; phiMax = 0;
795 for(Int_t iDigit=0; iDigit<rp->GetMultiplicity(); iDigit++) {
796 absid = rp->GetAbsId()[iDigit];
797 edigi = rp->GetEnergiesList()[iDigit];
798 if(wlog > 1.0) w = TMath::Max( 0., wlog + TMath::Log(edigi/e));
799 else w = edigi; // just energy
800
801 r = t->GetTable(absid);
802 xcog += w*(Double_t(r->iPhi) + 0.5);
803 wtot += w;
804 if(edigi > edigiMax) {
805 edigiMax = edigi;
806 rMax = (*r);
807 }
808 if(phiMin > r->iPhi) phiMin = r->iPhi;
809 if(phiMax < r->iPhi) phiMax = r->iPhi;
810 }
811 xcog /= wtot;
812 phiSize = phiMax - phiMin;
813 // printf("phiSize %i \n", phiSize);
814}
16d3c94d 815/* unused now
816TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfTowers(const char *name)
817{
818 //
819 // ESD: towers information was saved to pseudo clusters
820 //
821 gROOT->cd();
822 TH1::AddDirectory(1);
823 new TH1F("00_EmcalMultiplicity", "multiplicity of EMCAL ", 201, -0.5, 200.5); // number of pseudo RP
824
825 new TH1F("01_EnergyOf", "energy of ", 1000, 0.0, 100.);
826
827 TList *l = u::MoveHistsToList(Form("ListOfHists%s",name, " - ESD"), kFALSE);
828 u::AddToNameAndTitleToList(l, name, name);
829
830 return l;
831}
832*/
833
834void AliEMCALRecPointsQaESDSelector::FitEffMassHist()
835{
836 TH1* h = (TH1*)fLofHistsRP->At(8);
837 AliEMCALCell::FitHist(h, GetName(), "draw");
838}
839
840void AliEMCALRecPointsQaESDSelector::PrintInfo()
841{
842 // Service routine
0fc11500 843 printf("\n %i Entrie(s) | Option(s) |%s| \n", GetOptsArray().GetEntries(), fRunOpts.Data());
844 for(int i=0; i<nAnaOpt; i++) {
845 if(GetKeyOptsValue(i)) printf(" %i |%s| \n", i, anaOpt[i]);
846 }
847
16d3c94d 848 TList *l[2] = {fLofHistsPC, fLofHistsRP};
849 printf("\n");
850 for(int i=0; i<2; i++){
851 TH1F *h = (TH1F*)l[i]->At(2);
852 printf(" %s \t: %i \n", h->GetTitle(), int(h->GetEntries()));
853 }
0fc11500 854 printf(" fDistEff %f fW0 %f fSlopePhiShift %f \n", fDistEff, fW0, fSlopePhiShift);
855}
856
857void AliEMCALRecPointsQaESDSelector::SetMomentum(Double_t p)
858{ // Jul 9, 2007
859 fPmom = p;
860 // SetName(Form("%s_p_%f5.1", GetName(), p));
16d3c94d 861}
862
0fc11500 863
16d3c94d 864AliEMCALFolder* AliEMCALRecPointsQaESDSelector::CreateEmcalFolder(const Int_t it)
865{
866 AliEMCALFolder* newFolder = new AliEMCALFolder(it); // folder for iteration #it
867 if(it>1) {
868 fEMCALOld = fEMCAL;
869 AliEMCALCalibCoefs* tabOldOut = fEMCALOld->GetCCOut();
870 AliEMCALCalibCoefs* tabNewIn = new AliEMCALCalibCoefs(*tabOldOut);
871 tabNewIn->SetName(AliEMCALFolder::fgkCCinName.Data());
872 newFolder->Add(tabNewIn);
873 }
874 fEmcalPool->Add(newFolder);
875 fEMCAL = newFolder;
876
877 return fEMCAL;
878}
879
880AliEMCALFolder* AliEMCALRecPointsQaESDSelector::GetEmcalOldFolder(const Int_t nsm)
881{
882 AliEMCALFolder* folder=0;
0fc11500 883 if(fEmcalPool) folder = (AliEMCALFolder*)fEmcalPool->FindObject(Form("EMCAL_%2.2i",nsm));
16d3c94d 884 return folder;
885}
886
887
888void AliEMCALRecPointsQaESDSelector::SetEmcalFolder(AliEMCALFolder* folder)
889{
890 fEMCAL = folder;
891 fEmcalPool->Add(fEMCAL);
892}
893
894void AliEMCALRecPointsQaESDSelector::SetEmcalOldFolder(AliEMCALFolder* folder)
895{
896 fEMCALOld = folder;
897 fEmcalPool->Add(fEMCALOld);
898}
899
900void AliEMCALRecPointsQaESDSelector::Browse(TBrowser* b)
901{
902 if(fESD) b->Add(fESD);
16d3c94d 903 if(fChain) b->Add(fChain);
904 if(fEmcalPool) b->Add(fEmcalPool);
0fc11500 905 if(geo) b->Add(geo);
906 if(fCellsInfo) b->Add(fCellsInfo);
907 //
908 if(fLofHistsPC) b->Add(fLofHistsPC);
909 if(fLofHistsRP) b->Add(fLofHistsRP);
910 if(fLKineVsRP) b->Add(fLKineVsRP);
911 if(fLShowerProfile) b->Add(fLShowerProfile);
16d3c94d 912 // if(u) b->Add(u);
913}
914
915Bool_t AliEMCALRecPointsQaESDSelector::IsFolder() const
916{
917 if(fESD || fLofHistsRP || fEmcalPool) return kTRUE;
918 return kFALSE;
919}
920
0fc11500 921void AliEMCALRecPointsQaESDSelector::Save(Int_t ver, const char *optIO)
922{ // Aug 3, 2007
923 TString dir("/home/pavlinov/ALICE/SHISHKEBAB/RF/CALIB/"); // Root directory for saving
924 TString nf=dir;
925 if(GetKeyOptsValue(kPROF)) {
926 nf += "PROF/PROFILE_";
927 nf += ver;
928 nf += ".root";
929 TFile f(nf.Data(), optIO);
930 if(f.IsOpen()) {
931 this->Write();
932 f.Close();
933 printf("<I> Save selectort to file |%s| : optIO %s \n",nf.Data(), optIO);
934 } else {
935 printf("<W> File %s exits ! Increase version : ver %i !\n", nf.Data(), ver);
936 }
937 } else {
938 printf("No PROF option -> no saving now !\n");
939 }
940}
941
942AliEMCALRecPointsQaESDSelector* AliEMCALRecPointsQaESDSelector::ReadSelector(const char* nf)
943{
944 AliEMCALRecPointsQaESDSelector* selector=0;
945
946 TH1::AddDirectory(0);
947 TFile f(nf,"READ");
948 if(f.IsOpen()) {
949 TObject *o = f.Get("AliEMCALRecPointsQaESDSelector");
950 if(o) selector = dynamic_cast<AliEMCALRecPointsQaESDSelector *>(o);
951 }
952 printf("<I> read selector %p : file |%s| \n", selector, nf);
953 return selector;
954}
955
16d3c94d 956void AliEMCALRecPointsQaESDSelector::ReadAllEmcalFolders()
957{
958 if(fEmcalPool==0) {
0fc11500 959 fEmcalPool = new TFolder("PoolOfEMCAL","");
16d3c94d 960 for(Int_t it=1; it<=10; it++){
961 AliEMCALFolder* fold = AliEMCALFolder::Read(Form("EMCALFOLDER_It%i_fit.root",it), "READ");
962 if(fold) fEmcalPool->Add(fold);
963 }
964 }
965}
966
967void AliEMCALRecPointsQaESDSelector::PictVsIterNumber(const Int_t ind, const Int_t nsm)
968{
969 // Jun 27, 2007 - unfinished; which picture is the best
970 if(ind<0 || ind>5) return;
971 gROOT->cd();
972 TH1::AddDirectory(1);
973
974 Int_t itMax = 10, it=0;
975 map <int, char*> indName;
976 indName[0] = "eff.mass";
977 indName[3] = "mass of #pi_{0}";
978 indName[4] = "resolution of #pi_{0}";
979 indName[5] = "chi^{2}/ndf";
980
981 TH1F *hout = new TH1F(indName[ind], indName[ind], itMax, 0.5, double(itMax)+0.5);
982
983 TH1::AddDirectory(0);
984 Double_t content, error;
0fc11500 985 TList* col = (TList*)fEmcalPool->GetListOfFolders();
986 for(Int_t i=0; i<col->GetSize(); i++) { // cycle on EMCAL folders
987 AliEMCALFolder* folder = static_cast<AliEMCALFolder*>(col->At(i));
16d3c94d 988 if(folder==0) continue;
989 it = folder->GetIterationNumber();
990
991 AliEMCALSuperModule* sm = folder->GetSuperModule(nsm);
992 if(sm==0) continue;
993
994 TList* l = sm->GetHists();
995 if(l==0) continue;
996
997 TH1F *hin = (TH1F*)l->At(ind);
998 if(hin==0) continue;
999
1000 if(ind !=0 ) {
1001 content = hin->GetMean();
1002 error = hin->GetRMS();
1003 } else {
1004 sm->FitEffMassHist();
1005 TF1 *f = (TF1*)hin->GetListOfFunctions()->At(0);
1006 content = error = -1.;
1007 if(f) {
1008 // content = f->GetParameter(1);
1009 //error = f->GetParameter(2);
1010 content = f->GetParameter(2);
1011 error = f->GetParError(2);
1012 }
1013 }
1014
1015 if(content > 0.0) {
1016 hout->SetBinContent(it, content);
1017 hout->SetBinError(it, error);
1018 printf(" it %i content %f +/- %f \n", it, content, error);
1019 }
1020 }
1021
1022 u::DrawHist(hout,2);
1023 hout->SetMinimum(0.0);
1024
1025}
0fc11500 1026
1027TH1F* AliEMCALRecPointsQaESDSelector::FitHistOfRecPointEnergy(const char *opt)
1028{
1029 TH1::AddDirectory(0);
1030
1031 TString OPT(opt);
1032 OPT.ToUpper();
1033
1034 Int_t ind = 4, ind2 = 16;
1035 if(GetKeyOptsValue(kIDEAL)) {
1036 ind = 16; // Jul 12, 2007
1037 ind2 = 4;
1038 } else if(GetKeyOptsValue(kCORR1)) {
1039 ind = 18;
1040 ind2 = 4;
1041 }
1042 TH1F *hold = (TH1F*)fLofHistsRP->At(ind), *h=0;
1043 if(hold == 0) return 0;
1044 if(hold->GetEntries() <10.) return hold;
1045
1046 if(OPT.Contains("CLONE")) {
1047 TString newName(Form("C_%s",hold->GetName()));
1048 h = (TH1F*)hold->Clone(newName.Data());
1049 printf(" Clone hist %s -> |%s|%s| \n",hold->GetName(),h->GetName(),h->GetTitle());
1050 } else {
1051 h = hold;
1052 }
1053 TH1F* h2 = (TH1F*)fLofHistsRP->At(ind2);
1054
1055 Double_t xmax = h->GetXaxis()->GetXmax(), xmin = 0.4*xmax;;
1056 TF1 *g = u::Gausi("fRecPointE", xmin, xmax, h);
1057 g->SetLineColor(kRed);
1058 gStyle->SetOptFit(111);
1059
1060 h->Fit(g,"N","", xmin, xmax);
1061 printf(" (1) xmin %f : xmax %f \n", xmin, xmax);
1062
1063 xmin = g->GetParameter(1) - 4.*g->GetParameter(2);
1064 xmax = g->GetParameter(1) + 4.*g->GetParameter(2);
1065 h->Fit(g,"Q+","", xmin, xmax);
1066 u::DrawHist(h2,1, 1, "same",2);
1067 printf(" (2) xmin %f : xmax %f \n", xmin, xmax);
1068
1069 return h;
1070}
1071
1072TCanvas *AliEMCALRecPointsQaESDSelector::Linearity(TList *l, int ifun)
1073{ //Jul 10, 2007
1074 if(l==0) {
1075 printf("<E> AliEMCALRecPointsQaESDSelector::Linearity :TList is zero ! Bye ! \n");
1076 return 0;
1077 }
1078 Double_t p[9]={0.5, 1., 3., 5., 10., 20., 30., 50., 100.}, ep[9];
1079 // see macro rHits.C ->defineSampleFraction
1080 TH1F* hErecOverEin = new TH1F("hErecOverEin","Ratio E_{rec}/E_{#gamma} vs E_{#gamma}", 101,0.5,101.5);
1081 // Fill hist
1082 TArrayD invRat(9), eInvRat(9), erec(9),derec(9), residual(9);
1083 TArrayD res(9), eres(9); // resolution
1084 Int_t np=9;
1085 for(Int_t var=0; var<9; var++){
1086 TH1F *h = (TH1F*)l->At(var);
1087 TF1 *f = (TF1*)h->GetListOfFunctions()->At(0);
1088 Double_t mean = f->GetParameter(1);
1089 Double_t emean = f->GetParError(1);
1090
1091 Int_t bin = hErecOverEin->FindBin(p[var]);
1092 hErecOverEin->SetBinContent(bin, mean/p[var]);
1093 hErecOverEin->SetBinError(bin, emean/p[var]);
1094 //
1095 invRat[var] = p[var]/mean;
1096 eInvRat[var] = emean/p[var]*invRat[var];
1097 erec[var] = mean;
1098 derec[var] = emean;
1099 // Resolution in %
1100 res[var] = 100.*f->GetParameter(2)/p[var];
1101 eres[var] = 100.*f->GetParError(2)/p[var];
1102 ep[var] = 0.0;
1103 }
1104
1105 TCanvas *c = new TCanvas("Linearity","Linearity", 20,20, 700, 500);
1106 gStyle->SetOptStat(0);
1107 if(0) { // E(rec) / E(gamma)
1108 u::DrawHist(hErecOverEin,2);
1109 hErecOverEin->SetMinimum(0.984);
1110 hErecOverEin->SetMaximum(1.001);
1111 }
1112 Int_t markerColor=1;
1113 char *fun="", *optFit="";
1114 TF1 *f = 0;
1115 if(0) {
1116 if(ifun==-5) {
1117 fun = "fun5"; optFit="R+";
1118 f = new TF1(fun,"([0]*(x-7.5)*(x-7.5))*exp([1]+[2]*x+[3]*x*x)+1.", 0.0, 100.);
1119 f->FixParameter(0, 1.95380e-05);
1120 // f->SetParameter(0, 1.95380e-05);
1121 // f->SetParameter(1, 1.0);
1122 } else if(ifun==-4) {
1123 fun = "fun4"; optFit="R+";
1124 f = new TF1(fun,"([0]*(x-7.5)*(x-7.5))*([1]+[2]*abs(x)+[3]*(x-[4])*(x-[4]))+1.", 0.0, 100.);
1125 f->FixParameter(0, 1.95380e-05);
1126 f->SetParameter(4, 50.);
1127 // f = new TF1(fun,"exp(([0]+[1]*x)*(x-[2])*(x-[3]))", 0.0, 100.);
1128 //f = new TF1(fun,"([0]+[1]*x)*(x-[2])*(x-[3])", 0.0, 100.);
1129 //f->SetParameter(0, 1.);
1130 //f->SetParameter(1, 1.);
1131
1132 // f->SetParameter(2, 5.);
1133 //f->SetParLimits(2, 3.,8.);
1134 //f->SetParameter(3, 10.);
1135 //f->SetParLimits(3, 9.,15.);
1136 //f->SetParameter(2, 2.e-4);
1137 } else if(ifun==-3) {
1138 fun = "fun3"; optFit="R+";
1139 f = new TF1(fun,"[0]+[1]*x+[2]*x*x+[3]*x*x*x", 0.0, 11.);
1140 } else if(ifun==-2) {
1141 fun = "fun2"; optFit="R+";
1142
1143 /*
1144 f = new TF1(fun,"[0]*(x-7.5)+[1]*(x-7.5)*(x-7.5)+[2]", 0.0, 10.1);
1145 f->SetParameter(0, 5.98727e-04);
1146 f->SetParameter(1, 2.12045e-04);
1147 f->SetParameter(2, 1.);
1148 */
1149 f = new TF1(fun,"[0]+[1]*x+[2]*x*x", 9.0, 100.1);
1150 } else if(ifun==-1) {
1151 fun = "fun0"; optFit="R+";
1152 f = new TF1(fun,"[0]", 0.0, 100.1);
1153 } else if(ifun>=2) {
1154 fun = Form("pol%i",ifun);
1155 optFit="+";
1156 }
1157 TGraphErrors *gr = u::DrawGraphErrors(np, erec.GetArray(),invRat.GetArray(),derec.GetArray(),eInvRat.GetArray(),
1158 // markerColor,21+markerColor,"AP", " Ratio E_{#gamma}/E_{Rec} ", "E_{Rec} "," E_{#gamma}/E_{Rec}",
1159 markerColor,21+markerColor,"AP", " Ratio E_{#gamma}/E_{Rec}^{corr} ", "E_{Rec}^{corr} "," E_{#gamma}/E_{Rec}^{corr}",
1160 ifun, optFit, fun);
1161 gr->GetHistogram()->SetAxisRange(0.0,100.);
1162 // double xmi=0.999, xma=1.017;
1163 double xmi=0.995, xma=1.005;
1164 gr->GetHistogram()->SetMaximum(xma);
1165 gr->GetHistogram()->SetMinimum(xmi);
1166 gr->GetHistogram()->SetTitleOffset(1.4,"y");
1167 if(ifun==0) {
1168 f = new TF1("fres", "AliEMCALHistoUtilities::EnergyCorrectionForGamma_1(x)", 0., 101.);
1169 f->Draw("same");
1170 }
1171 }
1172 TLine *line = new TLine(0.0,1.0, 100.,1.0);
1173 line->Draw();
1174
1175 if(0) {
1176 c->Clear();
1177 for(int i=0; i<9; i++) {
1178 residual[i] = 100.*(invRat[i] - u::GetCorrectionCoefficientForGamma_1(erec[i])); // in percent
1179 printf(" erec %f : residual %5.3f \n", erec[i], residual[i]);
1180 }
1181 markerColor = 2;
1182 TGraphErrors *gr2=u::DrawGraphErrors(np, erec.GetArray(),residual.GetArray(),derec.GetArray(),eInvRat.GetArray(),
1183 markerColor,21+markerColor,"AP"," residual in %, rec.point level", "E_{Rec} "," residual in %",
1184 -1, "", 0);
1185 gr2->GetHistogram()->SetAxisRange(0.0,100.);
1186 gr2->GetHistogram()->SetMaximum(0.2);
1187 gr2->GetHistogram()->SetMinimum(-0.1);
1188 line = new TLine(0.0,0.0, 101.,0.0);
1189 line->Draw();
1190 TLatex *lat = u::Lat("linearity better 0.2% after correction",20., 0.15, 12, 0.06, 1);
1191 if(lat);
1192 }
1193 if(1) {
1194 TString latexName;
1195 c->Clear();
1196 gStyle->SetOptFit(111);
1197 markerColor = 1;
1198 ifun = -11;
1199 f = u::GetResolutionFunction("FRES2", latexName);
1200 f->SetNpx(10000);
1201 f->SetLineColor(kBlack);
1202
1203 TGraphErrors *gr3=u::DrawGraphErrors(np, p,res.GetArray(), ep, eres.GetArray(),
1204 markerColor,21+markerColor,"AP"," resolution in %, rec.point level","E_{#gamma} "," resolution in %",
1205 ifun, "+", f->GetName());
1206 gr3->GetHistogram()->SetAxisRange(0.0,101.);
1207 gr3->GetHistogram()->SetMaximum(14.);
1208 gr3->GetHistogram()->SetMinimum(0.0);
1209 gr3->SetMarkerSize(1.5);
1210
1211 TLatex *lat = u::Lat(latexName.Data(),82., 11., 12, 0.06, 1);
1212 if(lat);
1213 // Exp. data
1214 TF1 *fexp = new TF1(*f);
1215 fexp->SetName("fexp");
1216 fexp->SetParameter(0, 2.168);
1217 fexp->SetParameter(1, 8.818);
1218 fexp->SetLineWidth(1);
1219 fexp->SetLineColor(kRed);
1220 fexp->SetLineStyle(1);
1221 fexp->Draw("same");
1222
1223 TLegend *leg = new TLegend(0.21,0.36, 0.68,0.82);
1224 leg->AddEntry(gr3, "MC data", "P");
1225 leg->AddEntry(f, "fit of MC data", "L");
1226 TLegendEntry *ent3 = leg->AddEntry(fexp, "#splitline{fit of exp.data}{FNAL, Nov 2005}", "L");
1227 ent3->SetTextColor(fexp->GetLineColor());
1228 leg->Draw();
1229 }
1230 c->Update();
1231
1232 return c;
1233}
1234
1235TCanvas *AliEMCALRecPointsQaESDSelector::DrawKineVsRP(TList *l)
1236{ //Jul 25, 2007
1237 if(l==0) {
1238 printf("<W> AliEMCALRecPointsQaESDSelector::DrawKineVsRP : TList is zero ! \n");
1239 return 0;
1240 }
1241 TCanvas *c = new TCanvas("KineVsRP","KineVsRP", 20,20, 700, 500);
1242 gStyle->SetOptStat(1110);
1243 c->Divide(2,2);
1244
1245 c->cd(1);
1246 TH1F* h1 = (TH1F*)l->At(10);
1247 u::DrawHist(h1,2);
1248
1249 c->cd(2);
1250 TH1F* h2 = (TH1F*)l->At(11);
1251 u::DrawHist(h2,2);
1252
1253 c->cd(3);
1254 TH1F* h3 = (TH1F*)l->At(12);
1255 u::DrawHist(h3,2);
1256 u::DrawHist((TH1F*)l->At(14), 1,kRed,"same");
1257 u::DrawHist((TH1F*)l->At(15), 1,kGreen,"same");
1258
1259 c->cd(4);
1260 TH1F* h4 = (TH1F*)l->At(13);
1261 u::DrawHist(h4, 2);
1262 u::DrawHist((TH1F*)l->At(16), 1,kRed,"same");
1263 u::DrawHist((TH1F*)l->At(17), 1,kGreen,"same");
1264
1265 /*
1266 TH1F* h2 = (TH1F*)l->At(11);
1267 u::DrawHist(h2,2);
1268 */
1269
1270 c->Update();
1271 return c;
1272}
1273
1274TCanvas* AliEMCALRecPointsQaESDSelector::DrawMeanVsLog(TH2F *h2)
1275{
1276 // h - input histogramm : mean vds log coordinates
1277 if(h2==0) return 0;
1278
1279 TCanvas *c = new TCanvas("KineVsRP","KineVsRP", 20,20, 700, 500);
1280
1281 TH1F *hid1 = new TH1F("h1","h1 title ", h2->GetNbinsX(),
1282 h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
1283
1284 gROOT->cd();
1285 TH1::AddDirectory(1);
1286 TString newName;
1287 for(int ix=1; ix<=h2->GetNbinsX();ix++) {
1288 newName = "hd_"; newName += ix;
1289 TH1D *hd = h2->ProjectionY(newName.Data(),ix,ix);
1290 if(hd->Integral()>=4.) {
1291 // lhd->Add(hd);
1292 hid1->SetBinContent(ix, hd->GetMean());
1293 hid1->SetBinError(ix, hd->GetRMS());
1294 }
1295 }
1296 TF1 *f1 = new TF1("fcg", "0.5*TMath::SinH(x/[0])/TMath::SinH(0.5/[0])", -0.5, 0.5);
1297 f1->SetParameter(0,0.13);
1298 f1->SetLineColor(kRed);
1299 f1->SetLineWidth(1);
1300
1301 gStyle->SetOptStat(0);
1302 gStyle->SetOptFit(111);
1303 hid1->Fit(f1,"R+");
1304
1305 u::DrawHist(hid1,2);
1306
1307//u::DrawHist(h2,2);
1308
1309 c->Update();
1310 return c;
1311}
1312
1313TCanvas* AliEMCALRecPointsQaESDSelector::DrawPhiEtaAnglesDistribution(const char* gn)
1314{
1315 // Aug 6, 2007
1316 // Proper geometry should be defined already
1317 // dTheta distibution has two peaks which coresponding different cells inside module !
1318
1319 TCanvas *c = new TCanvas("Geometry","Geometry", 20,20, 700, 500);
1320 c->Divide(2,2);
1321
1322 if(geo==0) geo = AliEMCALGeometry::GetInstance(gn);
1323
1324 gROOT->cd();
1325 TH1::AddDirectory(1);
1326 TH1F *hDtheta = new TH1F("hDtheta","#Delta#theta in one SM", 60, -2.0, +1.0); // in degree
1327 TH2F *hDtheta2 = new TH2F("hDtheta2","#Delta#theta vs iEta of cell", 48, -0.5, 47.5, 60,-2.0,+1.0);
1328
1329 TH1F *hDphi = new TH1F("hDphi","#Delta#ph in one SM", 2000, -10.0, +10.0); // in degree
1330
1331 TVector3 vg3;
1332 AliEMCALCellInfo* t = GetCellsInfo();
1333 Double_t thetaCell=0., thetaModule=0.;
1334 Double_t phiCell=0., dphi=0.;
1335
1336 for(int absid=0; absid<12*24*4; absid++){
1337 cellInfo *r = t->GetTable(absid);
1338 geo->GetGlobal(absid, vg3);
1339
1340 thetaCell = vg3.Theta()*TMath::RadToDeg();
1341 thetaModule = 90. - 1.5*r->iEtam;
1342 hDtheta->Fill(thetaCell - thetaModule);
1343 hDtheta2->Fill(double(r->iEta), thetaCell - thetaModule);
1344
1345 phiCell = vg3.Phi()*TMath::RadToDeg();
1346 dphi = phiCell - 90.;
1347 hDphi->Fill(dphi);
1348 }
1349
1350 c->cd(1);
1351 u::DrawHist(hDtheta,2);
1352 c->cd(2);
1353 u::DrawHist(hDtheta2,2);
1354 c->cd(3);
1355 u::DrawHist(hDphi,2);
1356
1357 c->Update();
1358 return c;
1359}
1360
1361TCanvas* AliEMCALRecPointsQaESDSelector::DrawDeffVsEnergy()
1362{
1363 // Aug 31, 2007 - obsolete
1364 Double_t p[]={0.5, 1.0, 3., 5., 10., 20.,30.,50.,100.};
1365 Double_t ep[]={0., 0., 0., 0., 0., 0., 0., 0., 0.};
1366 // 2 pars
1367 Double_t deff[]={9.07515, 10.0835, 11.2032, 11.4229, 12.3578, 13.0332, 13.3281, 13.7910, 14.3319};
1368 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};
1369 //
1370 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};
1371 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};
1372 int np=sizeof(p)/sizeof(Double_t);
1373 printf("<I> AliEMCALRecPointsQaESDSelector::DrawDeffVsEnergy() | np %i \n", np);
1374
1375 TCanvas *c = new TCanvas("Deff","Deff", 20,20, 700, 500);
1376 c->Divide(2,1);
1377
1378 Int_t markerColor=1;
1379 c->cd(1);
1380 TF1 *fdeff= new TF1("fdeff","[0]+[1]*log(x)",0.4, 100.4);
1381 if(fdeff);
1382 TGraphErrors *gr = u::DrawGraphErrors(np, p,deff, ep, edeff,
1383 markerColor,21+markerColor,"AP"," D_{eff} vs E_{#gamma} ","E_{#gamma} "," D_{eff} in cm ",
1384 -1, "", 0);
1385 // -1, "+", fdeff->GetName());
1386 gr->GetHistogram()->SetMaximum(15.);
1387 gPad->SetLogx(1);
1388
1389 c->cd(2);
1390 TGraphErrors *gr2 = u::DrawGraphErrors(np, p,w0, ep, ew0,
1391 markerColor,21+markerColor,"AP"," w_{0} vs E_{#gamma} ","E_{#gamma} "," w_{0} ",
1392 -1, "", 0);
1393
1394 gr2->GetHistogram()->SetMaximum(5.);
1395 gPad->SetLogx(1);
1396
1397 c->Update();
1398 return c;
1399}
1400
1401TCanvas* AliEMCALRecPointsQaESDSelector::DrawDeffVsEnergy2(const char *opt)
1402{
1403 // Aug 28, 2008 - read pars and pars errors from files
1404 Double_t p[]={0.5, 1.0, 3., 5., 10., 20.,30.,50.,100.};
1405 Double_t ep[]={0., 0., 0., 0., 0., 0., 0., 0., 0.};
1406 // 2 pars
1407 Double_t deff[9], edeff[9], w0[9], ew0[9]; // max size now
1408 TString OPT(opt);
1409
1410 int np = sizeof(p)/sizeof(Double_t);
1411 printf("<I> AliEMCALRecPointsQaESDSelector::DrawDeffVsEnergy2() | np %i \n", np);
1412 ReadParsDeffAndW0("/data/r22b/ALICE/CALIB/FIT/", deff, edeff, w0, ew0, 1);
1413
1414 TCanvas *c = new TCanvas("Deff","Deff", 20,20, 700, 500);
1415 c->Divide(2,1);
1416
1417 TF1 *fdeff = 0, *fw0 = 0;
1418 TString optFit(""), funName("");
1419 if(OPT.Contains("fit1")) {
1420 fdeff= new TF1("fdeff","[0]+[1]*log(x)",0.1, 101.);
1421 fdeff->SetLineColor(kRed);
1422 fdeff->SetLineWidth(1);
1423
1424 /* good description - "[0]/(1.+exp([1]*x))"
1425 1 p0 4.82208e+00 9.93617e-04 -0.00000e+00 7.16648e-06
1426 2 p1 -7.79655e-01 4.07500e-03 -0.00000e+00 2.01009e-03
1427 better description - "[0]/(1.+exp([1]*(x-[2])))" (like the Woods-Saxon potential)
1428 1 p0 4.83713e+00 1.00437e-03 -6.98984e-10 4.90371e-07
1429 2 p1 -2.77970e-01 3.35587e-03 -1.79239e-09 2.67312e-07
1430 3 p2 4.41116e+00 8.72191e-02 1.82791e-04 1.55643e-05
1431 */
1432 fw0= new TF1("fw0","[0]/(1.+exp([1]*(x+[2])))",0.1, 101.);
1433 fw0->SetLineColor(kRed);
1434 fw0->SetLineWidth(1);
1435 fw0->SetParameter(0, 4.8);
1436 fw0->SetParameter(1, -2.77970e-01);
1437 fw0->SetParameter(2, 4.41116);
1438 optFit = "+";
1439 }
1440 if(fdeff) funName = fdeff->GetName();
1441
1442 Int_t markerColor=1;
1443 c->cd(1);
1444 gStyle->SetOptFit(111);
1445 TGraphErrors *gr = u::DrawGraphErrors(np, p,deff, ep, edeff,
1446 markerColor,21+markerColor,"AP"," D_{eff} vs E_{#gamma} ","E_{#gamma} "," D_{eff} in cm ",
1447 -1, optFit.Data(), funName.Data());
1448 gr->GetHistogram()->SetMaximum(15.);
1449 gPad->SetLogx(1);
1450 TLegend *leg1 = new TLegend(0.12,0.76, 0.70,0.90);
1451 TLegendEntry *le1 = leg1->AddEntry(fdeff, Form("%s",fdeff->GetTitle()), "lp");
1452 //TLegendEntry *le1 = leg1->AddEntry(fdeff, Form("%4.2f+%4.2f*log(E_{#gamma})",
1453 //fdeff->GetParameter(0),fdeff->GetParameter(1)), "lp");
1454 le1->SetTextColor(fdeff->GetLineColor());
1455 leg1->Draw();
1456
1457 c->cd(2);
1458 gStyle->SetOptFit(111);
1459 if(fw0) funName = fw0->GetName();
1460 TGraphErrors *gr2 = u::DrawGraphErrors(np, p,w0, ep, ew0,
1461 markerColor,21+markerColor,"AP"," w_{0} vs E_{#gamma} ","E_{#gamma} "," w_{0} ",
1462 -1, optFit.Data(), funName.Data());
1463
1464 gr2->GetHistogram()->SetMaximum(5.);
1465
1466 TLegend *leg2 = new TLegend(0.17,0.6, 0.99,0.72);
1467 TLegendEntry *le2 = leg2->AddEntry(fw0, Form("%s",fw0->GetTitle()), "lp");
1468 //TLegendEntry *le2 = leg2->AddEntry(fw0, Form("#frac{%4.2f}{1.+exp(%4.2f*(x+%4.2f)}",
1469 //fw0->GetParameter(0),fw0->GetParameter(1),fw0->GetParameter(2)), "lp");
1470 le2->SetTextColor(fw0->GetLineColor());
1471 leg2->Draw();
1472 //gPad->SetLogx(1);
1473
1474 c->Update();
1475 return c;
1476}
1477
1478void AliEMCALRecPointsQaESDSelector::ReadParsDeffAndW0
1479(const char *dirName, double *deff, double *edeff, double *w0, double *ew0, const Int_t pri)
1480{
1481 int strategy = 0, itmp=0;
1482 char line[100];
1483 for(int var=11; var<=19; var++){
1484 int ind = var -11;
1485 ifstream fin;
1486 TString fname = Form("%s/fitVar%iStrategy%i.txt", dirName, var, strategy);
1487 //printf(" open file %s \n", fname.Data());
1488 fin.open(fname.Data());
1489
1490 for(int i=1; i<=2; i++) { // skip to lines
1491 fin.getline(line,100).eof();
1492 if(pri>=2) printf("%s \n", line);
1493 }
1494
1495 fin >> itmp >> deff[ind] >> edeff[ind];
1496 fin >> itmp >> w0[ind] >> ew0[ind];
1497 if(pri>=1)
1498 printf(" %i def %f+/-%f : def %f+/-%f\n", ind, deff[ind],edeff[ind],w0[ind],ew0[ind]);
1499 fin.getline(line,100).eof(); // skip last line
1500 fin.close();
1501 }
1502}
1503
1504TCanvas* AliEMCALRecPointsQaESDSelector::DrawSpaceResolution()
1505{
1506 // Sep 4, 2007;
1507 // Space resolution vs optimised space resolution
1508 Double_t p[]={0.5, 1.0, 3., 5., 10., 20.,30.,50.,100.};
1509 Double_t ep[]={0., 0., 0., 0., 0., 0., 0., 0., 0.};
1510 Double_t spAng[]={0.1877, 0.1351, 0.08144, 0.06331, 0.05384, 0.04876, 0.04706, 0.04656, 0.04726};
1511 Double_t rmsSpAng[]={0.1033, 0.07685, 0.05235, 0.03992, 0.04012, 0.04257, 0.03472, 0.02814, 0.02784};
1512 Double_t spAngOpt[]={0.1733, 0.1311, 0.07961, 0.06401, 0.05347, 0.04618, 0.04288, 0.04, 0.03802};
1513 Double_t rmsSpAngOpt[]={0.09476, 0.07472, 0.05011, 0.04242, 0.04075, 0.04304, 0.03545, 0.02744, 0.02593};
1514
1515 int np=sizeof(p)/sizeof(Double_t);
1516 printf("<I> AliEMCALRecPointsQaESDSelector::DrawSpaceResolution() | np %i \n", np);
1517
1518 Double_t* eSpAng = new Double_t[np];
1519 Double_t* eSpAngOpt = new Double_t[np];
1520 Double_t C=TMath::Sqrt(8000.), C2 = 1000.*TMath::DegToRad();
1521 for(int i=0; i<np; i++){
1522 spAng[i] *= C2;
1523 rmsSpAng[i] *= C2;
1524 spAngOpt[i] *= C2;
1525 rmsSpAngOpt[i] *= C2;
1526
1527 eSpAng[i] = spAng[i]/C;
1528 eSpAngOpt[i] = spAngOpt[i]/C;
1529 }
1530
1531 TCanvas *c = new TCanvas("Deff","Deff", 20,20, 700, 500);
1532 //c->Divide(2,1);
1533
1534 c->cd(1);
1535 gStyle->SetOptFit(111);
1536
1537 TGraphErrors *gr = u::DrawGraphErrors(np, p,spAng, ep, eSpAng,
1538 kBlack, 21,"AP","Angle resolution (mrad) vs E_{#gamma} ","E_{#gamma} "," anlgle resolution (mrad) ",
1539 -1, "", 0);
1540 gr->GetHistogram()->SetMaximum(4.);
1541 gr->GetHistogram()->SetMinimum(0.6);
1542 gPad->SetLogx(1);
1543 gPad->SetLogy(1);
1544
1545 TF1 *fang = new TF1("fang","[0]+[1]/sqrt(x)",0.1, 101.);
1546 fang->SetLineColor(kRed);
1547 fang->SetLineWidth(1);
1548
1549 TGraphErrors *gr2 = u::DrawGraphErrors(np, p,spAngOpt, ep, eSpAngOpt,
1550 kRed,22,"P"," space vs E_{#gamma} ","E_{#gamma} "," anlgle resolution (mrad) ",
1551 -1, "+", fang->GetName());
1552 TLegend *leg1 = new TLegend(0.40,0.57, 0.92,0.87);
1553 TLegendEntry *le1 = leg1->AddEntry(gr, "Initial angle resolution", "p");
1554 le1->SetTextColor(gr->GetLineColor());
1555
1556 TLegendEntry *le2 = leg1->AddEntry(gr2, "Optimized angle resolution", "p");
1557 le2->SetTextColor(gr2->GetLineColor());
1558
1559 TLegendEntry *le3 = leg1->AddEntry(fang,
1560 Form("%3.2f+#frac{%4.2f}{#sqrt{E}}",fang->GetParameter(0),fang->GetParameter(1)), "lp");
1561 le3->SetTextColor(fang->GetLineColor());
1562
1563 leg1->Draw();
1564
1565 c->Update();
1566
1567 return c;
1568}
1569
1570void AliEMCALRecPointsQaESDSelector::ResetAllListOfHists()
1571{
1572 u::ResetListOfHists(fLofHistsPC);
1573 u::ResetListOfHists(fLofHistsRP);
1574 u::ResetListOfHists(fLKineVsRP);
1575 u::ResetListOfHists(fLShowerProfile);
1576}
1577
1578void AliEMCALRecPointsQaESDSelector::ReloadChain(Long64_t entry)
1579{ // unused now
1580 if(fChain) {
1581 fChain->LoadTree(entry);
1582 }
1583}
1584
1585void AliEMCALRecPointsQaESDSelector::GetInitialParsForFit(const Int_t var, Double_t &deff, Double_t &w0, Double_t &phislope, const int phiCase)
1586{
1587 // int phiCase=0; // 0 or 1
1588 phislope = 0.001;
1589 if(var==11) { // stay here
1590 deff = 9.07515;
1591 w0 = 3.44679;
1592 } else if(var==12) {
1593 deff = 1.00835e+01;
1594 w0 = 3.82714e+00;
1595 } else if(var==13) {
1596 deff = 1.12032e+01;
1597 w0 = 4.15035e+00;
1598 } else if(var==14) {
1599 deff = 1.14229e+01;
1600 w0 = 4.36650;
1601 } else if(var==15) {
1602 deff = 1.21361e+01;
1603 w0 = 4.72875;
1604 } else if(var==16) {
1605 deff = 1.28096e+01;
1606 w0 = 4.81753e+00;
1607 } else if(var==17) {
1608 deff = 1.33281e+01;
1609 w0 = 4.63289e+00;
1610 } else if(var==18) {
1611 deff = 1.37910e+01;
1612 w0 = 4.66568;
1613 } else if(var==19) {// Aug 15, 2007
1614 switch (phiCase) {
1615 case 0: // phislope was defined than deff and w0 were fixed
1616 deff = 1.43319e+01;
1617 w0 = 4.69279;
1618 phislope = 2.41419e-04;
1619 break;
1620 case 1: // all parameters free
1621 deff = 1.46327e+01;
1622 w0 = 4.68125;
1623 phislope = 8.95559e-04;
1624 break;
1625 }
1626 } else {
1627 printf("<E> var % i -> no definition ! \n", var);
1628 assert(0);
1629 }
1630 printf("<I> AliEMCALRecPointsQaESDSelector::GetInitialParForFit()\n deff %f w0 %f phislope %f \n",
1631 deff,w0, phislope);
1632}