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

#include <TF1.h>
#include <TH2F.h>
#include <TCanvas.h>  //Terminate()
#include <TChain.h>
#include <TBenchmark.h>
#include <TFile.h>    //docosmic()    
#include <AliSelector.h>      //base class
#include <AliESDEvent.h>           
#include <AliBitPacking.h> //HmpidPayload()
#include "AliHMPIDDigit.h" 
#include "AliHMPIDParam.h" 
#include "AliHMPIDCluster.h" 
#include "AliHMPIDReconstructor.h" //docosmic()

class AliHMPIDSelector : public AliSelector {
 public :
           AliHMPIDSelector():AliSelector(),fChain(0),fEsd(0),fCkovP(0),fMipXY(0),fDifX(0),fSigP(0) {for(Int_t i=0;i<5;i++) fProb[i]=0;}
  virtual ~AliHMPIDSelector()                                                                      {delete fEsd;}


  virtual Int_t   Version        () const {return 2;}
  virtual void    Begin          (TTree *) {}
  virtual void    SlaveBegin     (TTree *tree);
  virtual void    Init           (TTree *tree);
  virtual Bool_t  Notify         ()    {return kTRUE;}
  virtual Bool_t  Process        (Long64_t entry);
  virtual void    SetOption      (const char *option) { fOption = option; }
  virtual void    SetObject      (TObject *obj) { fObject = obj; }
  virtual void    SetInputList   (TList *input) {fInput = input;}
  virtual TList  *GetOutputList  () const { return fOutput; }
  virtual void    SlaveTerminate ();
  virtual void    Terminate      ();

 private: 
  TTree          *fChain ;   //!pointer to the analyzed TTree or TChain
  AliESDEvent    *fEsd ;     //!

  TH2F           *fCkovP,*fMipXY;                //!
  TH1F           *fDifX;                         //!
  TH2F           *fSigP;                         //!
  TH1F           *fProb[5];                      //!

  ClassDef(AliHMPIDSelector,0);  
};


