EMCAL/AliEMCALRecPointsQaESDSelector.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 /* $Log$ */
  17
  18 //_________________________________________________________________________
  19 // This is selector for making a set of histogramms from ESD for rec.points
  20 //*--  Authors: Aleksei Pavlinov (WSU)
  21 //*--  Feb 17, 2007 - May 23, 2007
  22 //*--  Pi0 calibration
  23 //_________________________________________________________________________
  24
  25 #include "AliEMCALRecPointsQaESDSelector.h"
  26 #include "AliEMCALFolder.h"
  27 #include "AliEMCALSuperModule.h"
  28 #include "AliEMCALCell.h"
  29 #include "AliEMCALPi0SelectionParam.h"
  30 #include "AliEMCALHistoUtilities.h"
  31 #include "AliEMCALCalibCoefs.h"
  32 #include "AliEMCALGeometry.h"
  33 #include "AliLog.h"
  34 #include "AliESD.h"
  35 #include "AliEMCALRecPoint.h"
  36
  37 #include <assert.h>
  38 #include <map>
  39 #include <string>
  40
  41 #include <TROOT.h>
  42 #include <TStyle.h>
  43 #include <TSystem.h>
  44 #include <TCanvas.h>
  45 #include <TVector3.h>
  46 #include <TLorentzVector.h>
  47 #include <TChain.h>
  48 #include <TFile.h>
  49 #include <TH1F.h>
  50 #include <TF1.h>
  51 #include <TMath.h>
  52 #include <TArrayS.h>
  53 #include <TBrowser.h>
  54 #include <TDataSet.h>
  55
  56 using namespace std;
  57
  58 typedef  AliEMCALHistoUtilities u;
  59
  60 ClassImp(AliEMCALRecPointsQaESDSelector)
  61
  62 Int_t  nmaxCell = 4*12*24*11;
  63
  64 AliEMCALRecPointsQaESDSelector::AliEMCALRecPointsQaESDSelector() :
  65   AliSelector(),
  66   fChain(0),
  67   fLofHistsPC(0),
  68   fLofHistsRP(0),
  69   fEmcalPool(0),
  70   fEMCAL(0),
  71   fEMCALOld(0)
  72 {
  73   //
  74   // Constructor. Initialization of pointers
  75   //
  76 }
  77
  78 AliEMCALRecPointsQaESDSelector::~AliEMCALRecPointsQaESDSelector()
  79 {
  80   //
  81   // Destructor
  82   //
  83
  84   // histograms are in the output list and deleted when the output
  85   // list is deleted by the TSelector dtor
  86 }
  87
  88 void AliEMCALRecPointsQaESDSelector::Begin(TTree* tree)
  89 {
  90   // Begin function
  91
  92   AliSelector::Begin(tree);
  93 }
  94
  95 void AliEMCALRecPointsQaESDSelector::Init(TTree* tree)
  96 {
  97   // Init function
  98
  99   AliSelector::Init(tree);
 100 }
 101
 102 void AliEMCALRecPointsQaESDSelector::SlaveBegin(TTree* tree)
 103 {
 104   // The SlaveBegin() function is called after the Begin() function.
 105   // When running with PROOF SlaveBegin() is called on each slave server.
 106   // The tree argument is deprecated (on PROOF 0 is passed).
 107
 108   printf("**** <I> Slave Begin \n **** ");
 109
 110   AliSelector::SlaveBegin(tree);
 111 }
 112
 113 void AliEMCALRecPointsQaESDSelector::InitStructure(Int_t it)
 114 {
 115   AliEMCALGeometry::GetInstance("SHISH_TRD1_CURRENT_2X2"); // initialize geometry just once
 116   //AliEMCALGeometry::GetInstance(""); // initialize geometry just once
 117
 118   fLofHistsPC = DefineHistsOfRP("PseudoCl");
 119   fLofHistsRP = DefineHistsOfRP("RP");
 120
 121   fEmcalPool = new TObjectSet("PoolOfEMCAL");
 122   if(it == 1) {
 123     fEMCAL     = new AliEMCALFolder(it); // folder for first itteration
 124     fEmcalPool->Add(fEMCAL);
 125   }
 126 }
 127
 128 Bool_t AliEMCALRecPointsQaESDSelector::Process(Long64_t entry)
 129 {
 130   // The Process() function is called for each entry in the tree (or possibly
 131   // keyed object in the case of PROOF) to be processed. The entry argument
 132   // specifies which entry in the currently loaded tree is to be processed.
 133   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 134   // to read either all or the required parts of the data. When processing
 135   // keyed objects with PROOF, the object is already loaded and is available
 136   // via the fObject pointer.
 137   //
 138   // This function should contain the "body" of the analysis. It can contain
 139   // simple or elaborate selection criteria, run algorithms on the data
 140   // of the event and typically fill histograms.
 141
 142   // WARNING when a selector is used with a TChain, you must use
 143   //  the pointer to the current TTree to call GetEntry(entry).
 144   //  The entry is always the local entry number in the current tree.
 145   //  Assuming that fTree is the pointer to the TChain being processed,
 146   //  use fTree->GetTree()->GetEntry(entry).
 147
 148   if (AliSelector::Process(entry) == kFALSE)
 149     return kFALSE;
 150
 151   // Check prerequisites
 152   if (!fESD)
 153   {
 154     AliDebug(AliLog::kError, "ESD branch not available");
 155     return kFALSE;
 156   }
 157   static pi0SelectionParam* rPar = GetEmcalFolder()->GetPi0SelectionParRow(0);
 158
 159   static Int_t nEmcalClusters, indOfFirstEmcalRP, nEmcalRP,nEmcalPseudoClusters;
 160   nEmcalClusters    = fESD->GetNumberOfEMCALClusters();
 161   indOfFirstEmcalRP = fESD->GetFirstEMCALCluster();
 162   u::FillH1(fLofHistsRP, 1, double(indOfFirstEmcalRP));
 163
 164   static AliESDCaloCluster *cl = 0;
 165   nEmcalRP = nEmcalPseudoClusters = 0;
 166   TList *l=0;
 167   double eDigi=0;
 168
 169   TClonesArray lvM1("TLorentzVector", 100);  // for convenience; 100 is enough now
 170   TArrayI      indLv(100);                   // index of RP
 171
 172   static TLorentzVector v, vcl;
 173   int nrp = 0; // # of RP for gg analysis
 174   for(int i=indOfFirstEmcalRP; i<indOfFirstEmcalRP+nEmcalClusters; i++) {
 175     cl = fESD->GetCaloCluster(i);
 176     if(cl->GetClusterType() == AliESDCaloCluster::kPseudoCluster) {
 177       nEmcalPseudoClusters++;
 178       l = fLofHistsPC;
 179     } else if(cl->GetClusterType() == AliESDCaloCluster::kClusterv1){
 180       nEmcalRP++;
 181       if(fEMCAL->GetIterationNumber()>1) {
 182         AliEMCALRecPoint *rp = AliEMCALFolder::GetRecPoint(cl, fEMCAL->GetCCFirst(), fEMCAL->GetCCIn(), fLofHistsRP);
 183         //        if(rp->GetPointEnergy()>=rPar->eOfRpMin && u::GetLorentzVectorFromRecPoint(v, rp)) {
 184         if(u::GetLorentzVectorFromRecPoint(v, rp)) { // comparing with RP
 185           new(lvM1[nrp]) TLorentzVector(v);
 186           indLv[nrp] = i;
 187           nrp++;
 188           // Conroling of recalibration
 189           u::GetLorentzVectorFromESDCluster(vcl, cl);
 190           u::FillH1(fLofHistsRP, 11, vcl.P()-v.P());
 191           u::FillH1(fLofHistsRP, 12, TMath::RadToDeg()*(vcl.Theta()-v.Theta()));
 192           u::FillH1(fLofHistsRP, 13, TMath::RadToDeg()*(vcl.Phi()-v.Phi()));
 193           l = 0; // no filling
 194         }
 195         delete rp;
 196       } else { // first iteration
 197         //        if(cl->E()>=rPar->eOfRpMin && u::GetLorentzVectorFromESDCluster(v, cl)) {
 198         if(u::GetLorentzVectorFromESDCluster(v, cl)) { // comparing with RP
 199         // cut 0.4 GeV may be high !
 200           new(lvM1[nrp]) TLorentzVector(v);
 201           indLv[nrp] = i;
 202           nrp++;
 203           l = fLofHistsRP;
 204         }
 205       }
 206
 207     } else {
 208       printf(" wrong cluster type : %i\n", cl->GetClusterType());
 209       assert(0);
 210     }
 211     u::FillH1(l, 2, double(cl->GetClusterType()));
 212
 213     u::FillH1(l, 3, double(cl->GetNumberOfDigits()));
 214     u::FillH1(l, 4, double(cl->E()));
 215     // Cycle on digits (towers)
 216     Short_t *digiAmpl  = cl->GetDigitAmplitude()->GetArray();
 217     Short_t *digiTime  = cl->GetDigitTime()->GetArray();
 218     Short_t *digiAbsId = cl->GetDigitIndex()->GetArray();
 219     for(int id=0; id<cl->GetNumberOfDigits(); id++) {
 220       eDigi = double(digiAmpl[id]) / 500.; // See AliEMCALClusterizerv1
 221       //      if(eDigi <= 0.0) { // sometimes it is happen
 222       //if(eDigi > 10.0 && cl->GetClusterType() == AliESDCaloCluster::kClusterv1) {
 223       // printf(" %i digiAmpl %i : %f \n", id, int(digiAmpl[id]), eDigi);
 224       //}
 225       u::FillH1(l, 5, eDigi);
 226       u::FillH1(l, 6, double(digiTime[id]));
 227       u::FillH1(l, 7, double(digiAbsId[id]));
 228       if(int(digiAbsId[id]) >= nmaxCell) {
 229         printf(" id %i :  digiAbsId[id] %i (%i) : %s \n",
 230                id, int(digiAbsId[id]), nmaxCell, l->GetName());
 231       }
 232     }
 233   }
 234   u::FillH1(fLofHistsRP, 0, double(nEmcalRP));
 235   u::FillH1(fLofHistsPC, 0, double(nEmcalPseudoClusters));
 236
 237   static TLorentzVector *lv1=0, *lv2=0, lgg;
 238   static Double_t mgg, pgg;
 239   mgg = pgg = 0.;
 240   nrp = lvM1.GetEntriesFast();
 241   if(nrp >= 2) {
 242     // eff.mass analysis
 243     for(int i1=0; i1<nrp-1; i1++){
 244       lv1 = (TLorentzVector*)lvM1.At(i1);
 245       for(int i2=i1+1; i2<nrp; i2++){
 246         lv2 = (TLorentzVector*)lvM1.At(i2);
 247         lgg = (*lv1) + (*lv2);
 248         mgg = lgg.M();  // eff.mass
 249         pgg = lgg.P();  // momentum
 250         u::FillH1(fLofHistsRP, 8, mgg);
 251
 252         if((mgg>=rPar->massGGMin && mgg<=rPar->massGGMax)) {// pi0 candidates
 253           if((pgg>=rPar->momPi0Min && pgg>=rPar->momPi0Min)) {
 254             fEMCAL->FillPi0Candidate(mgg,fESD->GetCaloCluster(indLv[i1]),fESD->GetCaloCluster(indLv[i2]));
 255             u::FillH1(fLofHistsRP, 9, pgg);
 256             u::FillH1(fLofHistsRP,10, lv1->P());
 257             u::FillH1(fLofHistsRP,10, lv2->P());
 258           }
 259         }
 260       }
 261     }
 262   }
 263
 264   lvM1.Delete();
 265
 266   if(nEmcalClusters != (nEmcalRP+nEmcalPseudoClusters))
 267     Info("Process","nEmcalClusters %i : RP %i + PC %i ",nEmcalClusters, nEmcalRP, nEmcalPseudoClusters);
 268
 269   return kTRUE;
 270 }
 271
 272 void AliEMCALRecPointsQaESDSelector::SlaveTerminate()
 273 {
 274   // The SlaveTerminate() function is called after all entries or objects
 275   // have been processed. When running with PROOF SlaveTerminate() is called
 276   // on each slave server.
 277
 278   AliSelector::SlaveTerminate();
 279
 280   // Add the histograms to the output on each slave server
 281   if (!fOutput)
 282   {
 283     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
 284     return;
 285   }
 286
 287   // fOutput->Add(fMultiplicity);
 288 }
 289
 290 void AliEMCALRecPointsQaESDSelector::Terminate()
 291 {
 292   // The Terminate() function is the last function to be called during
 293   // a query. It always runs on the client, it can be used to present
 294   // the results graphically or save the results to file.
 295
 296   AliSelector::Terminate();
 297   /*
 298   fMultiplicity = dynamic_cast<TH1F*> (fOutput->FindObject("fMultiplicity"));
 299
 300   if (!fMultiplicity)
 301   {
 302     AliDebug(AliLog::kError, Form("ERROR: Histogram not available %p", (void*) fMultiplicity));
 303     return;
 304   }
 305
 306   new TCanvas;
 307   fMultiplicity->DrawCopy();
 308   */
 309   PrintInfo();
 310 }
 311 //
 312 TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfRP(const char *name)
 313 {
 314   double ADCchannelEC = 0.00305; // ~3mev per adc count
 315
 316   gROOT->cd();
 317   TH1::AddDirectory(1);
 318   new TH1F("00_EmcalMultiplicity", "multiplicity of EMCAL  ", 201, -0.5, 200.5); // real and pseudo
 319   new TH1F("01_IndexOfFirstEmcal", "index of first emcal rec.points ", 201, -0.5, 200.5);
 320
 321   new TH1F("02_NumberOf", "number of  ", 6, -0.5, 5.5);
 322   new TH1F("03_NumberOfDigitsIn", "number of digits(towers) in rec.points ", 101, -0.5, 100.5);
 323   new TH1F("04_EnergyOf", "energy of ", 1000, 0.0, 100.);
 324
 325   int nmax=10000;
 326   double xmi = ADCchannelEC/2., xma = xmi + ADCchannelEC*nmax;
 327   // All energy(momentum) unit is GeV if don't notice
 328   new TH1F("05_DigitEnergyIn", "digit energy in ", nmax, xmi, xma);
 329   new TH1F("06_DigitTimeIn", "digit time in 10ps(0.01ns) ", 1000, 0.0, 3.e+3); // ns/100 = 10 ps
 330   new TH1F("07_DigitAbsIdIn", "digit abs id in ", nmaxCell, -0.5, double(nmaxCell)-0.5);
 331   new TH1F("08_EffMass", "effective mass of #gamma,#gamma(m_{#pi^{0}}=134.9766 MeV)", 100, 0.0, 0.5);
 332   new TH1F("09_MomOfPi0Candidate", "momentum of #pi^{0} candidates (0.085 <mgg<0.185)", 600, 0.0, 30.0);
 333   new TH1F("10_MomOfRpPi0Candidate", "momentum of RP for #pi^{0} candidates (0.085 <mgg<0.185)", 600, 0.0, 30.0);
 334   // Recalibration staf
 335   new TH1F("11_MomClESD-RpRecalib", "difference of momentum cl(ESD) - rp(Recalib)", 200, -1., +1.);
 336   new TH1F("12_ThetaClESD-RpRecalib", "difference of #theta cl(ESD) - rp(Recalib) in degree", 100, -0.05, +0.05);
 337   new TH1F("13_PhiClESD-RpRecalib", "difference of #phi cl(ESD) - rp(Recalib) in degree ", 100, -0.05, +0.05);
 338   // Digi
 339   new TH1F("14_EDigiRecalib", "energy of digits after recalibration", 2000, 0.0, 20.);
 340   AliEMCALGeometry* g = AliEMCALGeometry::GetInstance();
 341   new TH1F("15_AbsIdRecalib", "abs Id of digits after recalibration", g->GetNCells(),-0.5,Double_t(g->GetNCells())-0.5);
 342
 343   TList *l = u::MoveHistsToList(Form("ListOfHists%s",name), kFALSE);
 344   u::AddToNameAndTitleToList(l, name, name);
 345
 346   return l;
 347 }
 348
 349 /* unused now
 350 TList *AliEMCALRecPointsQaESDSelector::DefineHistsOfTowers(const char *name)
 351 {
 352   //
 353   // ESD: towers information was saved to pseudo clusters
 354   //
 355   gROOT->cd();
 356   TH1::AddDirectory(1);
 357   new TH1F("00_EmcalMultiplicity", "multiplicity of EMCAL  ", 201, -0.5, 200.5); // number of pseudo RP
 358
 359   new TH1F("01_EnergyOf", "energy of ", 1000, 0.0, 100.);
 360
 361   TList *l = u::MoveHistsToList(Form("ListOfHists%s",name, " - ESD"), kFALSE);
 362   u::AddToNameAndTitleToList(l, name, name);
 363
 364   return l;
 365 }
 366 */
 367
 368 void AliEMCALRecPointsQaESDSelector::FitEffMassHist()
 369 {
 370   TH1* h = (TH1*)fLofHistsRP->At(8);
 371   AliEMCALCell::FitHist(h, GetName(), "draw");
 372 }
 373
 374 void AliEMCALRecPointsQaESDSelector::PrintInfo()
 375 {
 376   // Service routine
 377   TList *l[2] = {fLofHistsPC, fLofHistsRP};
 378   printf("\n");
 379   for(int i=0; i<2; i++){
 380     TH1F *h = (TH1F*)l[i]->At(2);
 381     printf(" %s \t: %i \n", h->GetTitle(), int(h->GetEntries()));
 382   }
 383 }
 384
 385 AliEMCALFolder*  AliEMCALRecPointsQaESDSelector::CreateEmcalFolder(const Int_t it)
 386 {
 387   AliEMCALFolder* newFolder = new AliEMCALFolder(it); // folder for iteration #it
 388   if(it>1) {
 389     fEMCALOld = fEMCAL;
 390     AliEMCALCalibCoefs* tabOldOut = fEMCALOld->GetCCOut();
 391     AliEMCALCalibCoefs* tabNewIn = new AliEMCALCalibCoefs(*tabOldOut);
 392     tabNewIn->SetName(AliEMCALFolder::fgkCCinName.Data());
 393     newFolder->Add(tabNewIn);
 394   }
 395   fEmcalPool->Add(newFolder);
 396   fEMCAL = newFolder;
 397
 398   return fEMCAL;
 399 }
 400
 401 AliEMCALFolder* AliEMCALRecPointsQaESDSelector::GetEmcalOldFolder(const Int_t nsm)
 402 {
 403   AliEMCALFolder* folder=0;
 404   if(fEmcalPool) folder =  (AliEMCALFolder*)fEmcalPool->FindByName(Form("EMCAL_%2.2i",nsm));
 405   return folder;
 406 }
 407
 408
 409 void AliEMCALRecPointsQaESDSelector::SetEmcalFolder(AliEMCALFolder* folder)
 410 {
 411   fEMCAL = folder;
 412   fEmcalPool->Add(fEMCAL);
 413 }
 414
 415 void AliEMCALRecPointsQaESDSelector::SetEmcalOldFolder(AliEMCALFolder* folder)
 416 {
 417   fEMCALOld = folder;
 418   fEmcalPool->Add(fEMCALOld);
 419 }
 420
 421 void AliEMCALRecPointsQaESDSelector::Browse(TBrowser* b)
 422 {
 423   if(fESD)        b->Add(fESD);
 424   if(fLofHistsPC) b->Add(fLofHistsPC);
 425   if(fLofHistsRP) b->Add(fLofHistsRP);
 426   if(fChain)      b->Add(fChain);
 427   if(fEmcalPool)  b->Add(fEmcalPool);
 428   //  if(u) b->Add(u);
 429 }
 430
 431 Bool_t AliEMCALRecPointsQaESDSelector::IsFolder() const
 432 {
 433   if(fESD || fLofHistsRP || fEmcalPool) return kTRUE;
 434   return kFALSE;
 435 }
 436
 437 void AliEMCALRecPointsQaESDSelector::ReadAllEmcalFolders()
 438 {
 439   if(fEmcalPool==0) {
 440     fEmcalPool = new TObjectSet("PoolOfEMCAL");
 441     for(Int_t it=1; it<=10; it++){
 442       AliEMCALFolder* fold = AliEMCALFolder::Read(Form("EMCALFOLDER_It%i_fit.root",it), "READ");
 443       if(fold) fEmcalPool->Add(fold);
 444     }
 445   }
 446 }
 447
 448 void AliEMCALRecPointsQaESDSelector::PictVsIterNumber(const Int_t ind, const Int_t nsm)
 449 {
 450   // Jun 27, 2007 - unfinished; which picture is the best
 451   if(ind<0 || ind>5) return;
 452   gROOT->cd();
 453   TH1::AddDirectory(1);
 454
 455   Int_t itMax = 10, it=0;
 456   map <int, char*> indName;
 457   indName[0] = "eff.mass";
 458   indName[3] = "mass of #pi_{0}";
 459   indName[4] = "resolution of #pi_{0}";
 460   indName[5] = "chi^{2}/ndf";
 461
 462   TH1F *hout = new TH1F(indName[ind], indName[ind], itMax, 0.5, double(itMax)+0.5);
 463
 464   TH1::AddDirectory(0);
 465   Double_t content, error;
 466   for(Int_t i=0; i<fEmcalPool->GetListSize(); i++) { // cycle on EMCAL folders
 467     AliEMCALFolder* folder = static_cast<AliEMCALFolder*>(fEmcalPool->At(i));
 468     if(folder==0) continue;
 469     it = folder->GetIterationNumber();
 470
 471     AliEMCALSuperModule* sm = folder->GetSuperModule(nsm);
 472     if(sm==0) continue;
 473
 474     TList* l = sm->GetHists();
 475     if(l==0) continue;
 476
 477     TH1F *hin = (TH1F*)l->At(ind);
 478     if(hin==0) continue;
 479
 480     if(ind !=0 ) {
 481       content = hin->GetMean();
 482       error   = hin->GetRMS();
 483     } else {
 484       sm->FitEffMassHist();
 485       TF1 *f = (TF1*)hin->GetListOfFunctions()->At(0);
 486       content = error = -1.;
 487       if(f) {
 488         //        content = f->GetParameter(1);
 489         //error   = f->GetParameter(2);
 490         content = f->GetParameter(2);
 491         error   = f->GetParError(2);
 492       }
 493     }
 494
 495     if(content > 0.0) {
 496       hout->SetBinContent(it, content);
 497       hout->SetBinError(it, error);
 498       printf(" it %i content %f +/- %f \n", it, content, error);
 499     }
 500   }
 501
 502   u::DrawHist(hout,2);
 503   hout->SetMinimum(0.0);
 504
 505 }