EMCAL/AliEMCALHistoUtilities.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2002, 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 // This is a set of histogram
  20 // utilities for the EMCAL
  21 // to make some common
  22 // functions easier
  23 //
  24 //*-- Authors: J.L. Klay (LLNL) & Aleksei Pavlinov (WSU)
  25
  26 #include "AliEMCALHistoUtilities.h"
  27
  28 #include <iostream>
  29 #include <iomanip>
  30 #include <fstream>
  31
  32 #include <TROOT.h>
  33 #include <TPad.h>
  34 #include <TFile.h>
  35 #include <TH1.h>
  36 #include <TH2.h>
  37 #include <TProfile.h>
  38 #include <THnSparse.h>
  39 #include <TF1.h>
  40 #include <TGraph.h>
  41 #include <TGraphErrors.h>
  42 #include <TLatex.h>
  43 #include <TChain.h>
  44 #include <TList.h>
  45 #include <TObjArray.h>
  46 #include <TObjString.h>
  47 #include <TRegexp.h>
  48 #include <TString.h>
  49 #include <TLorentzVector.h>
  50 #include <Gtypes.h> // color, line style and so on
  51 #include <TArrayF.h>
  52
  53 //#include "AliESDCaloCluster.h"
  54 //#include "AliEMCALRecPoint.h"
  55 //#include "AliRunLoader.h"
  56 #include "AliHeader.h"
  57 #include "AliGenEventHeader.h"
  58 #include "AliGenPythiaEventHeader.h"
  59
  60 using namespace std;
  61
  62 ClassImp(AliEMCALHistoUtilities)
  63
  64 AliEMCALHistoUtilities::AliEMCALHistoUtilities(const char *name, const char *tit) : TNamed(name,tit)
  65 {
  66   // constructor
  67 }
  68
  69 AliEMCALHistoUtilities::~AliEMCALHistoUtilities()
  70 {
  71         //destructor
  72 }
  73
  74 TList* AliEMCALHistoUtilities::MoveHistsToList(const char* name, Bool_t putToBrowser, Bool_t setOwner)
  75 {
  76   // Move HIST to list
  77   gROOT->cd();
  78   TIter nextHist(gDirectory->GetList());
  79   TList *list = new TList;
  80   list->SetName(name);
  81   if(setOwner) list->SetOwner(setOwner);
  82   TObject *objHist=0;
  83   // Move from gROOT to list
  84   while((objHist=nextHist())){
  85     //objHist->Dump();
  86     if (!objHist->InheritsFrom("TH1")) continue;
  87     ((TH1*)objHist)->SetDirectory(0);
  88     list->Add(objHist);
  89   }
  90   // Clear gROOT
  91   gDirectory->GetList()->Clear();
  92
  93   if(putToBrowser) gROOT->GetListOfBrowsables()->Add((TObject*)list);
  94   return list;
  95 }
  96
  97 void AliEMCALHistoUtilities::FillH1(TList *l, Int_t ind, Double_t x, Double_t w, Double_t error)
  98 {
  99   //fill 1d histogram
 100   static TH1* hid=0;
 101   static int  bin=0;
 102   if(l == 0) return;
 103
 104   if(ind>=0 && ind < l->GetSize()){
 105     hid = dynamic_cast<TH1 *>(l->At(ind));
 106     if (hid==0) return;
 107
 108     if(error <= 0.0) { // standard way
 109       hid->Fill(x,w);
 110     } else {
 111       bin = hid->FindBin(x);
 112       hid->SetBinContent(bin,w);
 113       hid->SetBinError(bin,error);
 114     }
 115   }
 116 }
 117
 118 void AliEMCALHistoUtilities::FillH2(TList *l, Int_t ind, Double_t x, Double_t y, Double_t w)
 119 {
 120   //fill 2d histogram
 121   static TH2* hid=0;
 122   if(l == 0) return;
 123   if(ind>=0 && ind < l->GetSize()){
 124     hid = dynamic_cast<TH2 *>(l->At(ind));
 125     if(hid) hid->Fill(x,y,w);
 126   }
 127 }
 128
 129 void AliEMCALHistoUtilities::FillHProf(TList *l, Int_t ind, Double_t x, Double_t y, Double_t w)
 130 {
 131   // fill profile histogram
 132   static TProfile* h=0;
 133   if(l == 0) return;
 134   if(ind>=0 && ind < l->GetSize()){
 135     h = dynamic_cast<TProfile *>(l->At(ind));
 136     if(h) h->Fill(x,y,w);
 137   }
 138 }
 139
 140 void AliEMCALHistoUtilities:: FillHnSparse(TList *l, Int_t ind, Double_t* x, Double_t w)
 141 {
 142   // Nov 02,2010: fill THnSparse hist
 143   static THnSparse* hsp=0;
 144   if(l==0 || x==0) return;
 145   if(ind>=0 && ind < l->GetSize()){
 146     hsp = dynamic_cast<THnSparse *>(l->At(ind));
 147     if(hsp) hsp->Fill(x,w);
 148   }
 149 }
 150
 151 int AliEMCALHistoUtilities::SaveListOfHists(TList *mylist,const char* name,Bool_t kSingleKey,const char* opt)
 152 {
 153   //write histograms to file
 154   printf(" Name of out file |%s|\n", name);
 155   int save = 0;
 156   if(mylist && mylist->GetSize() && strlen(name)){
 157     TString nf(name);
 158     if(nf.Contains(".root") == kFALSE) nf += ".root";
 159     TFile file(nf.Data(),opt);
 160     TIter nextHist(mylist);
 161     TObject* objHist=0;
 162     int nh=0;
 163     if(kSingleKey) {
 164        file.cd();
 165        mylist->Write(mylist->GetName(),TObject::kSingleKey);
 166        mylist->ls();
 167        save = 1;
 168     } else {
 169       while((objHist=nextHist())) { // loop over list
 170         if(objHist->InheritsFrom("TH1")) {
 171           TH1* hid = (TH1*)objHist;
 172           file.cd();
 173           hid->Write();
 174           nh++;
 175           printf("Save hist. %s \n",hid ->GetName());
 176         }
 177       }
 178       printf("%i hists. save to file -> %s\n", nh,file.GetName());
 179       if(nh>0) save = 1;
 180     }
 181     file.Close();
 182   } else {
 183     printf("AliEMCALHistoUtilities::SaveListOfHists : N O  S A V I N G \n");
 184   }
 185   return save;
 186 }
 187
 188 void AliEMCALHistoUtilities::AddToNameAndTitle(TNamed *h, const char *name, const char *title)
 189 {
 190   // Oct 15, 2007
 191   if(h==0) return;
 192   if(name  && strlen(name))  h->SetName(Form("%s%s",h->GetName(),name));
 193   if(title && strlen(title)) h->SetTitle(Form("%s%s",h->GetTitle(),title));
 194 }
 195
 196 void AliEMCALHistoUtilities::AddToNameAndTitleToList(TList *l, const char *name, const char *title)
 197 {
 198   // Oct 15, 2007
 199   if(l==0) return;
 200   if(name || title) {
 201     for(int i=0; i<l->GetSize(); i++) {
 202       TObject *o = l->At(i);
 203       if(o->InheritsFrom("TNamed")) {
 204         TNamed *h = dynamic_cast<TNamed *>(o);
 205         AddToNameAndTitle(h, name, title);
 206       }
 207     }
 208   }
 209 }
 210
 211 void AliEMCALHistoUtilities::ResetListOfHists(TList *l)
 212 {
 213   // Oct 15, 2007
 214   if(l==0) return;
 215   for(int i=0; i<l->GetSize(); i++) {
 216     TH1F* h = (TH1F*)l->At(i);
 217     h->Reset();
 218   }
 219 }
 220
 221 void AliEMCALHistoUtilities::Titles(TH1 *hid, const char *titx,const char *tity)
 222 {
 223   if(hid) {
 224     hid->SetXTitle(titx);
 225     hid->SetYTitle(tity);
 226   } else {
 227     printf("<W> TAliasPAI::titles() -> hid is 0 !\n");
 228   }
 229 }
 230
 231 TList* AliEMCALHistoUtilities::CreateProjectionsX(TList *l, const Int_t ind, const Char_t* name)
 232 {
 233   // Sep 4 - seperate list for projections
 234   TH2F *hid = dynamic_cast<TH2F *>(l->At(ind));
 235   if(hid == 0) {
 236     printf("<E> Object at ind %i is %s : should ne TH2F \n",
 237            ind, l->At(ind)->ClassName());
 238     return 0;
 239   }
 240
 241   TList *lpt = new TList;
 242   lpt->SetName(Form("%s%s", hid->GetName(), name));
 243   l->Add(lpt);
 244
 245   TAxis*  yax   = hid->GetYaxis();
 246   TString tity  = yax->GetTitle();
 247   tity.ReplaceAll(" ","");
 248   TString name1 = hid->GetName();
 249   TString nam   = name1(3,name1.Length()-3);
 250   TString tit1  = hid->GetTitle(), tit=tit1;
 251
 252   TH1::AddDirectory(0);
 253
 254   TH1D *hd = hid->ProjectionX(Form("00_%sProx",nam.Data()), 1, yax->GetNbins());
 255   tit += Form(" : %5.2f < %s < %5.2f (GeV/c)", yax->GetXmin(),tity.Data(),yax->GetXmax());
 256   hd->SetTitle(tit.Data());
 257   lpt->Add(hd);
 258
 259   for(Int_t iy=1; iy<=yax->GetNbins(); iy++){
 260     tit = tit1;
 261     tit += Form(" : %5.2f < %s  < %5.2f (GeV/c)", yax->GetBinLowEdge(iy), tity.Data(), yax->GetBinUpEdge(iy));
 262     hd = hid->ProjectionX(Form("%2.2i_%sProx%i",iy, nam.Data(),iy),iy,iy);
 263     hd->SetTitle(tit.Data());
 264     lpt->Add(hd);
 265   }
 266
 267   return lpt;
 268 }
 269
 270 TLatex *AliEMCALHistoUtilities::Lat(const char *text, Float_t x,Float_t y, Int_t align, Float_t tsize, short tcolor)
 271 {
 272   // Oct 15, 2007
 273   double y1=y;
 274   TLatex *latex = new TLatex;
 275   latex->SetTextAlign(align);
 276   latex->SetTextSize(tsize);
 277   latex->SetTextColor(tcolor);
 278   latex->DrawLatex(x, y1, text);
 279   printf("<I> AliEMCALHistoUtilities::lat() -> %s gPad->GetLogy() %i\n",
 280   text, gPad->GetLogy());
 281   return latex;
 282 }
 283
 284 TGraph *AliEMCALHistoUtilities::DrawGraph(Int_t n, Double_t *x, Double_t *y, Int_t markerColor,
 285 Int_t markerStyle, const char* opt, const char* tit, const char* xTit,const char* yTit, Int_t ifun,
 286 const char *optFit, const char *fun)
 287 {
 288   /* Drawing options
 289   chopt='L' :  A simple polyline between every points is drawn
 290   chopt='F' :  A fill area is drawn ('CF' draw a smooth fill area)
 291   chopt='A' :  Axis are drawn around the graph
 292   chopt='C' :  A smooth Curve is drawn
 293   chopt='*' :  A Star is plotted at each point
 294   chopt='P' :  Idem with the current marker
 295   chopt='B' :  A Bar chart is drawn at each point
 296   chopt='1' :  ylow=rwymin
 297   chopt='X+' : The X-axis is drawn on the top side of the plot.
 298   chopt='Y+' : The Y-axis is drawn on the right side of the plot.
 299
 300     Fitting options
 301    The list of fit options is given in parameter option.
 302       option = "W"  Set all errors to 1
 303              = "U" Use a User specified fitting algorithm (via SetFCN)
 304              = "Q" Quiet mode (minimum printing)
 305              = "V" Verbose mode (default is between Q and V)
 306              = "B" Use this option when you want to fix one or more parameters
 307                    and the fitting function is like "gaus","expo","poln","landau".
 308              = "R" Use the Range specified in the function range
 309              = "N" Do not store the graphics function, do not draw
 310              = "0" Do not plot the result of the fit. By default the fitted function
 311                    is drawn unless the option"N" above is specified.
 312              = "+" Add this new fitted function to the list of fitted functions
 313                    (by default, any previous function is deleted)
 314              = "C" In case of linear fitting, not calculate the chisquare
 315                     (saves time)
 316              = "F" If fitting a polN, switch to minuit fitter
 317              = "ROB" In case of linear fitting, compute the LTS regression
 318                      coefficients (robust(resistant) regression), using
 319                      the default fraction of good points
 320                "ROB=0.x" - compute the LTS regression coefficients, using
 321                            0.x as a fraction of good points
 322
 323   */
 324   printf("AliEMCALHistoUtilities::drawGraph started \n");
 325   printf("Drawing opt |%s| : Fitting opt |%s|\n", opt, optFit);
 326
 327     TGraph *gr = new TGraph(n, x, y);
 328     gr->SetMarkerColor(markerColor);
 329     gr->SetLineColor(markerColor);
 330     gr->SetMarkerStyle(markerStyle);
 331     gr->SetTitle(tit);
 332     gr->Draw(opt);
 333
 334     TString ctmp(opt);
 335     if(ctmp.Contains("A")) {
 336        gr->GetHistogram()->SetXTitle(xTit);
 337        gr->GetHistogram()->SetYTitle(yTit);
 338     }
 339     ctmp = optFit;
 340     if(ifun>=0) {
 341       TString sf("pol"); sf += ifun;
 342       gr->Fit(sf.Data(),optFit);
 343       printf("\n ** Fit by Polynomial of degree %i : %s **\n", ifun, sf.Data());
 344     } else if(!ctmp.Contains("no",TString::kIgnoreCase)){
 345       printf("\n ** ifun %i : %s **\n", ifun, fun);
 346       gr->Fit(fun, optFit);
 347     }
 348
 349     return gr;
 350 }
 351
 352 TGraphErrors *AliEMCALHistoUtilities::DrawGraphErrors(const Int_t n,Double_t *x,Double_t *y,Double_t *ex,
 353 Double_t *ey, Int_t markerColor,  Int_t markerStyle, const char* opt, const char* tit,
 354 const char* xTit,const char* yTit, Int_t ifun, const char *optFit, const char *fun)
 355 {
 356   // Oct 15, 2007
 357   printf("AliEMCALHistoUtilities::drawGraphErrors started \n");
 358   printf("Drawing opt |%s| : ifun %i: Fitting opt |%s|, fun |%s|\n",
 359          opt, ifun, optFit, fun);
 360
 361   TGraphErrors *gr = new TGraphErrors(n, x,y,ex,ey);
 362   gr->SetMarkerColor(markerColor);
 363   gr->SetLineColor(markerColor);
 364   gr->SetMarkerStyle(markerStyle);
 365   if(tit&&strlen(tit)>0) gr->SetTitle(tit);
 366
 367   TString ctmp(opt);
 368   if(ctmp.Contains("A")) {
 369      gr->GetHistogram()->SetXTitle(xTit);
 370      gr->GetHistogram()->SetYTitle(yTit);
 371   }
 372   if(ifun>0) {
 373     if(ifun != 999) {
 374       TString sf("pol"); sf += ifun;
 375       gr->Fit(sf.Data(),optFit);
 376       printf("\n ** Fit by Polynomial of degree %i : %s **\n", ifun, sf.Data());
 377     } else {
 378       gr->Fit(fun, optFit);
 379       printf("\n ** Fit by %s **\n", fun);
 380     }
 381   } else {
 382     if(strlen(optFit)) {
 383       printf("\n ** ifun %i : %s **\n", ifun, fun);
 384       gr->Fit(fun, optFit);
 385     }
 386   }
 387
 388   gr->Draw(opt);
 389
 390   return gr;
 391 }
 392
 393 TF1* AliEMCALHistoUtilities::GetResolutionFunction(const char *opt, TString &latexName)
 394 {
 395   // Oct 15, 2007
 396   TString sopt(opt);
 397   sopt.ToUpper();
 398   TF1 *fres=0;
 399   if      (sopt.Contains("FRES1")) {
 400     fres = new TF1("fres","[0]+[1]/sqrt(x)", 0.0, 101.);
 401     latexName = "#frac{#sigma_{E}}{E} = A+#frac{B}{#sqrt{E}}";
 402   } else if(sopt.Contains("FRES2")) {
 403     fres = new TF1("fres","sqrt([0]*[0]+[1]*[1]/x)", 0.0, 101.);
 404     latexName = "#sqrt{A^{2}+#frac{B^{2}}{E}}";
 405   }
 406   if(fres) {
 407     fres->SetParName(0,"A");
 408     fres->SetParName(1,"B");
 409
 410     fres->SetParameter(0, 2.0);
 411     fres->SetParameter(1, 6.6);
 412     fres->SetLineWidth(2);
 413     fres->SetLineColor(kRed);
 414   }
 415   return fres;
 416 }
 417
 418 void AliEMCALHistoUtilities::InitChain(TChain *chain, const char* nameListOfFiles, Int_t nFileMax)
 419 {
 420   // Read name of files from text file with nameListOfFiles and added to the chain
 421   if(chain==0 || nameListOfFiles==0) return;
 422
 423   ifstream fin;
 424   fin.open(nameListOfFiles);
 425   if (!fin.is_open()) {
 426     cout << "Input file "<<nameListOfFiles<<" cannot be opened.\n";
 427     return;
 428   }
 429
 430   Int_t nfiles = 0; // number of files in chain
 431   char nf[200];     // name of current file
 432   while (!fin.getline(nf,200).eof()) {
 433     if(!fin.good()) break;
 434     chain->Add(nf);
 435     nfiles++;
 436     cout<<nfiles<<" "<<nf<<endl;
 437     if(nFileMax && nfiles>=nFileMax) break;
 438   }
 439   fin.close();
 440   //
 441   cout << " \n ********** <I> Accepted file "<< nfiles << "********* \n"<<endl;
 442   //  chainEsd->Print();
 443   //  chain->Lookup();
 444 }
 445
 446 //AliRunLoader* AliEMCALHistoUtilities::InitKinematics(const Int_t nev, const char* galiceName)
 447 //{
 448 //  // Oct 15, 2007
 449 //  // nev == 0 - new file
 450 //  static AliRunLoader *rl = 0;
 451 //
 452 //  if((rl == 0 || nev==0) && galiceName) {
 453 //    //printf("<I> AliEMCALHistoUtilities::InitKinematics() : nev %i : rl %p : %s (IN)\n",
 454 //      // nev, rl, galiceName);
 455 //    if(rl)  {
 456 //      rl->UnloadgAlice();
 457 //      delete rl;
 458 //    }
 459 //    rl = AliRunLoader::Open(galiceName,AliConfig::GetDefaultEventFolderName(),"read");
 460 //    rl->LoadgAlice(); // obligatory
 461 //    //printf("<I> AliEMCALHistoUtilities::InitKinematics() : nev %i : rl %p : %s (OUT)\n",
 462 //       //nev, rl, galiceName);
 463 //  }
 464 //  if(rl) {
 465 //    rl->GetEvent(nev);
 466 //    rl->LoadKinematics();
 467 //    /* Get what you need after that
 468 //      rl->LoadHits();
 469 //      AliStack *stack=rl->Stack();
 470 //      AliESDCaloCluster * clus = esd->GetCaloCluster(n);
 471 //      Int_t label = clus->GetLabel(); // what is this
 472 //      TParticle *primary=stack->Particle(label);
 473 //    */
 474 //  }
 475 //  return rl;
 476 //}
 477
 478 //AliRunLoader* AliEMCALHistoUtilities::GetRunLoader(const Int_t nev, const Char_t* galiceName,
 479 //                                     const Char_t* eventFolderName, AliRunLoader* rlOld)
 480 //{
 481 //  // Nov 26, 2007
 482 //  // nev == 0 - new file
 483 //  AliRunLoader *rl = 0;
 484 //
 485 //  if(nev==0 && galiceName) {
 486 //    printf("<I> AliEMCALHistoUtilities::GetLoader() : nev %i : %s (IN)\n",
 487 //         nev, galiceName);
 488 //    if(rlOld) {
 489 //      rlOld->UnloadgAlice();
 490 //      delete rlOld;
 491 //    }
 492 //    TString folderName(eventFolderName);
 493 //    if(folderName.Length() == 0) folderName = AliConfig::GetDefaultEventFolderName();
 494 //    rl = AliRunLoader::Open(galiceName, folderName.Data(),"read");
 495 //    rl->LoadgAlice(); // obligatory
 496 //    printf("<I> AliEMCALHistoUtilities::GetLoader() : nev %i : %s : %s (OUT)\n",
 497 //         nev, galiceName, folderName.Data());
 498 //  } else {
 499 //    rl = rlOld;
 500 //  }
 501 //  if(rl) rl->LoadgAlice(); // obligatory
 502 //  // if(rl) rl->GetEvent(nev);
 503 //  return rl;
 504 //}
 505
 506 Double_t AliEMCALHistoUtilities::GetMomentum(const char* nameListOfFiles)
 507 {
 508   // Get momentum value from string  - like /....100GEV/..
 509   Double_t p=-1; // negative if undefined
 510   TString st(nameListOfFiles);
 511   if(st.Length()>=5) {
 512     st.ToUpper();
 513     Ssiz_t ind1 = st.Index("GEV"), ind2=0;
 514     if(ind1>0) {
 515       for(Int_t i=2; i<=10; i++) {
 516         ind2 = st.Index("/",ind1-i);
 517         if(ind2>0 && ind2<ind1) break;
 518       }
 519       TString mom  = st(ind2+1, ind1-ind2-1);
 520       if(mom.IsFloat()) p = mom.Atof();
 521       printf(" dir |%s| : mom |%s| : p %f : ind2,1 %i->%i\n", st.Data(), mom.Data(), p, ind2, ind1);
 522     }
 523   }
 524   return p;
 525 }
 526
 527 int AliEMCALHistoUtilities::ParseString(const TString &topt, TObjArray &Opt)
 528 {
 529   // Moved from AliEMCALGeometry
 530   // Feb 06, 2006
 531   Ssiz_t begin, index, end, end2;
 532   begin = index = end = end2 = 0;
 533   TRegexp separator("[^ ;,\\t\\s/]+");
 534   while ( (begin < topt.Length()) && (index != kNPOS) ) {
 535     // loop over given options
 536     index = topt.Index(separator,&end,begin);
 537     if (index >= 0 && end >= 1) {
 538       TString substring(topt(index,end));
 539       Opt.Add(new TObjString(substring.Data()));
 540     }
 541     begin += end+1;
 542   }
 543   return Opt.GetEntries();
 544 }
 545
 546 // Analysis utilites
 547 //Bool_t AliEMCALHistoUtilities::GetLorentzVectorFromESDCluster(TLorentzVector &v, const AliESDCaloCluster* cl)
 548 //{
 549 //  // May 8, 2007
 550 //  static Double_t e=0.0;
 551 //  static Float_t pos[3];
 552 //  static TVector3 gpos;
 553 //  if(cl==0) return kFALSE;
 554 //
 555 //  e = Double_t(cl->E());
 556 //  if(e<=0.0) {
 557 //    printf(" negative cluster energy : %f \n", e);
 558 //    return kFALSE;
 559 //  }
 560 //  cl->GetPosition(pos);
 561 //  gpos.SetXYZ(Double_t(pos[0]), Double_t(pos[1]), Double_t(pos[2]));
 562 //  gpos.SetMag(e);
 563 //  v.SetVectM(gpos, 0.0);
 564 //
 565 //  return kTRUE;
 566 //}
 567
 568 //Bool_t AliEMCALHistoUtilities::GetLorentzVectorFromRecPoint(TLorentzVector &v, const AliEMCALRecPoint  *rp)
 569 //{
 570 //  // Jun 20, 2007
 571 //  static Double_t e=0.0;
 572 //  static TVector3 gpos;
 573 //  if(rp==0) return kFALSE;
 574 //
 575 //  e = Double_t(rp->GetPointEnergy());
 576 //  if(e<=0.0) {
 577 //    printf(" negative rec.point energy : %f \n", e);
 578 //    return kFALSE;
 579 //  }
 580 //  rp->GetGlobalPosition(gpos);
 581 //  gpos.SetMag(e);
 582 //  v.SetVectM(gpos, 0.0);
 583 //
 584 //  return kTRUE;
 585 //}
 586 //
 587 //// Drawing:
 588 //
 589 void AliEMCALHistoUtilities::DrawHist(TH1* hid,int lineWidth,int lineColor,const char* opt, int lineStyle)
 590 {
 591   // Oct 15, 2007
 592   TString sopt;
 593   sopt = opt;
 594   if(!hid) return;
 595   printf(" Hist. %s : option |%s| \n", hid->GetName(), opt);
 596   if(hid && hid->GetEntries()>=1.) {
 597     if(lineWidth) hid->SetLineWidth(lineWidth);
 598     if(lineColor) hid->SetLineColor(lineColor);
 599     if(lineStyle) hid->SetLineStyle(lineStyle);
 600     if(sopt.Contains("stat",TString::kIgnoreCase)) hid->SetStats(kTRUE);
 601     hid->Draw(opt);
 602   } else {
 603     if   (strcmp(opt,"empty")==0) hid->Draw();
 604     else printf(" has fewer entries %i or hid is zero\n", (int) hid->GetEntries());
 605   }
 606 }
 607
 608 //
 609 //// Fitting:
 610 //
 611 TF1* AliEMCALHistoUtilities::Gausi(const char *addName,double xmi,double xma,double N,double mean,double sig,double width)
 612 {
 613   // Fit by gaus where first parameter is the number of events under ga
 614   // Oct 15, 2007
 615   TString name("gi");
 616   name += addName;
 617   TF1 *f = new TF1(name.Data(), Gi, xmi, xma, 4);
 618   f->SetParNames("INTEGRAL","MEAN","SIGMA","WIDTH");
 619
 620   f->SetParameter(0,N);
 621   f->SetParameter(1,mean);
 622   f->SetParameter(2,sig);
 623
 624   f->FixParameter(3,width); // width of histogramm bin
 625   return f;
 626 }
 627
 628 TF1* AliEMCALHistoUtilities::Gausi(const char *addName, double xmi, double xma, TH1 *h)
 629 {
 630   // You can change parameters after that if you don't like this assumption
 631   if(h) return Gausi(addName, xmi, xma, h->Integral(), h->GetMean(),h->GetRMS(), h->GetBinWidth(1));
 632   else  return 0;
 633 }
 634
 635 TF1* AliEMCALHistoUtilities::GausiPol2(const char *addName,double xmi,double xma, TF1 *g, TF1* bg)
 636 {
 637   // Fit by gaus where first parameter is the number of events under ga
 638   TString name("giPol2");
 639   name += addName;
 640   TF1 *f = new TF1(name.Data(), GiPol2, xmi, xma, 7);
 641   f->SetParNames("INTEGRAL","MEAN","SIGMA","WIDTH","a0","a1","a2");
 642
 643   if(g) {
 644     for(int i=0; i<4; i++) f->SetParameter(i, g->GetParameter(i));
 645     f->FixParameter(3,g->GetParameter(3));
 646   }
 647
 648   if(bg) {
 649     for(int i=4; i<7; i++) f->SetParameter(i, bg->GetParameter(i+4));
 650   }
 651   f->SetLineColor(kRed);
 652   return f;
 653 }
 654
 655 Double_t AliEMCALHistoUtilities::Gi(Double_t *x, Double_t *par)
 656 {
 657   // First parameter is integral (number events under gaus)
 658   // Forth parameter (par[3]) is width of histogram bin
 659   // gaus(0) is a substitute for [0]*exp(-0.5*((x-[1])/[2])**2)
 660
 661   static Double_t c = TMath::Sqrt(2.*TMath::Pi()), y=0.0, f=0.0; // sqrt(2.*pi)
 662
 663   y  = (x[0]-par[1])/par[2];
 664   f  = par[0]*par[3]/(par[2]*c) * TMath::Exp(-0.5*y*y);
 665
 666   return f;
 667 }
 668
 669 Double_t AliEMCALHistoUtilities::GiPol2(Double_t *x, Double_t *par)
 670 {
 671   // First parameter is integral (number events under gaus)
 672   // Forth parameter (par[3]) is width of histogram bin
 673   // gaus(0) is a substitute for [0]*exp(-0.5*((x-[1])/[2])**2)
 674   // + pol2 -> 7 parameters
 675   static Double_t c = TMath::Sqrt(2.*TMath::Pi()), y=0.0, f=0.0; // sqrt(2.*pi)
 676
 677   y  = (x[0]-par[1])/par[2];
 678   f  = par[0]*par[3]/(par[2]*c) * TMath::Exp(-0.5*y*y);
 679
 680   f += par[4] + par[5]*x[0] + par[6]*x[0]*x[0];
 681
 682   return f;
 683 }
 684
 685 // Calibration stuff
 686 Double_t AliEMCALHistoUtilities::GetCorrectionCoefficientForGamma1(const Double_t eRec)
 687 {
 688   // Correction to rec.energy - Jul 15, 2007
 689   // E(gamma) = corr * E(eRec);
 690   /* corr = p0*(eRec-7.5)+p1*(eRec-7.5)*(eRec-7.5)+p2; 0.0<eRec<10.0
 691    1  p0           6.07157e-05   1.15179e-04  -0.00000e+00   1.20997e-03
 692    2  p1           1.50019e-04   3.13566e-05  -0.00000e+00   1.22531e-02
 693    3  p2           9.99019e-01   4.08251e-04  -0.00000e+00   1.40606e-03
 694
 695      corr = p3 + p4*eRec + p5*eRec*eRec
 696    1  p3           9.97135e-01   5.31970e-04   1.37962e-09   1.68120e-08
 697    2  p4           3.15740e-04   2.53371e-05   1.11475e-11   1.74192e-04
 698    3  p5          -1.35383e-06   2.19495e-07  -5.82864e-13   4.52277e-02
 699   */
 700   static Double_t p0=6.07157e-05, p1=1.50019e-04, p2=9.99019e-01;
 701   static Double_t p3=9.97135e-01, p4=3.15740e-04, p5=-1.35383e-06;
 702   static Double_t corr=1.0;
 703   if(eRec>=0.0 && eRec <=10.0) {
 704     corr = p0*(eRec-7.5) + p1*(eRec-7.5)*(eRec-7.5) + p2;
 705   } else if(eRec>10.0 && eRec <=105.0) {
 706     corr = p3 + p4*eRec + p5*eRec*eRec;
 707   }
 708   //printf(" eRec %f | corr %f \n", eRec, corr);
 709   return corr;
 710 }
 711
 712 Double_t AliEMCALHistoUtilities::GetCorrectedEnergyForGamma1(const Double_t eRec)
 713 {
 714   return GetCorrectionCoefficientForGamma1(eRec) * eRec;
 715 }
 716
 717 // Trigger
 718 TList* AliEMCALHistoUtilities::GetTriggersListOfHists(const Int_t scale, const Int_t nTrig, const Bool_t toBrowser)
 719 {
 720   // Oct 22, 2007 - trigger technical assurance
 721   gROOT->cd();
 722   TH1::AddDirectory(1);
 723
 724   new TH1F("00_hXpos2x2", "X coord. of max Amp 2x2",100, -500., +500.);
 725   new TH1F("01_hYpos2x2", "Y coord. of max Amp 2x2",100, -500., +500.);
 726   new TH1F("02_hZpos2x2", "Z coord. of max Amp 2x2",100, -500., +500.);
 727   new TH1F("03_hXposnxn", "X coord. of max Amp NXN",100, -500., +500.);
 728   new TH1F("04_hYposnxn", "Y coord. of max Amp NXN",100, -500., +500.);
 729   new TH1F("05_hZposnxn", "Z coord. of max Amp NXN",100, -500., +500.);
 730   // May 7, 2008 - jet trigger
 731   new TH1F("06_hJetTriggerPhi", "%phi of COG of jet trigger patch", 110, 80., 190.);
 732   new TH1F("07_hJetTriggerEta", "%eta of COG of jet trigger patch", 70, -0.7, +0.7);
 733   //
 734   new TH1F("08_hMaxAmp2x2", "max Amp 2x2", 1000, 0.0, pow(2.,14.));
 735   new TH1F("09_hAmpOutOf2x2", "Amp out of patch 2x2", 1000, 0.0, pow(2.,14.));
 736   new TH1F("10_hMaxAmpnxn", "max Amp NXN", 1000, 0.0, pow(2.,14.));
 737   new TH1F("11_hAmpOutOfnxn", "Amp out of patch nxn", 1000, 0.0, pow(2.,14.));
 738   // May 7, 2008 - jet trigger
 739   for(Int_t i=0; i<nTrig; i++) {
 740     new TH1F(Form("%2.2i_hJetPatchAmp%2.2i",i+12,i), Form("jet patch amplitude : jet trig %i",i),
 741     1000, 0.0, pow(2.,14.));
 742   }
 743   // For checking
 744   Double_t maxEdigit=100., maxN=1000., maxPC = 200.;
 745   if(scale==1) {
 746     maxN  *= 10.;
 747     maxPC *= 10.;
 748   }
 749   Int_t ind = 12+nTrig;
 750   new TH1F(Form("%2.2i_hDigitsAmp",ind++), "amplitude of digits (PC)  ", 1001, -0.5, 1000.5);
 751   new TH1F(Form("%2.2i_hDigitsE",ind++), " energy of digits (PC)", 1000, 0.0, maxEdigit);
 752   new TH1F(Form("%2.2i_hNDigitsInPCs",ind++), " number of digits in PC's", 1000, 0.0, maxN);
 753   new TH1F(Form("%2.2i_hEinInPCs",ind++), " energy in PC's", 200, 0.0, maxPC);
 754
 755   return MoveHistsToList("TriggerLiOfHists", toBrowser);
 756 }
 757
 758 void AliEMCALHistoUtilities::FillTriggersListOfHists(TList *l, TArrayF *triggerPosition, TArrayF *triggerAmplitudes)
 759 {
 760   // Oct 22, 2007 - trigger technical assurance
 761   if(l==0) {
 762     printf("<E> FillTriggersListOfHists() : list of hists undefined. \n");
 763     return;
 764   }
 765   for(int i=0; i<triggerPosition->GetSize(); i++) {
 766     FillH1(l, i, double(triggerPosition->At(i)));
 767   }
 768
 769   for(int i=0; i<triggerAmplitudes->GetSize(); i++) {
 770     FillH1(l, triggerPosition->GetSize() + i, double(triggerAmplitudes->At(i)) );
 771   }
 772 }
 773
 774 // Jet(s) kinematics
 775 TList* AliEMCALHistoUtilities::GetJetsListOfHists(Int_t njet, Bool_t toBrowser)
 776 {
 777   // Oct 30, 2007
 778   gROOT->cd();
 779   TH1::AddDirectory(1);
 780   new TH1F("00_nJet", "number of jets in Pythia",5, -0.5, +4.5);
 781   int ic=1;
 782   for(Int_t ij=0; ij<njet; ij++)
 783   {
 784     new TH1F(Form("%2.2i_EtaOfJet%i",ic++,ij), Form("#eta of jet (%i)",ij),80, -2., 2.);
 785     new TH1F(Form("%2.2i_ThetaOfJet%i",ic++,ij), Form("#theta(degree) of jet (%i)",ij),
 786     180, 0.0, 180.);
 787     new TH1F(Form("%2.2i_PhiOfJet%i",ic++,ij), Form("#phi(degree) of jet (%i)",ij),
 788     180, 0.0, 360.);
 789     new TH1F(Form("%2.2i_PtOfJet%i",ic++,ij), Form(" p_{T} of jet (%i)",ij),
 790     200, 0.0, 200.);
 791     new TH1F(Form("%2.2i_EOfJet%i",ic++,ij), Form(" E of jet (%i)",ij),
 792     200, 0.0, 200.);
 793   }
 794
 795   return MoveHistsToList("JetLiOfHists", toBrowser);
 796 }
 797
 798 //void  AliEMCALHistoUtilities::FillJetKineListOfHists(TList *l, AliRunLoader* rl, TLorentzVector &goodJet)
 799 //{
 800 //  // Oct 30, 2007; Nov 07;
 801 //  // goodJet - output; only one jet with EMCAL acceptance
 802 //
 803 //  // Bad case
 804 //  goodJet.SetPxPyPzE(0., 0., 0., 0.);
 805 //
 806 //  if(l==0 || rl==0) return;
 807 //  //try to get trigger jet info
 808 //  AliHeader* alih = rl->GetHeader();
 809 //  if (alih == 0) return;
 810 //
 811 //  AliGenEventHeader * genh = alih->GenEventHeader();
 812 //  if (genh == 0) return;
 813 //
 814 //  AliGenPythiaEventHeader* genhpy = dynamic_cast<AliGenPythiaEventHeader *>(genh);
 815 //  if(genhpy==0) return; // Not Pythia
 816 //
 817 //  Int_t nj = genhpy->NTriggerJets();
 818 //  FillH1(l, 0, double(nj));
 819 //
 820 //  TLorentzVector* jets[4];
 821 //  nj = nj>4?4:nj;
 822 //
 823 //  int ic=1;
 824 //  for (Int_t i=0; i< nj; i++) {
 825 //     Float_t p[4];
 826 //     genhpy->TriggerJet(i,p);
 827 //     jets[i] = new TLorentzVector(p);
 828 //     FillH1(l, ic++, jets[i]->Eta());
 829 //     FillH1(l, ic++, jets[i]->Theta()*TMath::RadToDeg());
 830 //     FillH1(l, ic++, TVector2::Phi_0_2pi(jets[i]->Phi())*TMath::RadToDeg());
 831 //     FillH1(l, ic++, jets[i]->Pt());
 832 //     FillH1(l, ic++, jets[i]->E());
 833 //
 834 //     //printf(" %i pj %f %f %f %f : ic %i\n", i, p[0], p[1], p[2], p[3], ic);
 835 //     if(ic >= l->GetSize()) break;
 836 //  }
 837 //  if(nj==1) {
 838 //    Double_t eta=jets[0]->Eta(), phi=TVector2::Phi_0_2pi(jets[0]->Phi())*TMath::RadToDeg();
 839 //    if(TMath::Abs(eta)<0.5 && (phi>90.&&phi<180.)) {
 840 //      goodJet = (*jets[0]);
 841 //    }
 842 //  }
 843 //}