//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDSelector::SlaveBegin(TTree *tree)
{
// The SlaveBegin() function is called after the Begin() function.  When running with PROOF SlaveBegin() is called on each slave server.
// The tree argument is deprecated (on PROOF 0 is passed).

   Init(tree);

   TString option = GetOption();

   // create histograms on each slave server
   fCkovP    = new TH2F("CkovP" , "#theta_{c}, [rad];P, [GeV]", 150,   0,  7  ,500, 0, 1); 
   fSigP     = new TH2F("SigP"  ,"#sigma_{#theta_c}"          , 150,   0,  7  ,100, 0, 1e20);
   fMipXY    = new TH2F("MipXY" ,"mip position"               , 260,   0,130  ,252,0,126); 
   fDifX     = new TH1F("DifX" ,"diff"                       , 100, -2.5, 2.5);
      
   fProb[0] = new TH1F("PidE" ,"PID: e yellow #mu magenta"  ,100,0,1); fProb[0]->SetLineColor(kYellow);
   fProb[1] = new TH1F("PidMu","pid of #mu"                 ,100,0,1); fProb[1]->SetLineColor(kMagenta);
   fProb[2] = new TH1F("PidPi","PID: #pi red K green p blue",100,0,1); fProb[2]->SetLineColor(kRed);
   fProb[3] = new TH1F("PidK" ,"pid of K"                   ,100,0,1); fProb[3]->SetLineColor(kGreen);
   fProb[4] = new TH1F("PidP" ,"pid of p"                   ,100,0,1); fProb[4]->SetLineColor(kBlue);
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDSelector::Init(TTree *pTr)
{
// Called when the selector needs to initialize a new tree of chain. Typically here the branch addresses of the tree is set
// Init() will be called many times when running with PROOF.

  if ( !pTr )  return ;
  fChain = pTr ;
  fChain->SetBranchAddress("ESD", &fEsd) ;
  fChain->SetBranchStatus("*", 1);
  fChain->SetBranchStatus("fTracks.*", 1);
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Bool_t AliHMPIDSelector::Process(Long64_t entry)
{
  AliHMPIDParam *pParam=AliHMPIDParam::Instance();
  fChain->GetTree()->GetEntry(entry);
 
  for(Int_t iTrk=0;iTrk<fEsd->GetNumberOfTracks();iTrk++){
    AliESDtrack *pTrk=fEsd->GetTrack(iTrk);
     
//     if(pTrk->GetHMPIDsignal()==-1) continue;
     
    fCkovP->Fill(pTrk->GetP(),pTrk->GetHMPIDsignal()) ; 
    fSigP ->Fill(pTrk->GetP(),TMath::Sqrt(pTrk->GetHMPIDchi2()));
     
    Float_t xClu,yClu; Int_t q,np;  
    Float_t xPc,yPc;  
    pTrk->GetHMPIDmip(xClu,yClu,q,np);  
    fMipXY->Fill(xClu,yClu); //mip info
    Float_t xRad,yRad,th,ph;        
    pTrk->GetHMPIDtrk(xRad,yRad,th,ph); 
    Int_t iCh=pTrk->GetHMPIDcluIdx();iCh/=1000000;
    Double_t p1[3],n1[3]; pParam->Norm(iCh,n1); pParam->Lors2Mars(iCh,0,0,p1,AliHMPIDParam::kPc);      //point & norm  for PC
    if(pTrk->Intersect(p1,n1,fEsd->GetMagneticField())==kFALSE) continue;                              //try to intersect track with PC
    pParam->Mars2Lors   (iCh,p1,xPc,yPc);                                                              //TRKxPC position
    fDifX->Fill(xPc-xClu); //track info 
     
    Double_t pid[5];  pTrk->GetHMPIDpid(pid); for(Int_t i =0;i<5;i++) fProb[i]->Fill(pid[i]);
  }//tracks loop 
     
  return kTRUE;
}//Process()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDSelector::SlaveTerminate()
{
  // The SlaveTerminate() function is called after all entries or objects
  // have been processed. When running with PROOF SlaveTerminate() is called
  // on each slave server.
  
  // Add the histograms to the output on each slave server
  
  fOutput->Add(fCkovP);
  fOutput->Add(fSigP); 
  fOutput->Add(fMipXY);
  fOutput->Add(fDifX);
  
  for(Int_t i=0;i<5;i++) fOutput->Add(fProb[i]);
}//SlaveTerminate()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDSelector::Terminate()
{
// The last function to be called. It always runs on the client, it can be used to present
// the results graphically or save the results to file.
  fCkovP   = dynamic_cast<TH2F*>(fOutput->FindObject("CkovP")) ;
  fSigP    = dynamic_cast<TH2F*>(fOutput->FindObject("SigP")) ; 
  fMipXY   = dynamic_cast<TH2F*>(fOutput->FindObject("MipXY")) ;
  fDifX    = dynamic_cast<TH1F*>(fOutput->FindObject("DifX")) ;
  
  fProb[0] = dynamic_cast<TH1F*>(fOutput->FindObject("PidE")) ;
  fProb[1] = dynamic_cast<TH1F*>(fOutput->FindObject("PidMu")) ;
  fProb[2] = dynamic_cast<TH1F*>(fOutput->FindObject("PidPi")) ;
  fProb[3] = dynamic_cast<TH1F*>(fOutput->FindObject("PidK")) ;
  fProb[4] = dynamic_cast<TH1F*>(fOutput->FindObject("PidP")) ;

  Float_t n=1.292; //mean freon ref idx 
  TF1 *pPi=new TF1("RiPiTheo","acos(sqrt(x*x+[0]*[0])/(x*[1]))",1.2,7); pPi->SetLineWidth(1); pPi->SetParameter(1,n); 
  AliPID ppp;                 pPi->SetLineColor(kRed);   pPi->SetParameter(0,AliPID::ParticleMass(AliPID::kPion));    //mass
  TF1 *pK=(TF1*)pPi->Clone(); pK ->SetLineColor(kGreen); pK ->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon)); 
  TF1 *pP=(TF1*)pPi->Clone(); pP ->SetLineColor(kBlue);  pP ->SetParameter(0,AliPID::ParticleMass(AliPID::kProton)); 

  TCanvas *pC=new TCanvas("c1","ESD QA");pC->SetFillColor(10); pC->SetHighLightColor(10); pC->Divide(3,2);
  pC->cd(1); fCkovP->Draw(); pPi->Draw("same"); pK->Draw("same"); pP->Draw("same");   pC->cd(2); fMipXY->Draw();  pC->cd(3); fProb[0]->Draw(); fProb[1]->Draw("same"); 
  pC->cd(4); fSigP ->Draw();                                                          pC->cd(5); fDifX->Draw();   pC->cd(6); fProb[2]->Draw(); fProb[3]->Draw("same"); fProb[4]->Draw("same"); 
}//Terminate()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void loc()
{
  TChain* pChain =new TChain("esdTree");  pChain->Add("AliESDs.root");

  pChain->Process("AliHMPIDSelector.C+");	
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Int_t DateHeader(ifstream *pFile,Bool_t isPrint=0)
{
  Int_t iSize=-1;
  pFile->read((char*)&iSize,4);
  if(!isPrint)     
    pFile->seekg(16*4,ios::cur);
  else{
    Int_t w32=-1; 
                                Printf("");
                                Printf("Event size        %i bytes",iSize);                           //1  common DATE header 17 words
    pFile->read((char*)&w32,4); Printf("Event magic       0x%x"    ,w32);                             //2
    pFile->read((char*)&w32,4); Printf("Event head size   %i bytes",w32);                             //3  
    pFile->read((char*)&w32,4); Printf("Event version     0x%x"    ,w32);                             //4
    pFile->read((char*)&w32,4); Printf("Event type        %i (%s)" ,w32,(w32==7)? "physics":"SOR");   //5 
    pFile->read((char*)&w32,4); Printf("Run number        %i"      ,w32);                             //6
  
    pFile->read((char*)&w32,4); Printf("Event ID 1        %i"      ,w32);                             //7
    pFile->read((char*)&w32,4); Printf("Event ID 2        %i"      ,w32);                             //8
  
    pFile->read((char*)&w32,4); Printf("Trigger pattern 1 %i"      ,w32);                             //9
    pFile->read((char*)&w32,4); Printf("Trigger pattern 2 %i"      ,w32);                             //10
  
    pFile->read((char*)&w32,4); Printf("Detector pattern  %i"      ,w32);                             //11
  
    pFile->read((char*)&w32,4); Printf("Type attribute 1  %i"      ,w32);                             //12
    pFile->read((char*)&w32,4); Printf("Type attribute 2  %i"      ,w32);                             //13
    pFile->read((char*)&w32,4); Printf("Type attribute 3  %i"      ,w32);                             //14
  
    pFile->read((char*)&w32,4); Printf("LDC ID            %i"      ,w32);                             //15
    pFile->read((char*)&w32,4); Printf("GDC ID            %i"      ,w32);                             //16
    pFile->read((char*)&w32,4); TDatime time(w32); time.Print();                                      //17
  
    Printf("");
  }
  return iSize;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void HmpidHeader(ifstream *pFile,Bool_t isPrint=kFALSE)
{
// Prints hmpid trailer and returns number of words for this trailer  
  if(!isPrint) {pFile->seekg(15*4,ios::cur);return;}
  Int_t w32=-1;
  Printf("\nHMPID Header:");//private HEADER is 15 words
  for(Int_t i=1;i<=11;i++) { pFile->read((char*)&w32,4); Printf("Word #%2i=%12i meaningless",i,w32);}
                             pFile->read((char*)&w32,4); Printf("Word #12=%12i event counter",w32);   
  for(Int_t i=13;i<=15;i++){ pFile->read((char*)&w32,4); Printf("Word #%2i=%12i empty",i,w32);}
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void HmpidPayload(ifstream *pFile,Int_t iDdl,TObjArray *pDigAll)
{
// payload is 8 sequences with structure: WC36A8 then WC number of w32
  
  TClonesArray *pDig1=(TClonesArray*)pDigAll->At(iDdl/2); //get list of digits for requested chamber
  
  UInt_t w32=0;
  Int_t iDigCnt=pDig1->GetEntriesFast();
  for(Int_t row=1;row<=8;row++){  
    pFile->read((char*)&w32,4);  Int_t wc=AliBitPacking::UnpackWord(w32,16,31); Int_t ma=AliBitPacking::UnpackWord(w32, 0,15);    
    if(ma!=0x36a8) Printf("ERROR ERROR ERROR WRONG Marker=0x%x ERROR ERROR ERROR",ma);    
    for(Int_t i=1;i<=wc;i++){//words loop
      pFile->read((char*)&w32,4);
      if(w32&0x08000000) continue; //it's DILOGIC CW
      AliHMPIDDigit *pDig=new AliHMPIDDigit;
      pDig->Raw(w32,iDdl);
      new ((*pDig1)[iDigCnt++]) AliHMPIDDigit(*pDig);
    }//words loop 
  }//rows loop
}//HmpidPayload()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void docosmic(const char* name,Int_t iMaxEvt=9999999)
{
  gBenchmark->Start("Cosmic converter");
  ifstream in(name);

  
  Bool_t isPrint=kFALSE;
  
  TString rooName=name; rooName.Replace(1+rooName.First('.'),4,"root");
  Int_t ddl=0;    
  
  TFile *pOut=new TFile(rooName,"recreate");
  TTree *pTr=new TTree("cosmic","some for time being");
  
  
  TObjArray *pDigAll=new TObjArray(7); pDigAll->SetOwner(); for(Int_t i=0;i<7;i++) pDigAll->AddAt(new TClonesArray("AliHMPIDDigit")  ,i); pTr->Branch("Digs",&pDigAll,64000,0); 
  TObjArray *pCluAll=new TObjArray(7); pCluAll->SetOwner(); for(Int_t i=0;i<7;i++) pCluAll->AddAt(new TClonesArray("AliHMPIDCluster"),i); pTr->Branch("Clus",&pCluAll,64000,0);

  Int_t iEvt=0;
  while(1){      
    Int_t iSize=DateHeader(&in,      isPrint);  if(iSize==68) continue;  //DATE header 
    if(in.eof()) break;
    HmpidHeader (&in,      isPrint);    //HMPID header 
    HmpidPayload(&in,ddl+1,pDigAll);    //HMPID payload
    HmpidHeader (&in,      isPrint);    //HMPID header 
    HmpidPayload(&in,ddl  ,pDigAll);    //HMPID payload
    
    AliHMPIDReconstructor::Dig2Clu(pDigAll,pCluAll);
    pTr->Fill();
    for(Int_t i=0;i<7;i++){
      ((TClonesArray*)pDigAll->At(i))->Clear(); 
      ((TClonesArray*)pCluAll->At(i))->Clear(); 
    }
    
    if(!(iEvt%200)) Printf("Event %i processed",iEvt);
    iEvt++;
    if(iEvt==iMaxEvt) break;
  }//events loop
  
  pTr->Write();  pOut->Close();
  pDigAll->Delete(); pCluAll->Delete();
  in.close();
  
  Printf("Total %i events processed",iEvt);
  gBenchmark->Show("Cosmic converter");
}//docosmic()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void cosmic()
{
  TChain *pCh=new TChain("cosmic");  pCh->Add("test.root");

  
  TH1F *pDigQ=new TH1F("digQ","Digits QDC",500,0,4100);
  TH1F *pDigO=new TH1F("digO","Digits # per event",500,0,8000);
  TH2F *pDigM=new TH2F("digM","Digits map",500,0,131,500,0,127);
    
  TH1F *pCluQ   =new TH1F("cluQ","Clusters QDC",500,0,12100);
  TH1F *pCluQmax=new TH1F("cluQmax","Clusters Max QDC",500,0,12100); pCluQmax->SetLineColor(kRed);
  TH1F *pCluO=new TH1F("cluO","Clusters # per event",500,0,5000);
  TH2F *pCluM=new TH2F("cluM","Clusters map",500,0,131,500,0,127);
  
  TObjArray *pDigAll=0; pCh->SetBranchAddress("Digs",&pDigAll);
  TObjArray *pCluAll=0; pCh->SetBranchAddress("Clus",&pCluAll);

    
  for(Int_t iEvt=0;iEvt<pCh->GetEntries();iEvt++){
    pCh->GetEntry(iEvt);
    
    TClonesArray *pDigCh=(TClonesArray*)pDigAll->At(0);
    TClonesArray *pCluCh=(TClonesArray*)pCluAll->At(0);
    pDigO->Fill(pDigCh->GetEntriesFast());
    pCluO->Fill(pCluCh->GetEntriesFast());
    
    for(Int_t iDig=0;iDig<pDigCh->GetEntriesFast();iDig++){//digits loop
      AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigCh->UncheckedAt(iDig);
      pDigQ->Fill(pDig->Q());
      pDigM->Fill(pDig->LorsX(),pDig->LorsY());
    }//digits loop
    Float_t qmax=0;    
    for(Int_t iClu=0;iClu<pCluCh->GetEntriesFast();iClu++){//clusters loop
      AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluCh->UncheckedAt(iClu);
      pCluQ->Fill(pClu->Q());
      if(pClu->Q()>qmax) qmax=pClu->Q();
      pCluM->Fill(pClu->X(),pClu->Y());
    }//digits loop
    pCluQmax->Fill(qmax);
  }//entries loop
  
  TCanvas *pC=new TCanvas("comic","cosmic"); pC->Divide(2,3);
  
  pC->cd(1); pDigM->Draw(); pC->cd(2); pCluM->Draw();
  pC->cd(3); gPad->SetLogy(); pDigQ->Draw(); pC->cd(4); gPad->SetLogy(); pCluQ->Draw(); pCluQmax->Draw("same");
  pC->cd(5); gPad->SetLogy(); pDigO->Draw(); pC->cd(6); gPad->SetLogy(); pCluO->Draw();
}//cosmic()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Commit	Line	Data
d3da6dc4	1	#include <TF1.h>
	2	#include <TH2F.h>
	3	#include <TCanvas.h> //Terminate()
	4	#include <TChain.h>
	5	#include <TBenchmark.h>
da08475b	6	#include <TFile.h> //docosmic()
d3da6dc4	7	#include <AliSelector.h> //base class
6e06db1d	8	#include <AliESDEvent.h>
da08475b	9	#include <AliBitPacking.h> //HmpidPayload()
da08475b	10	#include "AliHMPIDDigit.h"
1cdf1f6e	11	#include "AliHMPIDParam.h"
da08475b	12	#include "AliHMPIDCluster.h"
da08475b	13	#include "AliHMPIDReconstructor.h" //docosmic()
d3da6dc4	14
	15	class AliHMPIDSelector : public AliSelector {
	16	public :
bd73b83f	17	AliHMPIDSelector():AliSelector(),fChain(0),fEsd(0),fCkovP(0),fMipXY(0),fDifX(0),fSigP(0) {for(Int_t i=0;i<5;i++) fProb[i]=0;}
d3da6dc4	18	virtual ~AliHMPIDSelector() {delete fEsd;}
	19
	20
	21	virtual Int_t Version () const {return 2;}
	22	virtual void Begin (TTree *) {}
	23	virtual void SlaveBegin (TTree *tree);
	24	virtual void Init (TTree *tree);
	25	virtual Bool_t Notify () {return kTRUE;}
	26	virtual Bool_t Process (Long64_t entry);
	27	virtual void SetOption (const char *option) { fOption = option; }
	28	virtual void SetObject (TObject *obj) { fObject = obj; }
	29	virtual void SetInputList (TList *input) {fInput = input;}
	30	virtual TList *GetOutputList () const { return fOutput; }
	31	virtual void SlaveTerminate ();
	32	virtual void Terminate ();
	33
	34	private:
	35	TTree *fChain ; //!pointer to the analyzed TTree or TChain
6e06db1d	36	AliESDEvent *fEsd ; //!
d3da6dc4	37
bd73b83f	38	TH2F fCkovP,fMipXY; //!
	39	TH1F *fDifX; //!
	40	TH2F *fSigP; //!
d3da6dc4	41	TH1F *fProb[5]; //!
	42
	43	ClassDef(AliHMPIDSelector,0);
	44	};
	45
	46
	47	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	48	void AliHMPIDSelector::SlaveBegin(TTree *tree)
	49	{
	50	// The SlaveBegin() function is called after the Begin() function. When running with PROOF SlaveBegin() is called on each slave server.
	51	// The tree argument is deprecated (on PROOF 0 is passed).
	52
	53	Init(tree);
	54
	55	TString option = GetOption();
	56
	57	// create histograms on each slave server
bd73b83f	58	fCkovP = new TH2F("CkovP" , "#theta_{c}, [rad];P, [GeV]", 150, 0, 7 ,500, 0, 1);
d3da6dc4	59	fSigP = new TH2F("SigP" ,"#sigma_{#theta_c}" , 150, 0, 7 ,100, 0, 1e20);
d3da6dc4	60	fMipXY = new TH2F("MipXY" ,"mip position" , 260, 0,130 ,252,0,126);
5582ed3a	61	fDifX = new TH1F("DifX" ,"diff" , 100, -2.5, 2.5);
b0f83594	62
d3da6dc4	63	fProb[0] = new TH1F("PidE" ,"PID: e yellow #mu magenta" ,100,0,1); fProb[0]->SetLineColor(kYellow);
	64	fProb[1] = new TH1F("PidMu","pid of #mu" ,100,0,1); fProb[1]->SetLineColor(kMagenta);
	65	fProb[2] = new TH1F("PidPi","PID: #pi red K green p blue",100,0,1); fProb[2]->SetLineColor(kRed);
	66	fProb[3] = new TH1F("PidK" ,"pid of K" ,100,0,1); fProb[3]->SetLineColor(kGreen);
	67	fProb[4] = new TH1F("PidP" ,"pid of p" ,100,0,1); fProb[4]->SetLineColor(kBlue);
	68	}
	69	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	70	void AliHMPIDSelector::Init(TTree *pTr)
	71	{
ca579b7d	72	// Called when the selector needs to initialize a new tree of chain. Typically here the branch addresses of the tree is set
ca579b7d	73	// Init() will be called many times when running with PROOF.
d3da6dc4	74
ca579b7d	75	if ( !pTr ) return ;
	76	fChain = pTr ;
	77	fChain->SetBranchAddress("ESD", &fEsd) ;
5582ed3a	78	fChain->SetBranchStatus("*", 1);
ca579b7d	79	fChain->SetBranchStatus("fTracks.*", 1);
d3da6dc4	80	}
	81	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	82	Bool_t AliHMPIDSelector::Process(Long64_t entry)
	83	{
1cdf1f6e	84	AliHMPIDParam *pParam=AliHMPIDParam::Instance();
d3da6dc4	85	fChain->GetTree()->GetEntry(entry);
5582ed3a	86
d3da6dc4	87	for(Int_t iTrk=0;iTrk<fEsd->GetNumberOfTracks();iTrk++){
5582ed3a	88	AliESDtrack *pTrk=fEsd->GetTrack(iTrk);
d3da6dc4	89
b0f83594	90	// if(pTrk->GetHMPIDsignal()==-1) continue;
d3da6dc4	91
5582ed3a	92	fCkovP->Fill(pTrk->GetP(),pTrk->GetHMPIDsignal()) ;
5582ed3a	93	fSigP ->Fill(pTrk->GetP(),TMath::Sqrt(pTrk->GetHMPIDchi2()));
d3da6dc4	94
5582ed3a	95	Float_t xClu,yClu; Int_t q,np;
	96	Float_t xPc,yPc;
	97	pTrk->GetHMPIDmip(xClu,yClu,q,np);
	98	fMipXY->Fill(xClu,yClu); //mip info
	99	Float_t xRad,yRad,th,ph;
	100	pTrk->GetHMPIDtrk(xRad,yRad,th,ph);
	101	Int_t iCh=pTrk->GetHMPIDcluIdx();iCh/=1000000;
	102	Double_t p1[3],n1[3]; pParam->Norm(iCh,n1); pParam->Lors2Mars(iCh,0,0,p1,AliHMPIDParam::kPc); //point & norm for PC
	103	if(pTrk->Intersect(p1,n1,fEsd->GetMagneticField())==kFALSE) continue; //try to intersect track with PC
	104	pParam->Mars2Lors (iCh,p1,xPc,yPc); //TRKxPC position
	105	fDifX->Fill(xPc-xClu); //track info
d3da6dc4	106
5582ed3a	107	Double_t pid[5]; pTrk->GetHMPIDpid(pid); for(Int_t i =0;i<5;i++) fProb[i]->Fill(pid[i]);
d3da6dc4	108	}//tracks loop
	109
	110	return kTRUE;
	111	}//Process()
	112	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	113	void AliHMPIDSelector::SlaveTerminate()
	114	{
	115	// The SlaveTerminate() function is called after all entries or objects
	116	// have been processed. When running with PROOF SlaveTerminate() is called
	117	// on each slave server.
	118
	119	// Add the histograms to the output on each slave server
	120
	121	fOutput->Add(fCkovP);
	122	fOutput->Add(fSigP);
	123	fOutput->Add(fMipXY);
bd73b83f	124	fOutput->Add(fDifX);
d3da6dc4	125
	126	for(Int_t i=0;i<5;i++) fOutput->Add(fProb[i]);
	127	}//SlaveTerminate()
	128	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	129	void AliHMPIDSelector::Terminate()
	130	{
ca579b7d	131	// The last function to be called. It always runs on the client, it can be used to present
ca579b7d	132	// the results graphically or save the results to file.
d3da6dc4	133	fCkovP = dynamic_cast<TH2F*>(fOutput->FindObject("CkovP")) ;
	134	fSigP = dynamic_cast<TH2F*>(fOutput->FindObject("SigP")) ;
	135	fMipXY = dynamic_cast<TH2F*>(fOutput->FindObject("MipXY")) ;
bd73b83f	136	fDifX = dynamic_cast<TH1F*>(fOutput->FindObject("DifX")) ;
d3da6dc4	137
	138	fProb[0] = dynamic_cast<TH1F*>(fOutput->FindObject("PidE")) ;
	139	fProb[1] = dynamic_cast<TH1F*>(fOutput->FindObject("PidMu")) ;
	140	fProb[2] = dynamic_cast<TH1F*>(fOutput->FindObject("PidPi")) ;
	141	fProb[3] = dynamic_cast<TH1F*>(fOutput->FindObject("PidK")) ;
	142	fProb[4] = dynamic_cast<TH1F*>(fOutput->FindObject("PidP")) ;
	143
	144	Float_t n=1.292; //mean freon ref idx
	145	TF1 pPi=new TF1("RiPiTheo","acos(sqrt(xx+[0][0])/(x[1]))",1.2,7); pPi->SetLineWidth(1); pPi->SetParameter(1,n);
	146	AliPID ppp; pPi->SetLineColor(kRed); pPi->SetParameter(0,AliPID::ParticleMass(AliPID::kPion)); //mass
	147	TF1 pK=(TF1)pPi->Clone(); pK ->SetLineColor(kGreen); pK ->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon));
	148	TF1 pP=(TF1)pPi->Clone(); pP ->SetLineColor(kBlue); pP ->SetParameter(0,AliPID::ParticleMass(AliPID::kProton));
	149
	150	TCanvas *pC=new TCanvas("c1","ESD QA");pC->SetFillColor(10); pC->SetHighLightColor(10); pC->Divide(3,2);
bd73b83f	151	pC->cd(1); fCkovP->Draw(); pPi->Draw("same"); pK->Draw("same"); pP->Draw("same"); pC->cd(2); fMipXY->Draw(); pC->cd(3); fProb[0]->Draw(); fProb[1]->Draw("same");
bd73b83f	152	pC->cd(4); fSigP ->Draw(); pC->cd(5); fDifX->Draw(); pC->cd(6); fProb[2]->Draw(); fProb[3]->Draw("same"); fProb[4]->Draw("same");
ca579b7d	153	}//Terminate()
d3da6dc4	154	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d3da6dc4	155	void loc()
d3da6dc4	156	{
ca579b7d	157	TChain* pChain =new TChain("esdTree"); pChain->Add("AliESDs.root");
d3da6dc4	158
	159	pChain->Process("AliHMPIDSelector.C+");
	160	}
	161	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
da08475b	162	Int_t DateHeader(ifstream *pFile,Bool_t isPrint=0)
	163	{
	164	Int_t iSize=-1;
	165	pFile->read((char*)&iSize,4);
	166	if(!isPrint)
	167	pFile->seekg(16*4,ios::cur);
	168	else{
	169	Int_t w32=-1;
	170	Printf("");
	171	Printf("Event size %i bytes",iSize); //1 common DATE header 17 words
	172	pFile->read((char*)&w32,4); Printf("Event magic 0x%x" ,w32); //2
	173	pFile->read((char*)&w32,4); Printf("Event head size %i bytes",w32); //3
	174	pFile->read((char*)&w32,4); Printf("Event version 0x%x" ,w32); //4
	175	pFile->read((char*)&w32,4); Printf("Event type %i (%s)" ,w32,(w32==7)? "physics":"SOR"); //5
	176	pFile->read((char*)&w32,4); Printf("Run number %i" ,w32); //6
	177
	178	pFile->read((char*)&w32,4); Printf("Event ID 1 %i" ,w32); //7
	179	pFile->read((char*)&w32,4); Printf("Event ID 2 %i" ,w32); //8
	180
	181	pFile->read((char*)&w32,4); Printf("Trigger pattern 1 %i" ,w32); //9
	182	pFile->read((char*)&w32,4); Printf("Trigger pattern 2 %i" ,w32); //10
	183
	184	pFile->read((char*)&w32,4); Printf("Detector pattern %i" ,w32); //11
	185
	186	pFile->read((char*)&w32,4); Printf("Type attribute 1 %i" ,w32); //12
	187	pFile->read((char*)&w32,4); Printf("Type attribute 2 %i" ,w32); //13
	188	pFile->read((char*)&w32,4); Printf("Type attribute 3 %i" ,w32); //14
	189
	190	pFile->read((char*)&w32,4); Printf("LDC ID %i" ,w32); //15
	191	pFile->read((char*)&w32,4); Printf("GDC ID %i" ,w32); //16
	192	pFile->read((char*)&w32,4); TDatime time(w32); time.Print(); //17
	193
	194	Printf("");
	195	}
	196	return iSize;
	197	}
	198	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	199	void HmpidHeader(ifstream *pFile,Bool_t isPrint=kFALSE)
	200	{
	201	// Prints hmpid trailer and returns number of words for this trailer
	202	if(!isPrint) {pFile->seekg(15*4,ios::cur);return;}
	203	Int_t w32=-1;
	204	Printf("\nHMPID Header:");//private HEADER is 15 words
	205	for(Int_t i=1;i<=11;i++) { pFile->read((char*)&w32,4); Printf("Word #%2i=%12i meaningless",i,w32);}
	206	pFile->read((char*)&w32,4); Printf("Word #12=%12i event counter",w32);
	207	for(Int_t i=13;i<=15;i++){ pFile->read((char*)&w32,4); Printf("Word #%2i=%12i empty",i,w32);}
	208	}
	209	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3c6274c1	210	void HmpidPayload(ifstream pFile,Int_t iDdl,TObjArray pDigAll)
da08475b	211	{
da08475b	212	// payload is 8 sequences with structure: WC36A8 then WC number of w32
3c6274c1	213
	214	TClonesArray pDig1=(TClonesArray)pDigAll->At(iDdl/2); //get list of digits for requested chamber
	215
da08475b	216	UInt_t w32=0;
3c6274c1	217	Int_t iDigCnt=pDig1->GetEntriesFast();
da08475b	218	for(Int_t row=1;row<=8;row++){
	219	pFile->read((char*)&w32,4); Int_t wc=AliBitPacking::UnpackWord(w32,16,31); Int_t ma=AliBitPacking::UnpackWord(w32, 0,15);
	220	if(ma!=0x36a8) Printf("ERROR ERROR ERROR WRONG Marker=0x%x ERROR ERROR ERROR",ma);
	221	for(Int_t i=1;i<=wc;i++){//words loop
	222	pFile->read((char*)&w32,4);
	223	if(w32&0x08000000) continue; //it's DILOGIC CW
	224	AliHMPIDDigit *pDig=new AliHMPIDDigit;
1cdf1f6e	225	pDig->Raw(w32,iDdl);
3c6274c1	226	new ((pDig1)[iDigCnt++]) AliHMPIDDigit(pDig);
da08475b	227	}//words loop
	228	}//rows loop
	229	}//HmpidPayload()
	230	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3c6274c1	231	void docosmic(const char* name,Int_t iMaxEvt=9999999)
da08475b	232	{
	233	gBenchmark->Start("Cosmic converter");
	234	ifstream in(name);
	235
	236
	237	Bool_t isPrint=kFALSE;
	238
	239	TString rooName=name; rooName.Replace(1+rooName.First('.'),4,"root");
	240	Int_t ddl=0;
	241
	242	TFile *pOut=new TFile(rooName,"recreate");
	243	TTree *pTr=new TTree("cosmic","some for time being");
	244
3c6274c1	245
	246	TObjArray *pDigAll=new TObjArray(7); pDigAll->SetOwner(); for(Int_t i=0;i<7;i++) pDigAll->AddAt(new TClonesArray("AliHMPIDDigit") ,i); pTr->Branch("Digs",&pDigAll,64000,0);
	247	TObjArray *pCluAll=new TObjArray(7); pCluAll->SetOwner(); for(Int_t i=0;i<7;i++) pCluAll->AddAt(new TClonesArray("AliHMPIDCluster"),i); pTr->Branch("Clus",&pCluAll,64000,0);
da08475b	248
3c6274c1	249	Int_t iEvt=0;
da08475b	250	while(1){
	251	Int_t iSize=DateHeader(&in, isPrint); if(iSize==68) continue; //DATE header
	252	if(in.eof()) break;
	253	HmpidHeader (&in, isPrint); //HMPID header
3c6274c1	254	HmpidPayload(&in,ddl+1,pDigAll); //HMPID payload
da08475b	255	HmpidHeader (&in, isPrint); //HMPID header
3c6274c1	256	HmpidPayload(&in,ddl ,pDigAll); //HMPID payload
da08475b	257
3c6274c1	258	AliHMPIDReconstructor::Dig2Clu(pDigAll,pCluAll);
da08475b	259	pTr->Fill();
3c6274c1	260	for(Int_t i=0;i<7;i++){
	261	((TClonesArray*)pDigAll->At(i))->Clear();
	262	((TClonesArray*)pCluAll->At(i))->Clear();
	263	}
	264
	265	if(!(iEvt%200)) Printf("Event %i processed",iEvt);
	266	iEvt++;
	267	if(iEvt==iMaxEvt) break;
	268	}//events loop
da08475b	269
3c6274c1	270	pTr->Write(); pOut->Close();
3c6274c1	271	pDigAll->Delete(); pCluAll->Delete();
da08475b	272	in.close();
3c6274c1	273
3c6274c1	274	Printf("Total %i events processed",iEvt);
da08475b	275	gBenchmark->Show("Cosmic converter");
	276	}//docosmic()
	277	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	278	void cosmic()
	279	{
3c6274c1	280	TChain *pCh=new TChain("cosmic"); pCh->Add("test.root");
da08475b	281
	282
	283	TH1F *pDigQ=new TH1F("digQ","Digits QDC",500,0,4100);
	284	TH1F *pDigO=new TH1F("digO","Digits # per event",500,0,8000);
	285	TH2F *pDigM=new TH2F("digM","Digits map",500,0,131,500,0,127);
	286
3c6274c1	287	TH1F *pCluQ =new TH1F("cluQ","Clusters QDC",500,0,12100);
3c6274c1	288	TH1F *pCluQmax=new TH1F("cluQmax","Clusters Max QDC",500,0,12100); pCluQmax->SetLineColor(kRed);
da08475b	289	TH1F *pCluO=new TH1F("cluO","Clusters # per event",500,0,5000);
	290	TH2F *pCluM=new TH2F("cluM","Clusters map",500,0,131,500,0,127);
	291
3c6274c1	292	TObjArray *pDigAll=0; pCh->SetBranchAddress("Digs",&pDigAll);
	293	TObjArray *pCluAll=0; pCh->SetBranchAddress("Clus",&pCluAll);
	294
	295
da08475b	296	for(Int_t iEvt=0;iEvt<pCh->GetEntries();iEvt++){
	297	pCh->GetEntry(iEvt);
	298
3c6274c1	299	TClonesArray pDigCh=(TClonesArray)pDigAll->At(0);
	300	TClonesArray pCluCh=(TClonesArray)pCluAll->At(0);
	301	pDigO->Fill(pDigCh->GetEntriesFast());
	302	pCluO->Fill(pCluCh->GetEntriesFast());
da08475b	303
3c6274c1	304	for(Int_t iDig=0;iDig<pDigCh->GetEntriesFast();iDig++){//digits loop
3c6274c1	305	AliHMPIDDigit pDig=(AliHMPIDDigit)pDigCh->UncheckedAt(iDig);
da08475b	306	pDigQ->Fill(pDig->Q());
	307	pDigM->Fill(pDig->LorsX(),pDig->LorsY());
	308	}//digits loop
1d4857c5	309	Float_t qmax=0;
3c6274c1	310	for(Int_t iClu=0;iClu<pCluCh->GetEntriesFast();iClu++){//clusters loop
3c6274c1	311	AliHMPIDCluster pClu=(AliHMPIDCluster)pCluCh->UncheckedAt(iClu);
da08475b	312	pCluQ->Fill(pClu->Q());
3c6274c1	313	if(pClu->Q()>qmax) qmax=pClu->Q();
da08475b	314	pCluM->Fill(pClu->X(),pClu->Y());
da08475b	315	}//digits loop
3c6274c1	316	pCluQmax->Fill(qmax);
da08475b	317	}//entries loop
	318
	319	TCanvas *pC=new TCanvas("comic","cosmic"); pC->Divide(2,3);
	320
	321	pC->cd(1); pDigM->Draw(); pC->cd(2); pCluM->Draw();
3c6274c1	322	pC->cd(3); gPad->SetLogy(); pDigQ->Draw(); pC->cd(4); gPad->SetLogy(); pCluQ->Draw(); pCluQmax->Draw("same");
3c6274c1	323	pC->cd(5); gPad->SetLogy(); pDigO->Draw(); pC->cd(6); gPad->SetLogy(); pCluO->Draw();
da08475b	324	}//cosmic()
da08475b	325	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++