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