EVE/EveDet/AliEveEMCALData.cxx

   1 //*********************************************************************
   2 // - AliEVE implementation -
   3 // Fill containers for visualisation of EMCAL data structures
   4 //    - read and store MC Hits
   5 //    - read and store digits from esds or runloader
   6 //    - read and store clusters from esds or runloader
   7 //
   8 // Author: Magali Estienne (magali.estienne@cern.ch)
   9 // June 30 2008
  10 //*********************************************************************
  11
  12 #include <TObject.h>
  13 #include <TEveUtil.h>
  14 #include <TTree.h>
  15 #include <TBranch.h>
  16 #include <TObjArray.h>
  17 #include <TRefArray.h>
  18 #include <TClonesArray.h>
  19 #include <TEvePointSet.h>
  20
  21 #include "AliRun.h"
  22 #include "AliRunLoader.h"
  23 #include "AliEMCALLoader.h"
  24 #include "AliESDEvent.h"
  25 #include "AliESDVertex.h"
  26 #include "AliEMCAL.h"
  27 #include "AliEMCALGeometry.h"
  28 #include "AliEMCALHit.h"
  29 #include "AliEMCALDigit.h"
  30 #include "AliEMCALRecPoint.h"
  31 #include "AliESDCaloCells.h"
  32 #include "AliESDCaloCluster.h"
  33
  34 #include "AliEveEMCALData.h"
  35 #include "AliEveEMCALSModule.h"
  36 #include "AliEveEMCALSModuleData.h"
  37
  38 ClassImp(AliEveEMCALData)
  39
  40 //______________________________________________________________________________
  41 AliEveEMCALData::AliEveEMCALData():
  42   TObject(),
  43   TEveRefCnt(),
  44   fEmcal(0x0),
  45   fGeom(0x0),
  46   fNode(0x0),
  47   fHMatrix(0),
  48   fTree(0x0),
  49   fESD(0x0),
  50   fNsm(12),
  51   fNsmfull(10),
  52   fNsmhalf(2),
  53   fSM(12),
  54   fSMfull(10),
  55   fSMhalf(2),
  56   fRunLoader(0),
  57   fDebug(0),
  58   fPoint(0)
  59 {
  60
  61   //
  62   // Constructor
  63   //
  64   //   for(Int_t i=0; i<12; i++)
  65   //     fSM[i] = 0x0;
  66   //   for(Int_t i=0; i<10; i++)
  67   //     fSMfull[i] = 0x0;
  68   //   fSMhalf[0] = fSMhalf[1] = 0x0;
  69   //  InitEMCALGeom();
  70   CreateAllSModules();
  71
  72
  73 }
  74
  75 //______________________________________________________________________________
  76 AliEveEMCALData::AliEveEMCALData(AliRunLoader* rl, TGeoNode* node, TGeoHMatrix* m):
  77   TObject(),
  78   TEveRefCnt(),
  79   fEmcal(0x0),
  80   fGeom(0x0),
  81   fNode(node),
  82   fHMatrix(m),
  83   fTree(0x0),
  84   fESD(0x0),
  85   fNsm(12),
  86   fNsmfull(10),
  87   fNsmhalf(2),
  88   fSM(12),
  89   fSMfull(10),
  90   fSMhalf(2),
  91   fRunLoader(rl),
  92   fDebug(0),
  93   fPoint(0)
  94 {
  95
  96   //
  97   // Constructor
  98   //
  99 //   for(Int_t i=0; i<12; i++)
 100 //     fSM[i] = 0x0;
 101 //   for(Int_t i=0; i<10; i++)
 102 //     fSMfull[i] = 0x0;
 103 //   fSMhalf[0] = fSMhalf[1] = 0x0;
 104   InitEMCALGeom(rl);
 105   CreateAllSModules();
 106
 107
 108 }
 109
 110 //______________________________________________________________________________
 111 AliEveEMCALData::~AliEveEMCALData()
 112 {
 113   //
 114   // Destructor
 115   //
 116
 117   DeleteSuperModules();
 118   delete fTree;
 119   delete fEmcal;
 120   delete fGeom;
 121   delete fNode;
 122   delete fHMatrix;
 123   delete fPoint;
 124 }
 125
 126 //______________________________________________________________________________
 127 void AliEveEMCALData::Reset()
 128 {
 129
 130 }
 131
 132 //______________________________________________________________________________
 133 AliEveEMCALData::AliEveEMCALData(const AliEveEMCALData &edata) :
 134   TObject(edata),
 135   TEveRefCnt(edata),
 136   fEmcal(edata.fEmcal),
 137   fGeom(edata.fGeom),
 138   fNode(edata.fNode),
 139   fHMatrix(edata.fHMatrix),
 140   fTree(edata.fTree),
 141   fESD(edata.fESD),
 142   fNsm(edata.fNsm),
 143   fNsmfull(edata.fNsmfull),
 144   fNsmhalf(edata.fNsmhalf),
 145   fSM(edata.fSM),
 146   fSMfull(edata.fSMfull),
 147   fSMhalf(edata.fSMhalf),
 148   fRunLoader(edata.fRunLoader),
 149   fDebug(edata.fDebug),
 150   fPoint(edata.fPoint)
 151 {
 152   //
 153   // Copy constructor
 154   //
 155   InitEMCALGeom(edata.fRunLoader);
 156   CreateAllSModules();
 157 }
 158
 159 //______________________________________________________________________________
 160 AliEveEMCALData& AliEveEMCALData::operator=(const AliEveEMCALData &edata)
 161 {
 162   //
 163   // Assignment operator
 164   //
 165
 166   if (this != &edata) {
 167
 168   }
 169
 170   return *this;
 171
 172 }
 173
 174 //______________________________________________________________________________
 175 void AliEveEMCALData::SetTree(TTree* tree)
 176 {
 177
 178   // Set digit-tree to be used for digit retrieval. Data is loaded on
 179   // demand.
 180
 181   fTree = tree;
 182
 183 }
 184
 185 //______________________________________________________________________________
 186 void AliEveEMCALData::SetESD(AliESDEvent* esd)
 187 {
 188   fESD = esd;
 189 }
 190
 191 //______________________________________________________________________________
 192 void AliEveEMCALData::SetNode(TGeoNode* node)
 193 {
 194   fNode = node;
 195 }
 196
 197 //______________________________________________________________________________
 198 void AliEveEMCALData::InitEMCALGeom(AliRunLoader* rl)
 199 {
 200 //   // Handle an individual point selection from GL.
 201
 202 //   fParent->SpawnEditor();
 203   fEmcal = (AliEMCAL*) rl->GetAliRun()->GetDetector("EMCAL");
 204   fGeom  = (AliEMCALGeometry*) fEmcal->GetGeometry();
 205
 206 }
 207
 208 //______________________________________________________________________________
 209 void AliEveEMCALData::GetGeomInfo(Int_t id, Int_t &iSupMod, Double_t& x, Double_t& y, Double_t& z)
 210 {
 211   // Get geometrical information from hit/digit/cluster absolute id
 212
 213   Int_t iTower  =  0 ;
 214   Int_t iIphi   =  0 ;
 215   Int_t iIeta   =  0 ;
 216   Int_t iphi    =  0 ;
 217   Int_t ieta    =  0 ;
 218
 219   //Geometry methods
 220   fGeom->GetCellIndex(id,iSupMod,iTower,iIphi,iIeta);
 221   //Gives SuperModule and Tower numbers
 222   fGeom->RelPosCellInSModule(id, x, y, z);
 223
 224 }
 225
 226 //______________________________________________________________________________
 227 void  AliEveEMCALData::CreateAllSModules()
 228 {
 229
 230   // Create all fNsm super modules
 231   for(Int_t sm = 0; sm < fNsm; sm++)
 232     CreateSModule(sm);
 233
 234 }
 235
 236 //______________________________________________________________________________
 237 void  AliEveEMCALData::CreateSModule(Int_t sm)
 238 {
 239   // Create super-module-data for sm if it does not exist already.
 240   if(fSM[sm] == 0) fSM[sm] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
 241   if(fSMfull[sm] == 0 && sm < 10) fSMfull[sm] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
 242   if(fSMhalf[sm-10] == 0 && sm > 10) fSMhalf[sm-10] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
 243 }
 244
 245 //______________________________________________________________________________
 246 void AliEveEMCALData::DropAllSModules()
 247 {
 248   // Drop data of all existing sectors.
 249
 250   for (Int_t sm = 0; sm < fNsm; sm++) {
 251     if (fSM[sm] != 0)
 252       fSM[sm]->DropData();
 253   }
 254 }
 255
 256 //______________________________________________________________________________
 257 void AliEveEMCALData::DeleteSuperModules()
 258 {
 259   //
 260   // delete all super module data
 261   //
 262
 263   for (Int_t sm = 0; sm < fNsm; sm++)
 264     {
 265       fSM[sm] = 0;
 266       delete fSM[sm];
 267     }
 268
 269   for(Int_t smf = 0; smf < fNsmfull; smf++)
 270     {
 271       fSMfull[smf] = 0;
 272       delete fSMfull[smf];
 273     }
 274
 275   for(Int_t smh = 0; smh < fNsmhalf; smh++)
 276     {
 277       fSMhalf[smh] = 0;
 278       delete fSMhalf[smh];
 279     }
 280
 281 }
 282
 283 //______________________________________________________________________________
 284 void AliEveEMCALData::LoadHits(TTree* t)
 285 {
 286
 287   // With TTreeH
 288
 289   // These are global coordinates !
 290   char form[1000];
 291   const char *selection = "";
 292   const char *varexp = "fX:fY:fZ";
 293   sprintf(form,"EMCAL Hits '%s'", selection);
 294   fPoint = new TEvePointSet(form);
 295
 296   TEvePointSelector ps(t, fPoint, varexp, selection);
 297   ps.Select();
 298
 299   if (fPoint->Size() == 0) {
 300     Warning("emcal_hits", Form("No hits match '%s'", selection));
 301     delete fPoint;
 302     //    return 0;
 303   }
 304
 305
 306 //   TObjArray *harr=NULL;
 307 //   TBranch *hbranch=t->GetBranch("EMCAL");
 308 //   hbranch->SetAddress(&harr);
 309
 310 //   if(hbranch->GetEvent(0)) {
 311 //     for(Int_t ih = 0; ih < harr->GetEntriesFast(); ih++) {
 312 //       AliEMCALHit* hit =(AliEMCALHit*)harr->UncheckedAt(ih);
 313 //       if(hit != 0){
 314 //      cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
 315 //      fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
 316 //       }
 317 //     }
 318 //   }
 319 //   //********************************
 320 //   // To be completed !!!
 321 //   //********************************
 322
 323 }
 324
 325 //______________________________________________________________________________
 326 void AliEveEMCALData::LoadHitsFromEMCALLoader(AliEMCALLoader* emcl)
 327 {
 328
 329   // Need to know the local position in each sm
 330
 331   // With EMCALLoader
 332   AliEMCALHit* hit;
 333
 334   //Fill array of hits
 335   TClonesArray *hits = (TClonesArray*)emcl->Hits();
 336
 337   //Get hits from the list
 338   for(Int_t ihit = 0; ihit< hits->GetEntries();ihit++){
 339     //cout<<">> idig "<<idig<<endl;
 340     hit = static_cast<AliEMCALHit *>(hits->At(ihit)) ;
 341
 342     if(hit != 0){
 343       cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
 344     }
 345   }
 346
 347   //********************************
 348   // To be completed !!!
 349   //********************************
 350
 351
 352
 353 }
 354
 355 //______________________________________________________________________________
 356 void AliEveEMCALData::LoadDigits(TTree *t)
 357 {
 358   //
 359   // load digits from the TreeD
 360   //
 361
 362   TClonesArray *digits = 0;
 363   t->SetBranchAddress("EMCAL", &digits);
 364   t->GetEntry(0);
 365
 366   Int_t nEnt = digits->GetEntriesFast();
 367   cout << "nEnt: " << nEnt << endl;
 368   AliEMCALDigit * dig;
 369
 370   Float_t amp   = -1 ;
 371   Int_t id      = -1 ;
 372   Int_t iSupMod =  0 ;
 373   Double_t x, y, z;
 374
 375   for (Int_t idig = 0; idig < nEnt; idig++)
 376     {
 377       dig = static_cast<AliEMCALDigit *>(digits->At(idig));
 378
 379       if(dig != 0) {
 380         id   = dig->GetId() ; //cell (digit) label
 381         amp  = dig->GetAmp(); //amplitude in cell (digit)
 382
 383         GetGeomInfo(id,iSupMod,x,y,z);
 384
 385         fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
 386 //      fSM[iSupMod]->SaveDigit(dig);
 387 //      if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
 388 //      if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
 389       }
 390       else {
 391         cout << "Digit object empty" << endl;
 392         return;
 393       }
 394     } // end loop digits
 395   cout << "after loop on digits !" << endl;
 396 }
 397
 398 //______________________________________________________________________________
 399 void AliEveEMCALData::LoadDigitsFromEMCALLoader(AliEMCALLoader* emcl)
 400 {
 401
 402   //
 403   // load digits from EMCAL Loader
 404   //
 405
 406   AliEMCALDigit* dig;
 407
 408   //Fill array of digits
 409   TClonesArray *digits = (TClonesArray*)emcl->Digits();
 410
 411   //Get digits from the list
 412
 413   Float_t amp   = -1 ;
 414   Int_t id      = -1 ;
 415   Int_t iSupMod =  0 ;
 416   Double_t x, y, z;
 417
 418    for(Int_t idig = 0; idig< digits->GetEntries();idig++){
 419
 420      dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
 421
 422      if(dig != 0){
 423        if(fDebug>1) cout << "Digit info " << dig->GetId() << " " << dig->GetAmp() << endl;
 424        id   = dig->GetId() ; //cell (digit) label
 425        amp  = dig->GetAmp(); //amplitude in cell (digit)
 426
 427        GetGeomInfo(id,iSupMod,x,y,z);
 428
 429        fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
 430      }
 431       else {
 432         cout << "Digit object empty" << endl;
 433         return;
 434       }
 435    } // end loop on digits
 436
 437 }
 438
 439 //______________________________________________________________________________
 440 void AliEveEMCALData::LoadDigitsFromESD()
 441 {
 442   //
 443   // Get digit information from esd
 444   //
 445
 446   AliESDCaloCells &cells= *(fESD->GetEMCALCells());
 447   Int_t ncell = cells.GetNumberOfCells() ;
 448   Int_t type = cells.GetType();
 449   Float_t pos[3] ;
 450   Int_t iSupMod =  0 ;
 451   Double_t x, y, z;
 452
 453   // Extract digit information from the ESDs
 454   for (Int_t icell=  0; icell <  ncell; icell++)
 455     {
 456       Int_t id   = cells.GetCellNumber(icell);
 457       Float_t amp  = cells.GetAmplitude(icell);
 458
 459       GetGeomInfo(id,iSupMod,x,y,z);
 460
 461       fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
 462       if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
 463       if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
 464
 465     } // end loop cells
 466 }
 467
 468 //______________________________________________________________________________
 469 void AliEveEMCALData::LoadRecPoints(TTree* t)
 470 {
 471   //*************************************************
 472   // To be improved !!!!!
 473   // Size and shape of cluster to be implemented
 474   //
 475   //*************************************************
 476
 477   // From TTreeR
 478   TObjArray *carr=NULL;
 479   TBranch *cbranch=t->GetBranch("EMCALECARP");
 480   cbranch->SetAddress(&carr);
 481
 482   if(cbranch->GetEvent(0)) {
 483     for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
 484       AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
 485       if(rp){
 486         cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
 487         Int_t iSupMod = rp->GetSuperModuleNumber();
 488         Float_t amp = rp->GetEnergy();
 489         TVector3 lpos;
 490         rp->GetLocalPosition(lpos);
 491
 492         fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
 493       }
 494     }
 495   }
 496
 497 }
 498
 499 //______________________________________________________________________________
 500 void AliEveEMCALData::LoadRecPointsFromEMCALLoader(AliEMCALLoader* emcl)
 501 {
 502   //*************************************************
 503   // To be improved !!!!!
 504   // Size and shape of cluster to be implemented
 505   //
 506   //*************************************************
 507
 508   // From EMCALLoader
 509   AliEMCALRecPoint* rp;
 510
 511   //Fill array of clusters
 512   TClonesArray *clusters = (TClonesArray*)emcl->RecPoints();
 513
 514   //Get clusters from the list
 515   for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
 516
 517     rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
 518
 519     if(rp){
 520        cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
 521        Int_t iSupMod = rp->GetSuperModuleNumber();
 522        Float_t amp = rp->GetEnergy();
 523        TVector3 lpos;
 524        rp->GetLocalPosition(lpos);
 525
 526        fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
 527     }
 528   }
 529
 530 }
 531
 532 //______________________________________________________________________________
 533 void AliEveEMCALData::LoadRecPointsFromESD()
 534 {
 535   Int_t iSupMod =  0 ;
 536   Double_t x, y, z;
 537   Int_t iSM =  0 ;
 538   Int_t iT  =  0 ;
 539   Int_t iIp   =  0 ;
 540   Int_t iIe   =  0 ;
 541   Int_t ip    =  0 ;
 542   Int_t ie    =  0 ;
 543   Double_t xd, yd, zd;
 544   Float_t pos[3] ;
 545
 546   // Get reconstructed vertex position
 547   AliESDVertex* primVertex =(AliESDVertex*) fESD->GetVertex();
 548   Double_t vertex_position[3] ;
 549   primVertex->GetXYZ(vertex_position) ;
 550
 551   //Get the CaloClusters
 552   //select EMCAL clusters only
 553   TRefArray * caloClusters  = new TRefArray();
 554   fESD->GetEMCALClusters(caloClusters);
 555   Int_t nclus = caloClusters->GetEntries();
 556   cout << "nclus: " << nclus << endl;
 557
 558   if(!caloClusters) return;
 559
 560   for (Int_t iclus =  0; iclus <  nclus; iclus++)
 561     {
 562       AliESDCaloCluster *clus = (AliESDCaloCluster *) caloClusters->At(iclus) ;
 563       //Get the cluster info
 564
 565       Float_t energy = clus->E() ;
 566       Int_t   eneInt = (Int_t) energy*500+0.5;
 567       Float_t disp   = clus->GetClusterDisp() ;
 568       Int_t   iprim  = clus->GetLabel();
 569
 570       clus->GetPosition(pos) ; // Global position
 571       TVector3 vpos(pos[0],pos[1],pos[2]) ;
 572       TLorentzVector p4 ;
 573       TVector3 p3;
 574       clus->GetMomentum(p4,vertex_position);
 575       p3.SetXYZ(p4[0],p4[1],p4[2]);
 576       Float_t eta = p3.Eta();
 577       Float_t phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
 578
 579       Int_t mult = clus->GetNCells() ;
 580       if(fDebug>1) {
 581         cout << "In cluster: " << iclus << ", ncells: " << mult << ", energy: " <<
 582           eneInt << ", disp: " << disp << endl;
 583         cout << "Cluster " << iclus << ", eta: " << eta << ", phi: " << phi << endl;
 584       }
 585
 586       Int_t clusId = 0;
 587       fGeom->GetAbsCellIdFromEtaPhi(eta,phi,clusId);
 588       if(fDebug>1) {
 589         cout << "Abs Cluster Id: " << clusId << ", xc: " << pos[0] <<
 590           ", yc: " << pos[1] << ", zc: " << pos[2] << endl;
 591       }
 592
 593       GetGeomInfo(clusId,iSupMod,x,y,z);
 594
 595       //******** Not used yet but will come  ********
 596       Int_t     digMult = clus->GetNCells() ;
 597       UShort_t *digID   = clus->GetCellsAbsId() ;
 598       for(Int_t i=0; i<digMult; i++){
 599         fGeom->RelPosCellInSModule(digID[i], xd, yd, zd);
 600         //Geometry methods
 601         fGeom->GetCellIndex(digID[i],iSM,iT,iIp,iIe);
 602         //Gives SuperModule and Tower numbers
 603
 604       } // end digit loop
 605       //*********************************************
 606
 607       fSM[iSupMod]->RegisterCluster(iSM,energy,x,y,z);
 608
 609     } // end cluster loop
 610 }
 611
 612 //______________________________________________________________________________
 613 AliEveEMCALSModuleData* AliEveEMCALData::GetSModuleData(Int_t sm)
 614 {
 615   //
 616   // return super module data
 617   //
 618
 619   if (sm < 0 || sm > fNsm)
 620     {
 621       printf("The number of super modules must be lower or equal to %d",fNsm);
 622       return 0;
 623     }
 624
 625   return fSM[sm];
 626 }
 627
 628 //______________________________________________________________________________
 629 void AliEveEMCALData::LoadRaw()
 630 {
 631
 632 }
 633