[u/mrichter/AliRoot.git] / HMPID / Hdisp.C

TCanvas *pAll=0;
AliRunLoader *gAL=0; AliLoader *gHL=0; AliESD *gEsd=0; TTree *gEsdTr=0; AliHMPID *gH=0;
Int_t gEvt=0; Int_t gMaxEvt=0;
TObjArray *pNmean;
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void Hdisp()
{//display events from files if any in current directory or simulated events
  pAll=new TCanvas("all","",1300,900); pAll->Divide(3,3,0,0);
//  pAll->ToggleEditor();
  pAll->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)",0,"","DoZoom(Int_t,Int_t,Int_t,TObject*)");

  OpenCalib();  
  if(gSystem->IsFileInIncludePath("galice.root")){// tries to open session
    if(gAlice) delete gAlice;                                               //in case we execute this in aliroot delete default AliRun object 
    gAL=AliRunLoader::Open();                                                                    //try to open galice.root from current dir 
    gAL->LoadgAlice();                                                                           //take new AliRun object from galice.root   
    gHL=gAL->GetDetectorLoader("HMPID");  gH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");  //get HMPID object from galice.root
    gMaxEvt=gAL->GetNumberOfEvents()-1;
    gHL->LoadHits(); gHL->LoadSDigits(); gHL->LoadDigits(); gHL->LoadRecPoints();

    TFile *pEsdFl=TFile::Open("AliESDs.root"); gEsdTr=(TTree*) pEsdFl->Get("esdTree"); gEsdTr->SetBranchAddress("ESD", &gEsd);
    pAll->cd(7); TButton *pBtn=new TButton("Next","ReadEvt()",0,0,0.2,0.1);   pBtn->Draw();
                 TButton *pHitBtn=new TButton("Print hits","PrintHits()",0,0.2,0.3,0.3);   pHitBtn->Draw();
                 TButton *pSdiBtn=new TButton("Print sdis","PrintSdis()",0,0.4,0.3,0.5);   pSdiBtn->Draw();
                 TButton *pDigBtn=new TButton("Print digs","PrintDigs()",0,0.6,0.3,0.7);   pDigBtn->Draw();
                 TButton *pCluBtn=new TButton("Print clus","PrintClus()",0,0.8,0.3,0.9);   pCluBtn->Draw();  
    ReadEvt();
  }else{
    pAll->cd(7); TButton *pBtn=new TButton("Next","SimEvt()",0,0,0.2,0.1);   pBtn->Draw(); 
    SimEvt();
  }      
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void ReadEvt()
{// Read curent event and display it assumes that session is alredy opened
  TClonesArray hits("AliHMPIDHit");
  if(gEvt>gMaxEvt) gEvt=0; if(gEvt<0) gEvt=gMaxEvt;
  
  gEsdTr->GetEntry(gEvt); gAL->GetEvent(gEvt); 
  ReadHits(&hits); gHL->TreeS()->GetEntry(0); gHL->TreeD()->GetEntry(0); gHL->TreeR()->GetEntry(0);
  
  
  DrawEvt(&hits,gH->DigLst(),gH->CluLst(),gEsd);
  gEvt++;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void SimEvt()
{
  TClonesArray hits("AliHMPIDHit"); 
  TClonesArray sdig("AliHMPIDDigit"); 
  TObjArray    digs(7); for(Int_t i=0;i<7;i++) digs.AddAt(new TClonesArray("AliHMPIDDigit"),i);
  TObjArray    clus(7); for(Int_t i=0;i<7;i++) clus.AddAt(new TClonesArray("AliHMPIDCluster"),i);
  AliESD esd;
  AliHMPIDDigit::fSigmas=4;
  AliHMPIDDigitizer::DoNoise(kFALSE);
  SimEsd(&esd);
  SimHits(&esd,&hits);
             AliHMPIDv1::Hit2Sdi(&hits,&sdig);                               
      AliHMPIDDigitizer::Sdi2Dig(&sdig,&digs);     
      AliHMPIDReconstructor::Dig2Clu(&digs,&clus);
      AliHMPIDTracker::Recon(&esd,&clus,pNmean);
  
  pAll->cd(3);  gPad->Clear(); TLatex txt;txt.DrawLatex(0.2,0.2,Form("Simulated event %i",gEvt));
  DrawEvt(&hits,&digs,&clus,&esd);  
  gEvt++;
}//SimEvt()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void SimEsd(AliESD *pEsd)
{
  TParticle part; TLorentzVector mom;
  for(Int_t iTrk=0;iTrk<100;iTrk++){//stack loop
    part.SetPdgCode(kProton);
    part.SetProductionVertex(0,0,0,0);
    Double_t eta= -0.2+gRandom->Rndm()*0.4; //rapidity is random [-0.2,+0.2]
    Double_t phi= gRandom->Rndm()*60.*TMath::DegToRad();   //phi is random      [ 0  , 60 ] degrees    
    mom.SetPtEtaPhiM(5,eta,phi,part.GetMass());
    part.SetMomentum(mom);
    AliESDtrack trk(&part);
    pEsd->AddTrack(&trk);
  }//stack loop  
}//EsdFromStack()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void SimHits(AliESD *pEsd, TClonesArray *pHits)
{
  const Int_t kCerenkov=50000050;
  const Int_t kFeedback=50000051;
  
  AliHMPIDRecon rec;
  
  Int_t hc=0; TVector2 pos;
  for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop
    AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
    Float_t xRa,yRa;
    Int_t ch=AliHMPIDTracker::IntTrkCha(pTrk,xRa,yRa);
    if(ch<0) continue; //this track does not hit HMPID
    Float_t beta = pTrk->GetP()/(TMath::Sqrt(pTrk->GetP()*pTrk->GetP()+0.938*0.938));
    Float_t ckov=TMath::ACos(1./(beta*1.292));

    Float_t theta,phi,xPc,yPc,; 
    pTrk->GetHMPIDtrk(xPc,yPc,theta,phi); 
    rec.SetTrack(xRa,yRa,theta,phi); 
    
    if(!AliHMPIDDigit::IsInDead(xPc,yPc)) new((*pHits)[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,iTrk,xPc,yPc);                 //mip hit
    /*
    for(int i=0;i<4;i++) { 
      Float_t x=gRandom->Rndm()*130;Float_t y=gRandom->Rndm()*126;
      if(!AliHMPIDDigit::IsInDead(x,y)) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,x,y); //bkg hits 4 per track
    }
    */
    Int_t nPhots = (Int_t)(20.*TMath::Power(TMath::Sin(ckov),2)/TMath::Power(TMath::Sin(TMath::ACos(1./1.292)),2));
    for(int i=0;i<nPhots;i++){
      rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
      if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());
    }                      //photon hits  
  }//tracks loop    
}//SimHits()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void ReadHits(TClonesArray *pHitLst)
{
  pHitLst->Delete();  Int_t cnt=0;
  for(Int_t iEnt=0;iEnt<gHL->TreeH()->GetEntries();iEnt++){       //TreeH loop
    gHL->TreeH()->GetEntry(iEnt);                                 //get current entry (prim)                
    for(Int_t iHit=0;iHit<gH->Hits()->GetEntries();iHit++){       //hits loop
      AliHMPIDHit *pHit = (AliHMPIDHit*)gH->Hits()->At(iHit);     //get current hit        
      new((*pHitLst)[cnt++]) AliHMPIDHit(*pHit);
    }//hits loop for this entry
  }//tree entries loop
  
}//ReadHits()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void DrawCh(Int_t iCh) 
{ 
  gPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5); 
  if(iCh>=0){TLatex txt; txt.SetTextSize(0.1); txt.DrawLatex(-5,-5,Form("%i",iCh));}
  
  for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++){
    TBox *pBox=new TBox(AliHMPIDDigit::fMinPcX[iPc],AliHMPIDDigit::fMinPcY[iPc],
                        AliHMPIDDigit::fMaxPcX[iPc],AliHMPIDDigit::fMaxPcY[iPc]);
    
    pBox->SetFillStyle(0);
    pBox->Draw();
    
  }//PC loop      
}//DrawCh()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void DrawEvt(TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
{//draws all the objects of current event in given canvas
  AliHMPIDRecon rec;  
  TPolyMarker *pTxC[7], *pRin[7]; TMarker *pMip,*pCko,*pFee,*pDig,*pClu;
  pMip=new TMarker; pMip->SetMarkerColor(kRed);  pMip->SetMarkerStyle(kOpenTriangleUp);
  pCko=new TMarker; pCko->SetMarkerColor(kRed);  pCko->SetMarkerStyle(kOpenCircle);
  pFee=new TMarker; pFee->SetMarkerColor(kRed);  pFee->SetMarkerStyle(kOpenDiamond);
  pDig=new TMarker; pDig->SetMarkerColor(kGreen);pDig->SetMarkerStyle(kOpenSquare);
  pClu=new TMarker; pClu->SetMarkerColor(kBlue); pClu->SetMarkerStyle(kStar);
  for(Int_t ch=0;ch<7;ch++){
    pTxC[ch]=new TPolyMarker; pTxC[ch]->SetMarkerStyle(2); pTxC[ch]->SetMarkerColor(kRed); pTxC[ch]->SetMarkerSize(3);
    pRin[ch]=new TPolyMarker; pRin[ch]->SetMarkerStyle(6); pRin[ch]->SetMarkerColor(kMagenta);
  }
  
  for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop to collect cerenkov rings and intersection points
    AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
    Int_t ch=pTrk->GetHMPIDcluIdx();
    if(ch<0) continue; //this track does not hit HMPID
    ch/=1000000; 
    Float_t th,ph,xPc,yPc; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  //get info on current track
    pTxC[ch]->SetNextPoint(xPc,yPc);                          //add new intersection point
    
    Float_t ckov=pTrk->GetHMPIDsignal();  Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
    if(ckov>0){
      Printf("theta %f phi %f ckov %f",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),ckov);
      rec.SetTrack(xPc,yPc,th,ph+TMath::Pi());
      TVector2 pos;
       Double_t allGapz=rec.fgkWinThick+0.5*rec.fgkRadThick+rec.fgkGapThick;   //to semplify, the (x,y) are calculted from (xPc,yPc) back to radiator in straight line (Bz=0)
      TVector3 dir(0,0,-1);TVector3 miprad(xPc,yPc,allGapz);
      AliHMPIDRecon::Propagate(dir,miprad,0);
      Double_t xRad=miprad.X();
      Double_t yRad=miprad.Y();
      for(int j=0;j<100;j++){ 
        rec.TracePhot(xRad,yRad,ckov,j*0.0628,pos);
       if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) pRin[ch]->SetNextPoint(pos.X(),pos.Y());
      }      
    }
  }//tracks loop
      
  Int_t totHit=pHitLst->GetEntriesFast(),totDig=0,totClu=0,totTxC=0;
  Int_t totCkov=0, totFeed=0,totMip=0;
  for(Int_t iCh=0;iCh<7;iCh++){//chambers loop    
    totTxC+=pTxC[iCh]->GetN();
    totDig+=((TClonesArray*)pDigLst->At(iCh))->GetEntriesFast();
    totClu+=((TClonesArray*)pCluLst->At(iCh))->GetEntriesFast();
    
    switch(iCh){
      case 6: pAll->cd(1); break; case 5: pAll->cd(2); break;
      case 4: pAll->cd(4); break; case 3: pAll->cd(5); break; case 2: pAll->cd(6); break;
                                  case 1: pAll->cd(8); break; case 0: pAll->cd(9); break;
    }
    gPad->SetEditable(kTRUE); gPad->Clear(); 
    DrawCh(iCh);
                           
    ((TClonesArray*)pDigLst->At(iCh))->Draw();               //draw digits
    
    totCkov=0;totFeed=0;totMip=0;
    for(Int_t iHit=0;iHit<pHitLst->GetEntriesFast();iHit++){ // Draw hits
      AliHMPIDHit *pHit=(AliHMPIDHit*)pHitLst->At(iHit);
      if(pHit->Ch()==iCh) pHit->Draw();
      if(pHit->Pid()==50000050) totCkov++;
      else if(pHit->Pid()==50000051) totFeed++;
      else totMip++;
    }    
           
    ((TClonesArray*)pCluLst->At(iCh))->Draw();              //draw clusters
                            pTxC[iCh]->Draw();              //draw intersections
                            pRin[iCh]->Draw("CLP");         //draw rings
    gPad->SetEditable(kFALSE);
  }//chambers loop
  
  pAll->cd(3);  gPad->Clear(); TLegend *pLeg=new TLegend(0.2,0.2,0.8,0.8);
                                        pLeg->SetHeader(Form("Event %i Total %i",gEvt,gMaxEvt+1));
                                        pLeg->AddEntry(pTxC[0],Form("TRKxPC %i"     ,totTxC),"p");
                                        pLeg->AddEntry(pMip   ,Form("Mip hits %i"   ,totMip),"p");    
                                        pLeg->AddEntry(pCko   ,Form("Ckov hits %i"  ,totCkov),"p");    
                                        pLeg->AddEntry(pFee   ,Form("Feed hits %i"  ,totFeed),"p");    
                                        pLeg->AddEntry(pDig   ,Form("Digs %i"       ,totDig),"p");    
                                        pLeg->AddEntry(pClu   ,Form("Clus %i"       ,totClu),"p");    
                                        pLeg->Draw();
}//DrawEvt()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void DoZoom(Int_t evt, Int_t px, Int_t py, TObject *obj)
{
  if(evt!=5 && evt!=6) return; //5- zoom in 6-zoom out
  const Int_t minZoom=64;
  const Int_t maxZoom=2;
  static Int_t zoom=minZoom; //zoom level
  if(evt==5&&zoom==maxZoom) return; 
  if(evt==6&&zoom==minZoom) return; 
  
 // if(!obj->IsA()->InheritsFrom("TPad")) return;  //current object is not pad
  TVirtualPad *pPad=gPad->GetSelectedPad();
  if(pPad->GetNumber()==3 || pPad->GetNumber()==7) return; //current pad is wrong

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