HMPID/Hqa.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TSystem.h>
   3 #include <TFile.h>
   4 #include <TStyle.h>
   5 #include <TTree.h>
   6 #include <TClonesArray.h>
   7 #include <TObjArray.h>
   8 #include <TF1.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TProfile.h>
  12 #include <TCanvas.h>
  13
  14 #include <AliHMPIDParam.h>
  15 #include <AliHMPIDHit.h>
  16 #include <AliHMPIDCluster.h>
  17 #include <AliHMPIDDigit.h>
  18 #endif
  19
  20 Int_t nEntries = 0;
  21
  22 TObjArray *CreateContainer(const char *classname,TTree *pTree);
  23 void Hits(Int_t mode,TTree *pTree=0x0);
  24 void Digs(Int_t mode, TTree *pTree=0x0);
  25 void Clus(Int_t mode, TTree *pTree=0x0);
  26 void Summary();
  27
  28 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  29 TObjArray *CreateContainer(const char *classname,TTree *pTree)
  30 {
  31   TObjArray *pOA=new TObjArray(AliHMPIDParam::kMaxCh+1);
  32   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
  33     TClonesArray *pCA=new TClonesArray(classname);
  34     pOA->AddAt(pCA,iCh);
  35   }
  36
  37   pOA->SetOwner(kTRUE);
  38
  39   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
  40     pTree->SetBranchAddress(Form("HMPID%i",iCh),&(*pOA)[iCh]);
  41   }
  42   return pOA;
  43 }
  44 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  45 TH1F *hHitQdc,*hHitQdcCh[AliHMPIDParam::kMaxCh+1]; TH2F *hHitMap[AliHMPIDParam::kMaxCh+1];
  46 Double_t fHitMean[AliHMPIDParam::kMaxCh+1],fHitErr[AliHMPIDParam::kMaxCh+1];
  47 void Hits(Int_t mode,TTree *pTree)
  48 {
  49   switch(mode){
  50     case 1:
  51       hHitQdc=new TH1F("HitQdc","Hit Qdc all chamber;QDC",4000,0,4000);
  52   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
  53         hHitMap[iCh]=new TH2F(Form("HitMap%i",iCh),Form("Ch%i;x_{Hit};y_{Hit}",iCh),
  54             Int_t(AliHMPIDParam::SizeAllX()/AliHMPIDParam::SizePadX()),0,AliHMPIDParam::SizeAllX(),
  55             Int_t(AliHMPIDParam::SizeAllY()/AliHMPIDParam::SizePadY()),0,AliHMPIDParam::SizeAllY());
  56         hHitQdcCh[iCh]=new TH1F(Form("HMPID%i",iCh),Form("Charge for HMPID%i",iCh),4000,0,4000);
  57       }
  58       break;
  59     case 2:
  60       if(pTree==0) return;
  61       TClonesArray *pHits=new TClonesArray("AliHMPIDHit");  pTree->SetBranchAddress("HMPID",&pHits);
  62       for(Int_t iEnt=0;iEnt<pTree->GetEntriesFast();iEnt++){//entries loop
  63         pTree->GetEntry(iEnt);
  64         for(Int_t iHit=0;iHit<pHits->GetEntriesFast();iHit++){//hits loop
  65           AliHMPIDHit *pHit = (AliHMPIDHit*)pHits->UncheckedAt(iHit);
  66           hHitMap[pHit->Ch()]->Fill(pHit->LorsX(),pHit->LorsY());
  67           hHitQdc->Fill(pHit->Q());
  68           hHitQdcCh[pHit->Ch()]->Fill(pHit->Q());
  69         }//hits loop
  70       }//entries loop
  71       delete pHits;
  72       break;
  73     case 3:
  74       TCanvas *c1=new TCanvas("HitCan","Hits",1280,800); c1->Divide(3,3);
  75
  76       for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
  77         if(iCh==6) c1->cd(1); if(iCh==5) c1->cd(2);
  78         if(iCh==4) c1->cd(4); if(iCh==3) c1->cd(5); if(iCh==2) c1->cd(6);
  79                               if(iCh==1) c1->cd(8); if(iCh==0) c1->cd(9);
  80         gStyle->SetPalette(1);
  81         hHitMap[iCh]->Draw("colz");
  82       }
  83       for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
  84         fHitMean[iCh] = 0;
  85         fHitErr[iCh]  = 0;
  86         if((Int_t)hHitQdcCh[iCh]->GetEntries()<nEntries) continue;
  87         c1->cd(3);hHitQdcCh[iCh]->Fit("expo","Q");
  88         TF1 *funcfit = (TF1*)hHitQdcCh[iCh]->FindObject("expo");
  89         fHitMean[iCh] = funcfit->GetParameter(1);
  90         fHitErr[iCh]  = funcfit->GetParError(1);
  91       }
  92       TPad *pad = (TPad*)c1->cd(3); hHitQdc->SetFillColor(5); pad->SetLogy(); hHitQdc->Fit("expo","Q");
  93       break;
  94   }
  95 }//Hits()
  96 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  97 TH1F     *hDigQ;
  98 TProfile *hDigHighQ,*hDigChEvt;
  99 void Digs(Int_t mode, TTree *pTree)
 100 {
 101   switch(mode){
 102     case 1:
 103       hDigHighQ=new TProfile("hDigHighQ","Highest charge in chamber  ",AliHMPIDParam::kMaxCh+1,AliHMPIDParam::kMinCh,AliHMPIDParam::kMaxCh+1);
 104       hDigChEvt=new TProfile("hDigChEvt","Chamber occupancy per event",AliHMPIDParam::kMaxCh+1,AliHMPIDParam::kMinCh,AliHMPIDParam::kMaxCh+1);
 105       hDigQ    =new TH1F("hDigQ        ","Charge of digits (ADC)     ",3000,0,3000);
 106       break;
 107     case 2:
 108       if(pTree==0) return;
 109       TObjArray *pLst=CreateContainer("AliHMPIDDigit",pTree); pTree->GetEntry(0);
 110       for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){//chambers loop
 111         TClonesArray *pDigs=(TClonesArray *)pLst->UncheckedAt(iCh);
 112         hDigChEvt->Fill(iCh,pDigs->GetEntriesFast()/(48.*80.*6.)*100.);
 113         Double_t highQ = 0;
 114         for(Int_t iDig=0;iDig<pDigs->GetEntriesFast();iDig++){//digits loop
 115           AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigs->UncheckedAt(iDig);
 116           hDigQ->Fill(pDig->Q());
 117           if(pDig->Q()>highQ) highQ = pDig->Q();
 118         }
 119         hDigHighQ->Fill(iCh,highQ);
 120       }
 121       delete pLst;
 122       break;
 123     case 3:
 124       TCanvas *c1=new TCanvas("DigQa","Digit Check",1280,800); c1->Divide(2,2);
 125       gStyle->SetOptFit(1);
 126       TPad *pad = (TPad*)c1->cd(1); pad->SetLogy(); hDigQ->Fit("expo","QN");
 127       c1->cd(2); hDigHighQ->Draw();
 128       c1->cd(3); hDigChEvt->Draw();
 129       break;
 130   }//switch
 131 }//Dig()
 132 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 133 TH1F *hCluEvt,*hCluChi2,*hCluFlg,*hCluSize,*hCluQ;
 134 void Clus(Int_t mode, TTree *pTree)
 135 {
 136   switch(mode){
 137     case 1:
 138       hCluEvt=new TH1F("CluPerEvt","Cluster multiplicity"   ,100,0,100);
 139       hCluChi2  =new TH1F("CluChi2"  ,"Chi2 "               ,1000,0,100);
 140       hCluFlg   =new TH1F("CluFlg"   ,"Cluster flag"        ,14,-1.5,12.5);                       hCluFlg->SetFillColor(5);
 141       hCluSize  =new TH1F("CluSize"  ,"Cluster size        ",100,0,100);
 142       hCluQ     =new TH1F("CluQ"     ,"Cluster charge (ADC)",1000,0,5000);
 143       break;
 144     case 2:
 145       if(pTree==0) return;
 146
 147       TObjArray *pLst=CreateContainer("AliHMPIDCluster",pTree); pTree->GetEntry(0);
 148       for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){//chambers loop
 149         TClonesArray *pClus=(TClonesArray *)pLst->UncheckedAt(iCh);
 150         hCluEvt->Fill(pClus->GetEntriesFast());
 151         for(Int_t iClu=0;iClu<pClus->GetEntriesFast();iClu++){//clusters loop
 152           AliHMPIDCluster *pClu=(AliHMPIDCluster*)pClus->UncheckedAt(iClu);
 153           hCluFlg->Fill(pClu->Status());
 154           hCluChi2->Fill(pClu->Chi2());
 155           hCluSize->Fill(pClu->Size());
 156           hCluQ->Fill(pClu->Q());
 157         }
 158       }
 159       delete pLst;
 160       break;
 161     case 3:
 162       TCanvas *c1=new TCanvas("CluComCan","Clusters in common",1280,800); c1->Divide(3,3);
 163       c1->cd(1); hCluEvt->SetFillColor(5);      hCluEvt->Draw();
 164       c1->cd(2); hCluChi2->SetFillColor(5);     hCluChi2->Draw();
 165       c1->cd(3); hCluFlg->SetFillColor(5);      hCluFlg->Draw();
 166       c1->cd(4); hCluSize->SetFillColor(5);     hCluSize->Draw();
 167       TPad *pad = (TPad*)c1->cd(5); hCluQ->SetFillColor(5); pad->SetLogy(); hCluQ->Draw();
 168       break;
 169   }//switch
 170 }//Clus()
 171 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 172 void Hqa()
 173 {
 174   AliHMPIDParam::Instance();
 175
 176   TFile *fh=0; if(gSystem->IsFileInIncludePath("HMPID.Hits.root"))      fh=TFile::Open("HMPID.Hits.root"     ,"read");if(fh->IsZombie()) fh=0;
 177   TFile *fd=0; if(gSystem->IsFileInIncludePath("HMPID.Digits.root"))    fd=TFile::Open("HMPID.Digits.root"   ,"read");if(fd->IsZombie()) fd=0;
 178   TFile *fc=0; if(gSystem->IsFileInIncludePath("HMPID.RecPoints.root")) fc=TFile::Open("HMPID.RecPoints.root","read");if(fc->IsZombie()) fc=0;
 179   if(fh==0 && fd==0 && fc==0){Printf("Nothing to do!"); return;}
 180   if(fh) Hits(1); if(fc) Clus(1);  if(fd) Digs(1);//book
 181   Int_t iEvt=0;
 182   while(1){
 183     TTree *th=0; if(fh) th=(TTree*)fh->Get(Form("Event%i/TreeH",iEvt));
 184     TTree *td=0; if(fd) td=(TTree*)fd->Get(Form("Event%i/TreeD",iEvt));
 185     TTree *tc=0; if(fc) tc=(TTree*)fc->Get(Form("Event%i/TreeR",iEvt));
 186
 187     Hits(2,th); Digs(2,td); Clus(2,tc); //fill
 188     if(th==0 && td==0 && tc==0) break;
 189     iEvt++;
 190     if(!(iEvt%50)) Printf("Event %i processed",iEvt);
 191   }
 192
 193   nEntries = iEvt;
 194
 195   if(fd) Clus(3);//plot everything
 196   if(fc) Digs(3);
 197   if(fh) Hits(3);
 198   Summary();
 199 }
 200 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 201 void Summary()
 202 {
 203   //info for hits...
 204   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
 205     Printf(" #################### Summary for HMPID %i#################### ",iCh);
 206     //info for hits...
 207     Printf("-----Summary of Hits----- ");
 208     if(fHitMean[iCh]==0) {
 209       Printf("gain %5.2f +/- %5.2f",fHitMean[iCh],fHitErr[iCh]);
 210     } else {
 211       Double_t gain = 1./TMath::Abs(fHitMean[iCh]);
 212       Double_t errgain = gain*gain*fHitErr[iCh];
 213       Printf("gain %5.2f +/- %5.2f",gain,errgain);
 214     }
 215     //info for digits...
 216     Printf("-----Summary of Digits-----");
 217     Printf(" Chamber %d with occupancy (%) %5.2f +/- %5.2f",iCh,hDigChEvt->GetBinContent(iCh+1),hDigChEvt->GetBinError(iCh+1));
 218     Printf(" Chamber %d with higest Q (ADC) %7.0f +/- %7.0f",iCh,hDigHighQ->GetBinContent(iCh+1),hDigHighQ->GetBinError(iCh+1));
 219   }
 220 }