HMPID/Hdisp.C

   1 TCanvas *pAll=0;
   2 AliRunLoader *gAL=0; AliLoader *gHL=0; AliESD *gEsd=0; TTree *gEsdTr=0; AliHMPID *gH=0;
   3 Int_t gEvt=0; Int_t gMaxEvt=0;
   4 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   5 void Hdisp()
   6 {//display events from files if any in current directory or simulated events
   7   pAll=new TCanvas("all","",1000,900); pAll->Divide(3,3,0,0);pAll->ToggleEditor();
   8   pAll->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)",0,"","DoZoom(Int_t,Int_t,Int_t,TObject*)");
   9
  10   if(gSystem->IsFileInIncludePath("galice.root")){// tries to open session
  11     if(gAlice) delete gAlice;                                               //in case we execute this in aliroot delete default AliRun object
  12     gAL=AliRunLoader::Open();                                                                    //try to open galice.root from current dir
  13     gAL->LoadgAlice();                                                                           //take new AliRun object from galice.root
  14     gHL=gAL->GetDetectorLoader("HMPID");  gH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");  //get HMPID object from galice.root
  15     gMaxEvt=gAL->GetNumberOfEvents()-1;
  16     gHL->LoadHits(); gHL->LoadSDigits(); gHL->LoadDigits(); gHL->LoadRecPoints();
  17
  18     TFile *pEsdFl=TFile::Open("AliESDs.root"); gEsdTr=(TTree*) pEsdFl->Get("esdTree"); gEsdTr->SetBranchAddress("ESD", &gEsd);
  19     pAll->cd(7); TButton *pBtn=new TButton("Next","ReadEvt()",0,0,0.2,0.1);   pBtn->Draw();
  20                  TButton *pHitBtn=new TButton("Print hits","PrintHits()",0,0.2,0.3,0.3);   pHitBtn->Draw();
  21                  TButton *pSdiBtn=new TButton("Print sdis","PrintSdis()",0,0.4,0.3,0.5);   pSdiBtn->Draw();
  22                  TButton *pDigBtn=new TButton("Print digs","PrintDigs()",0,0.6,0.3,0.7);   pDigBtn->Draw();
  23                  TButton *pCluBtn=new TButton("Print clus","PrintClus()",0,0.8,0.3,0.9);   pCluBtn->Draw();
  24     ReadEvt();
  25   }else{
  26     pAll->cd(7); TButton *pBtn=new TButton("Next","SimEvt()",0,0,0.2,0.1);   pBtn->Draw();
  27     SimEvt();
  28   }
  29 }
  30 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  31 void ReadEvt()
  32 {// Read curent event and display it assumes that session is alredy opened
  33   TClonesArray hits("AliHMPIDHit");
  34   if(gEvt>gMaxEvt) gEvt=0; if(gEvt<0) gEvt=gMaxEvt;
  35
  36   gEsdTr->GetEntry(gEvt); gAL->GetEvent(gEvt);
  37   ReadHits(&hits); gHL->TreeS()->GetEntry(0); gHL->TreeD()->GetEntry(0); gHL->TreeR()->GetEntry(0);
  38
  39
  40   DrawEvt(&hits,gH->DigLst(),gH->CluLst(),gEsd);
  41   gEvt++;
  42 }
  43 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  44 void SimEvt()
  45 {
  46   TClonesArray hits("AliHMPIDHit");
  47   TClonesArray sdig("AliHMPIDDigit");
  48   TObjArray    digs(7); for(Int_t i=0;i<7;i++) digs.AddAt(new TClonesArray("AliHMPIDDigit"),i);
  49   TObjArray    clus(7); for(Int_t i=0;i<7;i++) clus.AddAt(new TClonesArray("AliHMPIDCluster"),i);
  50   AliESD esd;
  51
  52   SimEsd(&esd);
  53   SimHits(&esd,&hits);
  54              AliHMPIDv1::Hit2Sdi(&hits,&sdig);
  55       AliHMPIDDigitizer::Sdi2Dig(&sdig,&digs);
  56   AliHMPIDReconstructor::Dig2Clu(&digs,&clus);
  57         AliHMPIDTracker::Recon(&esd,&clus);
  58
  59   pAll->cd(3);  gPad->Clear(); TLatex txt;txt.DrawLatex(0.2,0.2,Form("Simulated event %i",gEvt));
  60   DrawEvt(&hits,&digs,&clus,&esd);
  61   gEvt++;
  62 }//SimEvt()
  63 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  64 void SimEsd(AliESD *pEsd)
  65 {
  66   TParticle part; TLorentzVector mom;
  67   for(Int_t iTrk=0;iTrk<25;iTrk++){//stack loop
  68     part.SetPdgCode(kProton);
  69     part.SetProductionVertex(0,0,0,0);
  70     Double_t eta= -0.4+gRandom->Rndm()*0.8; //rapidity is random [-0.4,+0.4]
  71     Double_t phi= gRandom->Rndm()*1.4;      //phi is random      [ 0  , 80 ] degrees
  72     mom.SetPtEtaPhiM(2,eta,phi,part.GetMass());
  73     part.SetMomentum(mom);
  74     AliESDtrack trk(&part);
  75     pEsd->AddTrack(&trk);
  76   }//stack loop
  77 }//EsdFromStack()
  78 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  79 void SimHits(AliESD *pEsd, TClonesArray *pHits)
  80 {
  81   AliHMPIDRecon rec;
  82   const Int_t kCerenkov=50000050,kFeedback=50000051;
  83   Int_t hc=0; TVector2 pos;
  84   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop
  85     AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
  86     Float_t xRa,yRa;
  87     Int_t ch=AliHMPIDTracker::IntTrkCha(pTrk,xRa,yRa);
  88     if(ch<0) continue; //this track does not hit HMPID
  89     Float_t ckov=0.63;
  90
  91     Float_t th,ph,xPc,yPc,; pTrk->GetHMPIDtrk(xPc,yPc,th,ph); rec.SetTrack(xRa,yRa,th,ph);
  92
  93     if(!AliHMPIDDigit::IsInDead(xPc,yPc)) new((*pHits)[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,iTrk,xPc,yPc);                 //mip hit
  94     for(int i=0;i<4;i++)  new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,gRandom->Rndm()*130,gRandom->Rndm()*126); //bkg hits 4 per track
  95     for(int i=0;i<16;i++){
  96       rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
  97       new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());
  98     }                      //photon hits
  99   }//tracks loop
 100 }//SimHits()
 101 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 102 void ReadHits(TClonesArray *pHitLst)
 103 {
 104   pHitLst->Delete();  Int_t cnt=0;
 105   for(Int_t iEnt=0;iEnt<gHL->TreeH()->GetEntries();iEnt++){       //TreeH loop
 106     gHL->TreeH()->GetEntry(iEnt);                                 //get current entry (prim)
 107     for(Int_t iHit=0;iHit<gH->Hits()->GetEntries();iHit++){       //hits loop
 108       AliHMPIDHit *pHit = (AliHMPIDHit*)gH->Hits()->At(iHit);     //get current hit
 109       new((*pHitLst)[cnt++]) AliHMPIDHit(*pHit);
 110     }//hits loop for this entry
 111   }//tree entries loop
 112
 113 }//ReadHits()
 114 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 115 void DrawCh(Int_t iCh)
 116 {
 117   gPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5);
 118   TLatex txt; txt.SetTextSize(0.1); txt.DrawLatex(-5,-5,Form("%i",iCh));
 119
 120   for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++){
 121     TBox *pBox=new TBox(AliHMPIDDigit::fMinPcX[iPc],AliHMPIDDigit::fMinPcY[iPc],
 122                         AliHMPIDDigit::fMaxPcX[iPc],AliHMPIDDigit::fMaxPcY[iPc]);
 123
 124     pBox->SetFillStyle(0);
 125     pBox->Draw();
 126
 127   }//PC loop
 128 }//DrawCh()
 129 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 130 void DrawEvt(TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
 131 {//draws all the objects of current event in given canvas
 132
 133   AliHMPIDRecon rec;
 134   TPolyMarker *pTxC[7], *pRin[7]; TMarker *pMip,*pCko,*pFee,*pDig,*pClu;
 135   pMip=new TMarker; pMip->SetMarkerColor(kRed);  pMip->SetMarkerStyle(kOpenTriangleUp);
 136   pCko=new TMarker; pCko->SetMarkerColor(kRed);  pCko->SetMarkerStyle(kOpenCircle);
 137   pFee=new TMarker; pFee->SetMarkerColor(kRed);  pFee->SetMarkerStyle(kOpenDiamond);
 138   pDig=new TMarker; pDig->SetMarkerColor(kGreen);pDig->SetMarkerStyle(kOpenSquare);
 139   pClu=new TMarker; pClu->SetMarkerColor(kBlue); pClu->SetMarkerStyle(kStar);
 140   for(Int_t ch=0;ch<7;ch++){
 141     pTxC[ch]=new TPolyMarker; pTxC[ch]->SetMarkerStyle(2); pTxC[ch]->SetMarkerColor(kRed); pTxC[ch]->SetMarkerSize(3);
 142     pRin[ch]=new TPolyMarker; pRin[ch]->SetMarkerStyle(6); pRin[ch]->SetMarkerColor(kMagenta);
 143   }
 144
 145   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop to collect cerenkov rings and intersection points
 146     AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
 147     Int_t ch=pTrk->GetHMPIDcluIdx();
 148     if(ch<0) continue; //this track does not hit HMPID
 149     ch/=1000000;
 150     Float_t th,ph,xPc,yPc; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  //get info on current track
 151     pTxC[ch]->SetNextPoint(xPc,yPc);                          //add new intersection point
 152
 153     Float_t ckov=pTrk->GetHMPIDsignal();  Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
 154     if(ckov>0){
 155       rec.SetTrack(xPc,yPc,th,ph);
 156      TVector2 pos;  for(int j=0;j<100;j++){rec.TracePhot(ckov,j*0.0628,pos); pRin[ch]->SetNextPoint(pos.X(),pos.Y());}
 157     }
 158   }//tracks loop
 159
 160   Int_t totHit=pHitLst->GetEntriesFast(),totDig=0,totClu=0,totTxC=0;
 161   for(Int_t iCh=0;iCh<7;iCh++){//chambers loop
 162     totTxC+=pTxC[iCh]->GetN();
 163     totDig+=((TClonesArray*)pDigLst->At(iCh))->GetEntriesFast();
 164     totClu+=((TClonesArray*)pCluLst->At(iCh))->GetEntriesFast();
 165
 166     switch(iCh){
 167       case 6: pAll->cd(1); break; case 5: pAll->cd(2); break;
 168       case 4: pAll->cd(4); break; case 3: pAll->cd(5); break; case 2: pAll->cd(6); break;
 169                                   case 1: pAll->cd(8); break; case 0: pAll->cd(9); break;
 170     }
 171     gPad->SetEditable(kTRUE); gPad->Clear();
 172     DrawCh(iCh);
 173
 174     for(Int_t iHit=0;iHit<pHitLst->GetEntriesFast();iHit++){
 175       AliHMPIDHit *pHit=(AliHMPIDHit*)pHitLst->At(iHit);
 176       if(pHit->Ch()==iCh) pHit->Draw();
 177     }
 178
 179     ((TClonesArray*)pDigLst->At(iCh))->Draw();  //draw digits
 180     ((TClonesArray*)pCluLst->At(iCh))->Draw();  //draw clusters
 181                             pTxC[iCh]->Draw();  //draw intersections
 182                             pRin[iCh]->Draw();  //draw rings
 183     gPad->SetEditable(kFALSE);
 184   }//chambers loop
 185
 186   pAll->cd(3);  gPad->Clear(); TLegend *pLeg=new TLegend(0.2,0.2,0.8,0.8);
 187                                         pLeg->SetHeader(Form("Event %i Total %i",gEvt,gMaxEvt+1));
 188                                         pLeg->AddEntry(pTxC[0],Form("TRKxPC %i"     ,totTxC),"p");
 189                                         pLeg->AddEntry(pMip   ,Form("Mip hits %i"   ,totHit),"p");
 190                                         pLeg->AddEntry(pCko   ,Form("Ckov hits %i"  ,totHit),"p");
 191                                         pLeg->AddEntry(pFee   ,Form("Feed hits %i"  ,totHit),"p");
 192                                         pLeg->AddEntry(pDig   ,Form("Digs %i"       ,totDig),"p");
 193                                         pLeg->AddEntry(pClu   ,Form("Clus %i"       ,totClu),"p");
 194                                         pLeg->Draw();
 195 }//DrawEvt()
 196 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 197 void DoZoom(Int_t evt, Int_t px, Int_t py, TObject *obj)
 198 {
 199   if(evt!=5 && evt!=6) return; //5- zoom in 6-zoom out
 200   const Int_t minZoom=64;
 201   const Int_t maxZoom=2;
 202   static Int_t zoom=minZoom; //zoom level
 203   if(evt==5&&zoom==maxZoom) return;
 204   if(evt==6&&zoom==minZoom) return;
 205
 206  // if(!obj->IsA()->InheritsFrom("TPad")) return;  //current object is not pad
 207   TVirtualPad *pPad=gROOT->GetSelectedPad();
 208   if(pPad->GetNumber()==3 || pPad->GetNumber()==7) return; //current pad is wrong
 209
 210  // Printf("evt=%i (%i,%i) %s",evt,px,py,obj->GetName());
 211
 212   Float_t x=pPad->AbsPixeltoX(px); Float_t y=pPad->AbsPixeltoY(py);
 213
 214   if(evt==5){ zoom=zoom/2;     pPad->Range(x-zoom*2,y-zoom*2,x+zoom*2,y+zoom*2);} //zoom in
 215   else      { zoom=zoom*2;     pPad->Range(x-zoom*2,y-zoom*2,x+zoom*2,y+zoom*2);} //zoom out
 216   if(zoom==minZoom) pPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5);
 217   ((TCanvas *)gTQSender)->SetTitle(Form("zoom x%i",minZoom/zoom));
 218   pPad->Modified();
 219   pPad->Update();
 220 }
 221 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 222 void PrintHits()
 223 {
 224 //Prints a list of HMPID hits for a given event. Default is event number 0.
 225   Printf("List of HMPID hits for event %i",gEvt);
 226
 227   Int_t iTotHits=0;
 228   for(Int_t iPrim=0;iPrim<gHL->TreeH()->GetEntries();iPrim++){//prims loop
 229     gHL->TreeH()->GetEntry(iPrim);
 230     gH->Hits()->Print();
 231     iTotHits+=gH->Hits()->GetEntries();
 232   }
 233   Printf("totally %i hits for event %i",iTotHits,gEvt);
 234 }
 235 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 236 void PrintSdis()
 237 {
 238 //prints a list of HMPID sdigits  for a given event
 239   Printf("List of HMPID sdigits for event %i",gEvt);
 240   gH->SdiLst()->Print();
 241   Printf("totally %i sdigits for event %i",gH->SdiLst()->GetEntries(),gEvt);
 242 }
 243 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 244 void PrintDigs()
 245 {
 246 //prints a list of HMPID digits
 247   Printf("List of HMPID digits for event %i",gEvt);
 248   gH->DigLst()->Print();
 249   Int_t totDigs=0;  for(Int_t i=0;i<7;i++) {totDigs+=gH->DigLst(i)->GetEntries();}
 250   Printf("totally %i digits for event %i",totDigs,gEvt);
 251 }
 252 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 253 void PrintClus()
 254 {//prints a list of HMPID clusters  for a given event
 255   Printf("List of HMPID clusters for event %i",gEvt);
 256   gH->CluLst()->Print();
 257
 258   Int_t iCluCnt=0; for(Int_t iCh=0;iCh<7;iCh++) iCluCnt+=gH->CluLst(iCh)->GetEntries();
 259
 260   Printf("totally %i clusters for event %i",iCluCnt,gEvt);
 261 }