RICH/AliRICHReconstructor.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include "AliRICHReconstructor.h" //class header
  17 #include "AliRICH.h"              //Reconstruct(...)
  18 #include <AliRawReader.h>         //Reconstruct(...)
  19 #include <AliRun.h>               //ConvertDigits uses gAlice
  20 #include <AliESD.h>            //RichAna()
  21 #include <AliStack.h>          //RichAna()
  22 #include <TFile.h>             //RichAna()
  23 #include <TMinuit.h>              //Dig2Clu()
  24 #include <TParticle.h>            //TParticle()
  25 ClassImp(AliRICHReconstructor)
  26
  27 void AliRICHReconstructor::Reconstruct(AliRunLoader *pAL,AliRawReader* pRR)const
  28 {
  29 //Invoked  by AliReconstruction to convert raw digits from DDL files to clusters
  30 //Arguments: pAL - ALICE run loader pointer
  31 //           pRR - ALICE raw reader pointer
  32 //  Returns: none
  33   AliLoader *pRL=pAL->GetDetectorLoader("RICH");  AliRICH *pRich=(AliRICH*)pAL->GetAliRun()->GetDetector("RICH");//get pointers for RICH and RICH loader
  34
  35   AliRICHDigit dig; //tmp digit, raw digit will be converted to it
  36   TClonesArray *pDigList=new TClonesArray("AliRICHDigit"); Int_t iDigCnt=0; //tmp list of digits for single chamber only
  37
  38   Int_t iEvtN=0;
  39   while(pRR->NextEvent()){//events loop
  40     pAL->GetEvent(iEvtN++);
  41     pRL->MakeTree("R");  pRich->MakeBranch("R");
  42
  43     for(Int_t iChN=1;iChN<=7;iChN++){//chambers loop
  44       pRR->Select(AliRICHDigit::kRichRawId,2*iChN-2,2*iChN-1);//select only DDL files for the current chamber
  45       UInt_t w32=0;
  46       while(pRR->ReadNextInt(w32)){//raw records loop (in selected DDL files)
  47         UInt_t ddl=pRR->GetDDLID(); //returns 0,1,2 ... 13
  48         dig.Raw2Dig(ddl,w32);
  49         AliDebug(1,Form("Ch=%i DDL=%i raw=0x%x digit=(%3i,%3i,%3i,%3i) Q=%5.2f",iChN,ddl,w32,dig.Chamber(),dig.Sector(),dig.PadX(),dig.PadY(),dig.Qdc()));
  50         new((*pDigList)[iDigCnt++]) AliRICHDigit(dig); //add this digit to the tmp list
  51       }//raw records loop
  52       if(iDigCnt) Dig2Clu(pDigList,pRich->Clusters(iChN));//cluster finder for the current chamber if any digits present
  53       pRR->Reset();
  54       pDigList->Delete();  iDigCnt=0;//clean up list of digits for the current chamber
  55     }//chambers loop
  56     pRL->TreeR()->Fill();            //fill tree for current event
  57     pRL->WriteRecPoints("OVERWRITE");//write out clusters for current event
  58     pRich->ClustersReset();
  59   }//events loop
  60   pRL->UnloadRecPoints();
  61 }//Reconstruct raw data
  62 //__________________________________________________________________________________________________
  63 void AliRICHReconstructor::Dig2Clu(TClonesArray *pDigList,TClonesArray *pCluList)const
  64 {
  65 //Finds all clusters for a given digits list provided not empty. Currently digits list provided is a list of all digits for a single chamber.
  66 //Puts all found clusters in the given clusters list.
  67 //Arguments: pDigList - list of digits provided not empty
  68 //  Returns: none
  69   TMatrixF digMap(1,AliRICHParam::NpadsX(),1,AliRICHParam::NpadsY());  digMap=(Float_t)-1; //digit map for one chamber reseted to -1
  70   for(Int_t iDigN=0 ; iDigN < pDigList->GetEntriesFast() ; iDigN++){ //digits loop to fill digits map
  71     AliRICHDigit *pDig= (AliRICHDigit*)pDigList->At(iDigN);
  72     digMap( pDig->PadX(), pDig->PadY() )=iDigN;                     //(padx,pady) cell takes digit number
  73   }
  74
  75   AliRICHCluster clu; Int_t iCluCnt=0; //tmp cluster and cluster counter
  76
  77   for(Int_t iDigN=0;iDigN<pDigList->GetEntriesFast();iDigN++){//digits loop
  78     AliRICHDigit *pDig=(AliRICHDigit*)pDigList->At(iDigN);
  79     if(!(pDig=UseDig(pDig->PadX(),pDig->PadY(),pDigList,&digMap))) continue;  //this digit is already taken in FormCluster(), go after next digit
  80     FormCluster(&clu,pDig,pDigList,&digMap);                            //form cluster starting from this digit
  81     TMinuit *pMinuit=clu.Solve();                              //solve this cluster
  82
  83     if(pMinuit){//means cluster is solved into local maxima number of clusters, so add all of them in loop
  84       Double_t x=-1,y=-1,q=-1;TString str; Double_t b1,b2,b3; Int_t iErrFlg;//tmp to withdraw resulting parameters
  85       for(Int_t i=0;i<clu.Nlocmax();i++){//take resulting parameters
  86         pMinuit->mnpout(3*i  ,str, x, b1, b2, b3, iErrFlg);
  87         pMinuit->mnpout(3*i+1,str, y, b1, b2, b3, iErrFlg);
  88         pMinuit->mnpout(3*i+2,str, q, b1, b2, b3, iErrFlg);
  89         clu.Set(x,y,(Int_t)q);
  90         new((*pCluList)[iCluCnt++]) AliRICHCluster(clu);
  91       }
  92       delete pMinuit;
  93     }else//means cluster is solved as simple center of gravity cluster, add it
  94       new((*pCluList)[iCluCnt++]) AliRICHCluster(clu);
  95     clu.Reset();//make current cluster empty
  96   }//digits loop
  97 }//Dig2Clu()
  98 //__________________________________________________________________________________________________
  99 void  AliRICHReconstructor::FormCluster(AliRICHCluster *pClu,AliRICHDigit *pDig,TClonesArray *pDigList,TMatrixF *pDigMap)const
 100 {
 101 //Forms the initial cluster as a sum of all adjascent digits. Starts from the given digit
 102 //then calls itself recursevly  for all neighbours.
 103 //Arguments: pClu - pointer to cluster being formed
 104 //  Returns: none
 105   pClu->AddDigit(pDig);//take this digit in cluster and mark it as taken
 106
 107   Int_t x[4],y[4];
 108
 109   Int_t iNnei=AliRICHParam::PadNeighbours(pDig->PadX(),pDig->PadY(),x,y);//returns in x,y all possible neighbours of the given one
 110   for (Int_t i=0;i<iNnei;i++)
 111     if((pDig=UseDig(x[i],y[i],pDigList,pDigMap))) FormCluster(pClu,pDig,pDigList,pDigMap);
 112 }//FormCluster()
 113 //__________________________________________________________________________________________________
 114 void AliRICHReconstructor::CheckPR()
 115 {
 116 //Pattern recognition with stack particles
 117   TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition");
 118   TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex:Chamber:Particle");
 119 //  printf("\n\n");
 120 //  printf("Pattern Recognition done for event %5i",0);
 121   AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
 122   AliMagF * magf = gAlice->Field();
 123   AliTracker::SetFieldMap(magf,kTRUE);
 124   for(Int_t iEvtN=0;iEvtN<pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
 125     pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
 126     AliRICHTracker *tr = new AliRICHTracker();
 127     tr->RecWithStack(hn);
 128     AliInfoClass(Form("Pattern Recognition done for event %i \b",iEvtN));
 129 //    printf("\b\b\b\b\b%5i",iEvtN+1);
 130   }
 131   printf("\n\n");
 132   pFile->Write();pFile->Close();
 133 }
 134 //__________________________________________________________________________________________________
 135 void AliRICHReconstructor::RichAna(Int_t iNevMin,Int_t iNevMax,Bool_t askPatRec)
 136 {
 137   TFile *pFile=TFile::Open("AliESDs.root","read");
 138   if(!pFile || !pFile->IsOpen()) {AliInfoClass("ESD file not open.");return;}      //open AliESDs.root
 139   TTree *pTree = (TTree*) pFile->Get("esdTree");
 140   if(!pTree){AliInfoClass("ESD not found.");return;}                               //get ESD tree
 141   AliInfoClass("ESD found. Go ahead!");
 142
 143   AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
 144
 145   AliMagF * magf = gAlice->Field();
 146   AliTracker::SetFieldMap(magf,kTRUE);
 147   pRich->GetLoader()->GetRunLoader()->LoadHeader();
 148   pRich->GetLoader()->GetRunLoader()->LoadKinematics();
 149   TString var1 = "Pmod:Charge:TrackTheta:TrackPhi:MinX:MinY:ThetaCerenkov:NPhotons:";
 150   TString var2 = "ChargeMIP:Chamber:TOF:LengthTOF:prob1:prob2:prob3:";
 151   TString var3 = "ErrPar1:ErrPar2:ErrPar3:Th1:Th2:Th3:nPhotBKG:pdgCode";
 152   TString varList = var1+var2+var3;
 153
 154   Double_t dx,dy;
 155   Double_t hnvec[30];
 156
 157 //  TFile *pFileRA = new TFile("$(HOME)/RichAna.root","RECREATE","RICH Pattern Recognition");
 158   TFile *pFileRA = new TFile("./RichAna.root","RECREATE","RICH Pattern Recognition");
 159   TNtupleD *hn = new TNtupleD("hn","ntuple",varList);
 160   if(iNevMin<0) iNevMin=0;
 161   if(iNevMin>iNevMax) {iNevMin=0;iNevMax=0;}
 162   if(iNevMax==0) iNevMax=999999;
 163   if(pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents()<iNevMax) iNevMax = pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();
 164   AliESD *pESD=new AliESD;  pTree->SetBranchAddress("ESD", &pESD);
 165   for(Int_t iEvtN=iNevMin;iEvtN<iNevMax;iEvtN++) {
 166     pTree->GetEvent(iEvtN);
 167     pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
 168     pRich->GetLoader()->LoadRecPoints();
 169     pRich->GetLoader()->TreeR()->GetEntry(0);
 170     AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
 171 //Pattern recognition started
 172     if(pESD->GetNumberOfTracks()) {
 173       Int_t iNtracks=pESD->GetNumberOfTracks();
 174       AliInfoClass(Form("Start with %i tracks",iNtracks));
 175       for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
 176         if(iTrackN%100==0)AliInfoClass(Form("Track %i to be processed",iTrackN));
 177         AliRICHTracker *pTrRich = new AliRICHTracker();
 178         if(askPatRec==kTRUE) pTrRich->RecWithESD(pESD,pRich,iTrackN);
 179         AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
 180
 181         Int_t lab=TMath::Abs(pTrack->GetLabel());
 182         TParticle *pPart=pStack->Particle(lab);
 183         Int_t code=pPart->GetPdgCode();
 184
 185         pTrack->GetRICHdxdy(dx,dy);
 186         hnvec[0]=pTrack->GetP();
 187         hnvec[1]=pTrack->GetSign();
 188 //  cout << " Track momentum " << pTrack->GetP() << " charge " << pTrack->GetSign() << endl;
 189
 190         pTrack->GetRICHthetaPhi(hnvec[2],hnvec[3]);
 191         pTrack->GetRICHdxdy(hnvec[4],hnvec[5]);
 192         hnvec[6]=pTrack->GetRICHsignal();
 193         hnvec[7]=pTrack->GetRICHnclusters();
 194         hnvec[9]=pTrack->GetRICHcluster()/1000000;
 195         hnvec[8]=pTrack->GetRICHcluster()-hnvec[9]*1000000;
 196         hnvec[10]=pTrack->GetTOFsignal();
 197         hnvec[11]=pTrack->GetIntegratedLength();
 198         Double_t prob[5];
 199         pTrack->GetRICHpid(prob);
 200         hnvec[12]=prob[0]+prob[1]+prob[2];
 201         hnvec[13]=prob[3];
 202         hnvec[14]=prob[4];
 203         hnvec[15]=pTrRich->fErrPar[2];
 204         hnvec[16]=pTrRich->fErrPar[3];
 205         hnvec[17]=pTrRich->fErrPar[4];
 206         for(Int_t i=0;i<3;i++) {
 207           Double_t mass = AliRICHParam::fgMass[i+2];
 208           Double_t refIndex=AliRICHParam::RefIdxC6F14(AliRICHParam::MeanCkovEnergy());
 209           Double_t cosThetaTh = TMath::Sqrt(mass*mass+pTrack->GetP()*pTrack->GetP())/(refIndex*pTrack->GetP());
 210           hnvec[18+i]=0;
 211           if(cosThetaTh>=1) continue;
 212           hnvec[18+i]= TMath::ACos(cosThetaTh);
 213         }
 214         if(askPatRec==kTRUE) hnvec[21]=pTrRich->fnPhotBKG; else hnvec[21]=0;
 215         hnvec[22]=code;
 216         hn->Fill(hnvec);
 217       }
 218     }
 219     AliInfoClass(Form("Pattern Recognition done for event %i",iEvtN));
 220   }
 221   pFileRA->Write();pFileRA->Close();// close RichAna.root
 222   delete pESD;  pFile->Close();//close AliESDs.root
 223 }
 224 //__________________________________________________________________________________________________