EMCAL/macros/TestEMCALRecPoint.C

   1 // Test Macro, shows how to load RecPoints, and how can we get
   2 // some of the parameters and variables.
   3 // Author: Gustavo Conesa
   4
   5 void TestEMCALRecPoint()
   6 {
   7
   8   // Getting EMCAL Detector and Geometry.
   9
  10   AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
  11
  12   if (rl == 0x0)
  13     cout<<"Can not instatiate the Run Loader"<<endl;
  14
  15   rl->LoadgAlice();//Needed to get geometry
  16
  17   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
  18     (rl->GetDetectorLoader("EMCAL"));
  19
  20   //TGeoManager::Import("geometry.root");
  21   AliGeomManager::LoadGeometry("geometry.root");
  22   AliRun * alirun   = rl->GetAliRun(); // Needed to get Geometry
  23   AliEMCALGeometry * geom ;
  24   if(alirun){
  25     AliEMCAL * emcal  = (AliEMCAL*)alirun->GetDetector("EMCAL");
  26     geom = emcal->GetGeometry();
  27   }
  28   else {
  29     cout<<"alirun not available, instantiate"<<endl;
  30     geom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETE12SMV1") ;
  31   }
  32
  33   //Load RecPoints
  34   rl->LoadRecPoints("EMCAL");
  35   //Get maximum number of events
  36   Int_t maxevent =  rl->GetNumberOfEvents();
  37   cout<<"Number of events "<<maxevent<<endl;
  38
  39   Int_t iEvent     = -1 ;
  40   Int_t iprim      = -1 ;
  41   Float_t energy   = -1 ;
  42   TVector3 gpos ;
  43
  44   AliEMCALRecPoint * rp;
  45   for ( iEvent=0; iEvent<maxevent; iEvent++)
  46   {
  47     cout <<  " ======> Event " << iEvent << endl ;
  48     //Load Event
  49     rl->GetEvent(iEvent);
  50     //      AliStack *sta=rl->Stack();
  51     //Fill array of digits
  52     TObjArray *rpoints ;//= emcalLoader->RecPoints();
  53
  54     TTree *treeR = emcalLoader->TreeR();
  55     TBranch * branchR = treeR->GetBranch("EMCALECARP");
  56     branchR->SetAddress(&rpoints);
  57     branchR->GetEntry(0);
  58
  59     if(!rpoints->GetEntries()) continue;
  60     cout <<  " ======> Event " << iEvent << endl ;
  61
  62     cout<<">> Entries "<<rpoints->GetEntries()<<endl;
  63
  64     //Get recpoints  from the list
  65     for(Int_t irp = 0; irp< rpoints->GetEntries();irp++){
  66       rp = static_cast<AliEMCALRecPoint *>(rpoints->At(irp)) ;
  67       //rp = emcalLoader->RecPoint(irp);
  68       if(rp != 0){
  69         {
  70           energy  = rp->GetEnergy(); //cluster energy
  71           cout<<"Energy "<<energy<<" PointEnergy "<<rp->GetPointEnergy()<<endl;
  72           rp->GetGlobalPosition(gpos);//global ALICE xyz position
  73           Int_t  primMult  = 0;
  74           //Int_t *primInts =  rp->GetPrimaries(primMult);
  75           //for (Int_t ipr=0; ipr<primMult; ipr++)
  76           //    cout<<"primlist "<<ipr<<" index "<< primInts[ipr]<<endl;
  77           //iprim=rp->GetPrimaryIndex() ;
  78           //    cout<<"Selected primary index "<<iprim<<endl;
  79           //  if(iprim>0){
  80           //TParticle *primary=sta->Particle(iprim);
  81           //cout<<"Primary phi "<<primary->Phi()*180/TMath::Pi()<<" Reconstructed phi "<<gpos.Phi()*180/TMath::Pi()<<"; Energy "<<primary->Energy()<<endl;
  82           //h->Fill((gpos.Phi()-primary->Phi())*180/TMath::Pi());
  83           cout<<"rec point "<<irp<<"; Energy "<<energy<<" Phi "<<gpos.Phi()*180/TMath::Pi()<<" Eta "<<gpos.Eta()<<" iprim "<<iprim<<endl;
  84           //}
  85           Float_t lmb[2];
  86           rp->GetElipsAxis(lmb);
  87           cout<<"lmb0 "<<lmb[0]<<" lmb1 "<<lmb[1]<<endl;
  88         }
  89       }
  90       else
  91         cout<<"recpoint pointer 0x0"<<endl;
  92     }
  93   }
  94 }
  95