EMCAL/macros/TestEMCALData.C

   1 //Do a fast test of the different EMCAL data objects
   2 void TestEMCALData() {
   3
   4   TH1F* hE = new TH1F("hEmcalHits",    "Hits energy distribution in EMCAL",       100, 0., 10.) ;
   5   hE->Sumw2() ;
   6   TH1I* hM = new TH1I("hEmcalHitsMul", "Hits multiplicity distribution in EMCAL", 500, 0., 10000) ;
   7   hM->Sumw2() ;
   8
   9   TH1I * dA = new TH1I("hEmcalDigits",    "Digits amplitude distribution in EMCAL",    500, 0, 5000) ;
  10   dA->Sumw2() ;
  11   TH1I * dM = new TH1I("hEmcalDigitsMul", "Digits multiplicity distribution in EMCAL", 500, 0, 1000) ;
  12   dM->Sumw2() ;
  13
  14   TH1F * sE = new TH1F("hEmcalSDigits",    "SDigits energy distribution in EMCAL",    200, 0, 1000) ;
  15   sE->Sumw2() ;
  16   TH1I * sM = new TH1I("hEmcalSDigitsMul", "SDigits multiplicity distribution in EMCAL", 500, 0, 1000) ;
  17   sM->Sumw2() ;
  18
  19   TH1F * cE = new TH1F("hEmcalRecPoints",    "RecPoints energy distribution in EMCAL",    100, 0, 20) ;
  20   cE->Sumw2() ;
  21   TH1I * cM = new TH1I("hEmcalRecPointsMul", "RecPoints multiplicity distribution in EMCAL", 500, 0, 1000) ;
  22   cM->Sumw2() ;
  23
  24   AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
  25   if (rl == 0x0)
  26     cout<<"Can not instatiate the Run Loader"<<endl;
  27
  28   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
  29     (rl->GetDetectorLoader("EMCAL"));
  30
  31   //Load Hits
  32   rl->LoadHits("EMCAL");
  33   //Load Digits
  34   rl->LoadDigits("EMCAL");
  35   //Load SDigits
  36   rl->LoadSDigits("EMCAL");
  37   //Load RecPoints
  38   rl->LoadRecPoints("EMCAL");
  39
  40
  41   //one way to do it that works for all three
  42   AliEMCALHit* hit;
  43   AliEMCALDigit* dig;
  44   AliEMCALRecPoint* rp;
  45
  46   rl->GetEvent(0);
  47
  48   TClonesArray *hits = 0;
  49   //Fill array of digits
  50   TClonesArray *sdigits = emcalLoader->SDigits();
  51   //Fill array of digits
  52   TClonesArray *digits = emcalLoader->Digits();
  53   //Fill array of clusters
  54   //TObjArray *clusters = emcalLoader->RecPoints(); //It should work, need to FIX
  55   TObjArray *clusters = 0;
  56
  57   //Get hits from the list
  58
  59   //Hits are stored in different branches in the hits Tree,
  60   //first get the branch and then access the hits in the branch
  61   Int_t nHit = 0;
  62   TTree *treeH = emcalLoader->TreeH();
  63   if (treeH) {
  64     // TreeH exists, get the branch
  65     Int_t nTrack = treeH->GetEntries();  // TreeH has array of hits for every primary
  66     TBranch * branchH = treeH->GetBranch("EMCAL");
  67     branchH->SetAddress(&hits);
  68     //Now get the hits in this branch
  69     for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
  70       branchH->GetEntry(iTrack);
  71       Int_t nHit = hits->GetEntriesFast();
  72       for(Int_t ihit = 0; ihit< nHit;ihit++){
  73         hit = static_cast<AliEMCALHit *>hits->At(ihit);
  74         if(hit != 0){
  75           nHit++;
  76           hE->Fill(hit->GetEnergy());
  77           cout<<"Hit Info "<<hit->GetId()<<" ELoss "<<hit->GetEnergy()<<endl;
  78         }//hit?
  79       }//hit loop
  80     }// track loop
  81   }//treeH?
  82
  83   hM->Fill(nHit);
  84
  85   //Get digits from the list
  86   if(sdigits){
  87     sM->Fill(sdigits->GetEntries());
  88     for(Int_t isdig = 0; isdig< sdigits->GetEntries();isdig++){
  89       //cout<<">> idig "<<idig<<endl;
  90       dig = static_cast<AliEMCALDigit *>(sdigits->At(isdig)) ;
  91
  92       if(dig != 0){
  93         sE->Fill(dig->GetAmplitude());
  94         cout << "SDigit info " << dig->GetId() << " " << dig->GetAmplitude() << endl;
  95       }
  96     }
  97   }
  98   else printf("No Sdigits available\n");
  99
 100
 101   //Get digits from the list
 102   if(digits){
 103     dM->Fill(digits->GetEntries());
 104     for(Int_t idig = 0; idig< digits->GetEntries();idig++){
 105       //cout<<">> idig "<<idig<<endl;
 106       dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
 107
 108       if(dig != 0){
 109         dA->Fill(dig->GetAmplitude());
 110         cout << "Digit info " << dig->GetId() << " " << dig->GetAmplitude() << endl;
 111       }
 112     }
 113   }
 114   else printf("No digits available\n");
 115
 116   //Get clusters from the list
 117   TTree *treeR = emcalLoader->TreeR();
 118   TBranch * branchR = treeR->GetBranch("EMCALECARP");
 119   branchR->SetAddress(&clusters);
 120   branchR->GetEntry(0);
 121
 122   if(clusters){
 123     cM->Fill(clusters->GetEntries());
 124     for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
 125       //cout<<">> idig "<<idig<<endl;
 126       rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
 127
 128       if(rp != 0){
 129         cE->Fill(rp->GetEnergy());
 130         cout << "RecPoint info " << rp->GetAbsId(0) << " " << rp->GetEnergy() << endl;
 131       }
 132     }
 133   }else printf("No recpoints available\n");
 134
 135   /*
 136    //another way to do it
 137    TTree *hTree = rl->GetTreeH("EMCAL",false);
 138    TTree *dTree = rl->GetTreeD("EMCAL",false);
 139    TTree *sTree = rl->GetTreeS("EMCAL",false);
 140    TTree *cTree = rl->GetTreeR("EMCAL",false);
 141
 142    TObjArray *harr=NULL;
 143    TBranch *hbranch=hTree->GetBranch("EMCAL");
 144    hbranch->SetAddress(&harr);
 145
 146    TObjArray *darr=NULL;
 147    TBranch *dbranch=dTree->GetBranch("EMCAL");
 148    dbranch->SetAddress(&darr);
 149
 150    TObjArray *sarr=NULL;
 151    TBranch *sbranch=sTree->GetBranch("EMCAL");
 152    sbranch->SetAddress(&sarr);
 153
 154    TObjArray *carr=NULL;
 155    TBranch *cbranch=cTree->GetBranch("EMCALECARP");
 156    cbranch->SetAddress(&carr);
 157
 158
 159    if(hbranch->GetEvent(0)) {
 160    for(Int_t ih = 0; ih < harr->GetEntriesFast(); ih++) {
 161    hM->Fill(harr->GetEntriesFast());
 162    AliEMCALHit* hit =(AliEMCALHit*)harr->UncheckedAt(ih);
 163    if(hit != 0){
 164          hE->Fill(hit->GetEnergy());
 165          cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy()*10.5 << endl;
 166    }
 167    }
 168    }
 169
 170    if(dbranch->GetEvent(0)) {
 171    for(Int_t id = 0; id < darr->GetEntriesFast(); id++) {
 172    dM->Fill(darr->GetEntriesFast());
 173    AliEMCALDigit* dig =(AliEMCALDigit*)darr->UncheckedAt(id);
 174    if(dig != 0){
 175          dA->Fill(dig->GetAmp());
 176          cout << "Digit info " << dig->GetId() << " " << dig->GetAmp() << endl;
 177    }
 178    }
 179    }
 180
 181    if(sbranch->GetEvent(0)) {
 182    for(Int_t id = 0; id < sarr->GetEntriesFast(); id++) {
 183    sM->Fill(sarr->GetEntriesFast());
 184    AliEMCALDigit* sdig =(AliEMCALDigit*)sarr->UncheckedAt(id);
 185    if(sdig != 0){
 186          sE->Fill(sdig->GetAmp()/1.e+6);
 187          cout << "SDigit info " << sdig->GetId() << " " << sdig->GetAmp()/1.e+6 << endl;
 188    }
 189    }
 190    }
 191
 192    if(cbranch->GetEvent(0)) {
 193    for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
 194    cE->Fill(carr->GetEntriesFast());
 195    AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
 196    if(rp != 0){
 197          cE->Fill(rp->GetEnergy());
 198          cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
 199    }
 200    }
 201    }
 202
 203    */
 204
 205   TCanvas *chits = new TCanvas("chits","Hits",20,20,800,400);
 206   chits->Divide(2,1);
 207   chits->cd(1);
 208   hE->Draw();
 209   chits->cd(2);
 210   hM->Draw();
 211
 212   TCanvas *cdig = new TCanvas("cdig","Digits",20,40,800,400);
 213   cdig->Divide(2,1);
 214   cdig->cd(1);
 215   dA->Draw();
 216   cdig->cd(2);
 217   dM->Draw();
 218
 219
 220   TCanvas *csdig = new TCanvas("csdig","SDigits",20,60,800,400);
 221   csdig->Divide(2,1);
 222   csdig->cd(1);
 223   sE->Draw();
 224   csdig->cd(2);
 225   sM->Draw();
 226
 227   TCanvas *cclu = new TCanvas("cclu","Clusters",20,60,800,400);
 228   cclu->Divide(2,1);
 229   cclu->cd(1);
 230   cE->Draw();
 231   cclu->cd(2);
 232   cM->Draw();
 233
 234 }