HMPID/Hqa.C

   1 void Hqa()
   2 {
   3   gROOT->LoadMacro("Hdisp.C");
   4   HitQA();
   5 }//Hqa()
   6 void HitQA()
   7 {
   8   TFile *pFile = new TFile("Hqa.root","recreate");
   9   AliHMPIDCluster::DoCorrSin(kFALSE);
  10   AliHMPIDDigit::fSigmas=4;
  11   Int_t nEvts=10000;
  12   TLegend *lQ=new TLegend(0.5,0.5,0.9,0.9);
  13
  14   TH1F *hQ7eV  =new TH1F("hQ7eV"  ,"" ,300,-50,2000);     hQ7eV  ->SetLineColor(kRed);      lQ->AddEntry(hQ7eV  ,"Ckov 7 eV");  hQ7eV->SetStats(0);
  15   TH1F *hQ200eV=new TH1F("hQ200eV","" ,300,-50,2000);     hQ200eV->SetLineColor(kBlack);    lQ->AddEntry(hQ200eV,"mip 200 eV");
  16   TH1F *hQ500eV=new TH1F("hQ500eV","" ,300,-50,2000);     hQ500eV->SetLineColor(kCyan);     lQ->AddEntry(hQ500eV,"mip 500 eV");
  17   TH1F *hQ900eV=new TH1F("hQ900eV","" ,300,-50,2000);     hQ900eV->SetLineColor(kGreen);    lQ->AddEntry(hQ900eV,"mip 900 eV");
  18
  19   TH1F *hCluPerEvt=new TH1F("hCluPerEvt","# clusters per event",11,-0.5,10.5);
  20   TH1F *hCluChi2  =new TH1F("hChi2","Chi2 ",1000,0,100);
  21   TH1F *hCluFlg   =new TH1F("hCluFlg","Cluster flag",14,-1.5,12.5);                       hCluFlg->SetFillColor(5);
  22   TH1F *hCluRawSize= new TH1F("hCluRawSize","Raw cluster size ",100,0,100);
  23
  24   gStyle->SetOptStat(10);
  25   TH2F *pCluMapSi1  =new TH2F("cluMapSi1","Size 1 map"       ,1700,-10,160,1700,-10,160);
  26   TH2F *pCluMapLo0  =new TH2F("cluMNoLo0","Loc Max 0 map"    ,1700,-10,160,1700,-10,160);
  27   TH2F *pCluMapLo1  =new TH2F("cluMapLo1","Loc Max 1 map"    ,1700,-10,160,1700,-10,160);
  28   TH2F *pCluMapUnf  =new TH2F("cluMapUnf","Unfolded map"     ,1700,-10,160,1700,-10,160);
  29   TH2F *pCluMapEdg  =new TH2F("cluMapEdg","On edge map"      ,1700,-10,160,1700,-10,160);
  30   TH2F *pCluMapCoG  =new TH2F("cluMapCoG","CoG  map"         ,1700,-10,160,1700,-10,160);
  31   TH2F *pCluMapEmp  =new TH2F("cluMapEmp","undefined-empty"  ,1700,-10,160,1700,-10,160);
  32   TH2F *pCluMapNoLoc=new TH2F("cluMapNoLoc","no loc maxima"  ,1700,-10,160,1700,-10,160);
  33   TH2F *pCluMapAbn  =new TH2F("cluMapAbn","abnormal fit"     ,1700,-10,160,1700,-10,160);
  34   TH2F *pCluMapNot  =new TH2F("cluMapNot","Raw Clusters"     ,1700,-10,160,1700,-10,160);
  35   TH2F *pCluMapMax  =new TH2F("cluMapMax","N. locs excceds"  ,1700,-10,160,1700,-10,160);
  36
  37   TH1F *hHitCluDifX = new TH1F("hHitCluDifX" ,";entries;x_{Hit}-x_{Clu} [cm]"   ,1000,-1,1);          hHitCluDifX->Sumw2();    hHitCluDifX->SetFillColor(kYellow);
  38 //  TH2F *hHitCluDifXv= new TH2F("hHitCluDifXv",";x_{Hit};x_{Hit}-x_{Clu} [cm]"   ,500,-0.5,0.5,1000,-0.2,0.2);hHitCluDifXv->Sumw2();
  39   TProfile *hHitCluDifXv= new TProfile("hHitCluDifXv",";x_{Hit};x_{Hit}-x_{Clu} [cm]"   ,500,-0.5,0.5);
  40   TH1F *hHitCluDifY = new TH1F("hHitCluDifY" ,";entries;y_{Hit}-y_{Clu} [cm]"   ,1000,-1,1);          hHitCluDifY->Sumw2();    hHitCluDifY->SetFillColor(kYellow);
  41   TH2F *hHitCluDifXY= new TH2F("hHitCluDifXY",";x_{Hit}-x_{Clu};y_{Hit}-y_{Clu}",1000,-1,1,1000,-1,1);hHitCluDifXY->Sumw2();
  42   TH1F *hHitCluDifQ = new TH1F("hHitCluDifQ" ,";entries;(Q_{Clu}-Q_{Hit})/Q_{Hit}" ,200 ,-200,200);   hHitCluDifQ->Sumw2();    hHitCluDifQ->SetFillColor(kYellow);
  43
  44   TH2F *hHitMap= new TH2F("hHitMap",";x_{Hit};y_{Hit}",1700,-10,160,1700,-10,160);
  45
  46   Float_t e200=200e-9,e500=500e-9,e900=900e-9,e7=7e-9;//predefined  Eloss
  47
  48   AliHMPIDHit hit(0,0,kProton,0,0,0);
  49   for(int i=0;i<5000;i++){
  50     hQ200eV->Fill(hit.QdcTot(e200));  hQ500eV->Fill(hit.QdcTot(e500));   hQ900eV->Fill(hit.QdcTot(e900));  hQ7eV->Fill(hit.QdcTot(e7));
  51   }
  52   TClonesArray hits("AliHMPIDHit");  TClonesArray sdigs("AliHMPIDDigit");
  53   TObjArray digs(7); for(Int_t i=0;i<7;i++) digs.AddAt(new TClonesArray("AliHMPIDDigit"),i);
  54   TObjArray clus(7); for(Int_t i=0;i<7;i++) clus.AddAt(new TClonesArray("AliHMPIDCluster"),i);
  55
  56
  57   for(Int_t iEvt=0;iEvt<nEvts;iEvt++){//events loop
  58     if(iEvt%50==0)Printf("============> iEvt = %d ",iEvt);
  59
  60     Int_t ch,pid; Float_t e,hitx,hity,hitq;
  61 //    Int_t nHits=(type==999)?1:40;
  62     Int_t nHits=1;
  63     for(Int_t iHit=0;iHit<nHits;iHit++){//hits loop for the current event
  64       switch(iHit){
  65         case 0:  ch=0;pid=kProton;e=e200;
  66                  hitx=AliHMPIDDigit::SizePadX()*(6+gRandom->Rndm());
  67                  hity=AliHMPIDDigit::SizePadY()*(6+gRandom->Rndm());break; //mip ramdomly distributed in one pad
  68         case 1:  ch=0;pid=kProton;e=e200;hitx=0.4;hity=0.42;break; //mip in left-hand bottom coner of chamber 0
  69         case 2:  ch=0;pid=kProton;e=e200;hitx=0.4;hity=30  ;break; //mip on left edge of chamber 0
  70         case 3:  ch=0;pid=kProton;e=e200;hitx=40; hity=0.42;break; //mip on bottom edge of chamber 0
  71         default: ch=gRandom->Rndm()*6; pid=(gRandom->Rndm()>0.9)? kProton:50000050;
  72                   if(pid==kProton)
  73                     e=gRandom->Rndm()*900e-9;
  74                   else
  75                     e=5.5e-9+3e-9*gRandom->Rndm();
  76                   hitx=gRandom->Rndm()*AliHMPIDDigit::SizeAllX(); hity=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();break; //random hit
  77       }
  78       new(hits[iHit]) AliHMPIDHit(ch,e,pid,iHit,hitx,hity);
  79       hitq=e;
  80     }//hits loop
  81
  82     AliHMPIDv1::Hit2Sdi(&hits,&sdigs);
  83     AliHMPIDDigitizer::DoNoise(kFALSE);
  84     AliHMPIDDigitizer::Sdi2Dig(&sdigs,&digs);
  85     AliHMPIDReconstructor::Dig2Clu(&digs,&clus,kTRUE);
  86
  87 // From here normal procedure for QA
  88
  89     for(Int_t iHit=0;iHit<hits.GetEntriesFast();iHit++) {
  90       AliHMPIDHit *pHit = (AliHMPIDHit*)hits.UncheckedAt(iHit);
  91       hHitMap->Fill(pHit->LorsX(),pHit->LorsY());
  92     }
  93     Int_t kMaxCh=(nHits==1)?0:AliHMPIDDigit::kMaxCh;
  94     for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=kMaxCh;iCh++){//chambers loop
  95       TClonesArray *pDigs=(TClonesArray *)digs.UncheckedAt(iCh);
  96       TClonesArray *pClus=(TClonesArray *)clus.UncheckedAt(iCh);
  97
  98 //      if(pClus->GetEntriesFast()>nHits) {hits.Print();}
  99 //      if(pClus->GetEntriesFast()>nHits) {pDigs->Print();Printf("----S D I G I T S-------------");}
 100 //      if(pClus->GetEntriesFast()>nHits) {sdigs.Print();Printf("-------------------------------");}
 101       hCluPerEvt->Fill(pClus->GetEntriesFast());
 102       for(Int_t iClu=0;iClu<pClus->GetEntriesFast();iClu++){//clusters loop
 103         AliHMPIDCluster *pClu=(AliHMPIDCluster*)pClus->UncheckedAt(iClu);
 104         Float_t clux=pClu->X(); Float_t cluy=pClu->Y(); Float_t cluq=pClu->Q();
 105         hCluFlg->Fill(pClu->Status());
 106         hCluChi2->Fill(pClu->Chi2());
 107         hCluRawSize->Fill(pClu->Size());
 108
 109         switch(pClu->Status()){
 110 /*
 111     case        kFrm  : status="formed        "   ;break;
 112     case        kUnf  : status="unfolded (fit)"   ;break;
 113     case        kCoG  : status="coged         "   ;break;
 114     case        kLo1  : status="locmax 1 (fit)"   ;break;
 115     case        kMax  : status="exceeded (cog)"   ;break;
 116     case        kNot  : status="not done (cog)"   ;break;
 117     case        kEmp  : status="empty         "   ;break;
 118     case        kEdg  : status="edge     (fit)"   ;break;
 119     case        kSi1  : status="size 1   (cog)"   ;break;
 120     case        kNoLoc: status="no LocMax(fit)"   ;break;
 121     case        kAbn  : status="Abnormal fit  "   ;break;
 122 */
 123             case AliHMPIDCluster::kSi1:   pCluMapSi1->Fill(clux,cluy);  break;
 124             case AliHMPIDCluster::kLo1:   pCluMapLo1->Fill(clux,cluy);  break;
 125             case AliHMPIDCluster::kUnf:   pCluMapUnf->Fill(clux,cluy);  break;
 126             case AliHMPIDCluster::kMax:   pCluMapMax->Fill(clux,cluy);  break;
 127             case AliHMPIDCluster::kEdg:   pCluMapEdg->Fill(clux,cluy);  break;
 128             case AliHMPIDCluster::kCoG:   pCluMapCoG->Fill(clux,cluy);  break;
 129             case AliHMPIDCluster::kNoLoc: pCluMapNoLoc->Fill(clux,cluy);break;
 130             case AliHMPIDCluster::kAbn:   pCluMapAbn->Fill(clux,cluy);  break;
 131             case AliHMPIDCluster::kNot:   pCluMapNot->Fill(clux,cluy);  break;
 132             default:     pCluMapEmp->Fill(clux,cluy); break;
 133         }
 134
 135         hHitCluDifX->Fill(hitx-clux); hHitCluDifY->Fill(hity-cluy); hHitCluDifXY->Fill(hitx-clux,hity-cluy); hHitCluDifQ->Fill(100*(cluq-hitq)/hitq);
 136         // distorsion due to feedback photons
 137         Int_t pc,px,py;
 138         AliHMPIDDigit::Lors2Pad(hitx,hity,pc,px,py);
 139         Float_t padCenterX = AliHMPIDDigit::LorsX(pc,px);
 140         if(pClu->Size()>1)hHitCluDifXv->Fill(hitx-padCenterX,(hitx-clux));
 141         //
 142       }//clusters loop
 143     }//chambers loop
 144
 145     hits.Delete();  sdigs.Delete();  for(int i = 0;i<7;i++){((TClonesArray*)digs.At(i))->Delete();((TClonesArray*)clus.At(i))->Delete();}
 146   }//events loop
 147
 148   gStyle->SetPalette(1);
 149   TCanvas *pC2=new TCanvas("Digit canvas","Digit canvas",1280,800); pC2->Divide(3,3);
 150   pC2->cd(1);hHitCluDifX->Draw("hist");
 151   pC2->cd(2);hHitCluDifY->Draw("hist");
 152 //  pC2->cd(3);hHitCluDifXY->Draw("colz");
 153   // Draw CorrSin
 154   AliHMPIDCluster c;
 155   TPolyLine *pLine = new TPolyLine(500);
 156   for(Int_t i=0;i<500;i++) {
 157     Double_t x = 0 + i*AliHMPIDDigit::SizePadX()/500.;
 158     c.SetX(x);c.SetY(0);c.CorrSin();
 159     pLine->SetPoint(i,x-0.5*AliHMPIDDigit::SizePadX(),c.X()-x);
 160   }
 161   pC2->cd(3);hHitCluDifXv->Draw("colz");pLine->Draw("C");
 162   //
 163   pC2->cd(4);hHitCluDifQ->Draw("hist");
 164   pC2->cd(5);gPad->SetLogy(1);hCluFlg->Draw();
 165   pC2->cd(6);hCluChi2->Draw();
 166   pC2->cd(7);                 hCluRawSize->Draw();
 167   pC2->cd(8);                 hCluPerEvt->Draw();
 168   pC2->cd(9);                 hQ7eV->Draw(); hQ200eV->Draw("same"); hQ500eV->Draw("same"); hQ900eV->Draw("same"); lQ->Draw();
 169   TCanvas *pC1=new TCanvas("ClusterMaps","Cluster maps",1280,800); pC1->Divide(3,3);
 170   pC1->cd(1);  pCluMapSi1->Draw();  DrawCh(-1);
 171   pC1->cd(2);  pCluMapLo1->Draw();  DrawCh(-1);
 172   pC1->cd(3);  pCluMapUnf->Draw();  DrawCh(-1);
 173   pC1->cd(4);     hHitMap->Draw();  DrawCh(-1);
 174   pC1->cd(5);  pCluMapMax->Draw();  DrawCh(-1);
 175   pC1->cd(6);  pCluMapEdg->Draw();  DrawCh(-1);
 176   pC1->cd(7);  pCluMapNot->Draw();  DrawCh(-1);
 177   pC1->cd(8);  pCluMapNoLoc->Draw();DrawCh(-1);
 178   pC1->cd(9);  pCluMapEmp->Draw();  DrawCh(-1);
 179
 180   pC1->SaveAs("$HOME/HitMaps.png");  //?????
 181   pC2->SaveAs("$HOME/HitCluDif.gif");
 182
 183   pFile->Write();
 184   Printf("Digits - raw -digits conversion...");
 185
 186   AliHMPIDDigit d1,d2; Int_t ddl,r,d,a;UInt_t w32;
 187
 188   for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=AliHMPIDDigit::kMaxCh;iCh++)
 189   for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++)
 190   for(Int_t iPx=AliHMPIDDigit::kMinPx;iPx<=AliHMPIDDigit::kMaxPx;iPx++)
 191   for(Int_t iPy=AliHMPIDDigit::kMinPy;iPy<=AliHMPIDDigit::kMaxPy;iPy++){
 192     d1.Set(iCh,iPc,iPx,iPy,3040);   //set digit
 193     d1.Raw(w32,ddl,r,d,a);          //get raw word for this digit
 194     d2.Raw(w32,ddl);                //set another digit from that raw word
 195     if(d1.Compare(&d2)) {d1.Print(); d2.Print(); Printf("");}//compare them
 196   }
 197   Printf("OK");
 198 }//tst()