HMPID/AliHMPIDSelector.C

   1 #include <TF1.h>
   2 #include <TH2F.h>
   3 #include <TCanvas.h>  //Terminate()
   4 #include <TChain.h>
   5 #include <TBenchmark.h>
   6 #include <TFile.h>    //docosmic()
   7 #include <AliSelector.h>      //base class
   8 #include <AliESDEvent.h>
   9 #include <AliBitPacking.h> //HmpidPayload()
  10 #include "AliHMPIDDigit.h"
  11 #include "AliHMPIDParam.h"
  12 #include "AliHMPIDCluster.h"
  13 #include "AliHMPIDReconstructor.h" //docosmic()
  14
  15 class AliHMPIDSelector : public AliSelector {
  16  public :
  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;}
  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
  36   AliESDEvent    *fEsd ;     //!
  37
  38   TH2F           *fCkovP,*fMipXY;                //!
  39   TH1F           *fDifX;                         //!
  40   TH2F           *fSigP;                         //!
  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
  58    fCkovP    = new TH2F("CkovP" , "#theta_{c}, [rad];P, [GeV]", 150,   0,  7  ,500, 0, 1);
  59    fSigP     = new TH2F("SigP"  ,"#sigma_{#theta_c}"          , 150,   0,  7  ,100, 0, 1e20);
  60    fMipXY    = new TH2F("MipXY" ,"mip position"               , 260,   0,130  ,252,0,126);
  61    fDifX     = new TH1F("DifX" ,"diff"                       , 100, -2.5, 2.5);
  62
  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 {
  72 // Called when the selector needs to initialize a new tree of chain. Typically here the branch addresses of the tree is set
  73 // Init() will be called many times when running with PROOF.
  74
  75   if ( !pTr )  return ;
  76   fChain = pTr ;
  77   fChain->SetBranchAddress("ESD", &fEsd) ;
  78   fChain->SetBranchStatus("*", 1);
  79   fChain->SetBranchStatus("fTracks.*", 1);
  80 }
  81 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  82 Bool_t AliHMPIDSelector::Process(Long64_t entry)
  83 {
  84   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
  85   fChain->GetTree()->GetEntry(entry);
  86
  87   for(Int_t iTrk=0;iTrk<fEsd->GetNumberOfTracks();iTrk++){
  88     AliESDtrack *pTrk=fEsd->GetTrack(iTrk);
  89
  90 //     if(pTrk->GetHMPIDsignal()==-1) continue;
  91
  92     fCkovP->Fill(pTrk->GetP(),pTrk->GetHMPIDsignal()) ;
  93     fSigP ->Fill(pTrk->GetP(),TMath::Sqrt(pTrk->GetHMPIDchi2()));
  94
  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
 106
 107     Double_t pid[5];  pTrk->GetHMPIDpid(pid); for(Int_t i =0;i<5;i++) fProb[i]->Fill(pid[i]);
 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);
 124   fOutput->Add(fDifX);
 125
 126   for(Int_t i=0;i<5;i++) fOutput->Add(fProb[i]);
 127 }//SlaveTerminate()
 128 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 129 void AliHMPIDSelector::Terminate()
 130 {
 131 // The last function to be called. It always runs on the client, it can be used to present
 132 // the results graphically or save the results to file.
 133   fCkovP   = dynamic_cast<TH2F*>(fOutput->FindObject("CkovP")) ;
 134   fSigP    = dynamic_cast<TH2F*>(fOutput->FindObject("SigP")) ;
 135   fMipXY   = dynamic_cast<TH2F*>(fOutput->FindObject("MipXY")) ;
 136   fDifX    = dynamic_cast<TH1F*>(fOutput->FindObject("DifX")) ;
 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(x*x+[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);
 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");
 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");
 153 }//Terminate()
 154 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 155 void loc()
 156 {
 157   TChain* pChain =new TChain("esdTree");  pChain->Add("AliESDs.root");
 158
 159   pChain->Process("AliHMPIDSelector.C+");
 160 }
 161 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 210 void HmpidPayload(ifstream *pFile,Int_t iDdl,TObjArray *pDigAll)
 211 {
 212 // payload is 8 sequences with structure: WC36A8 then WC number of w32
 213
 214   TClonesArray *pDig1=(TClonesArray*)pDigAll->At(iDdl/2); //get list of digits for requested chamber
 215
 216   UInt_t w32=0;
 217   Int_t iDigCnt=pDig1->GetEntriesFast();
 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;
 225       pDig->Raw(w32,iDdl);
 226       new ((*pDig1)[iDigCnt++]) AliHMPIDDigit(*pDig);
 227     }//words loop
 228   }//rows loop
 229 }//HmpidPayload()
 230 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 231 void docosmic(const char* name,Int_t iMaxEvt=9999999)
 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
 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);
 248
 249   Int_t iEvt=0;
 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
 254     HmpidPayload(&in,ddl+1,pDigAll);    //HMPID payload
 255     HmpidHeader (&in,      isPrint);    //HMPID header
 256     HmpidPayload(&in,ddl  ,pDigAll);    //HMPID payload
 257
 258     AliHMPIDReconstructor::Dig2Clu(pDigAll,pCluAll);
 259     pTr->Fill();
 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
 269
 270   pTr->Write();  pOut->Close();
 271   pDigAll->Delete(); pCluAll->Delete();
 272   in.close();
 273
 274   Printf("Total %i events processed",iEvt);
 275   gBenchmark->Show("Cosmic converter");
 276 }//docosmic()
 277 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 278 void cosmic()
 279 {
 280   TChain *pCh=new TChain("cosmic");  pCh->Add("test.root");
 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
 287   TH1F *pCluQ   =new TH1F("cluQ","Clusters QDC",500,0,12100);
 288   TH1F *pCluQmax=new TH1F("cluQmax","Clusters Max QDC",500,0,12100); pCluQmax->SetLineColor(kRed);
 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
 292   TObjArray *pDigAll=0; pCh->SetBranchAddress("Digs",&pDigAll);
 293   TObjArray *pCluAll=0; pCh->SetBranchAddress("Clus",&pCluAll);
 294
 295
 296   for(Int_t iEvt=0;iEvt<pCh->GetEntries();iEvt++){
 297     pCh->GetEntry(iEvt);
 298
 299     TClonesArray *pDigCh=(TClonesArray*)pDigAll->At(0);
 300     TClonesArray *pCluCh=(TClonesArray*)pCluAll->At(0);
 301     pDigO->Fill(pDigCh->GetEntriesFast());
 302     pCluO->Fill(pCluCh->GetEntriesFast());
 303
 304     for(Int_t iDig=0;iDig<pDigCh->GetEntriesFast();iDig++){//digits loop
 305       AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigCh->UncheckedAt(iDig);
 306       pDigQ->Fill(pDig->Q());
 307       pDigM->Fill(pDig->LorsX(),pDig->LorsY());
 308     }//digits loop
 309     Float_t qmax=0;
 310     for(Int_t iClu=0;iClu<pCluCh->GetEntriesFast();iClu++){//clusters loop
 311       AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluCh->UncheckedAt(iClu);
 312       pCluQ->Fill(pClu->Q());
 313       if(pClu->Q()>qmax) qmax=pClu->Q();
 314       pCluM->Fill(pClu->X(),pClu->Y());
 315     }//digits loop
 316     pCluQmax->Fill(qmax);
 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();
 322   pC->cd(3); gPad->SetLogy(); pDigQ->Draw(); pC->cd(4); gPad->SetLogy(); pCluQ->Draw(); pCluQmax->Draw("same");
 323   pC->cd(5); gPad->SetLogy(); pDigO->Draw(); pC->cd(6); gPad->SetLogy(); pCluO->Draw();
 324 }//cosmic()
 325 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++