RICH/AliRICH.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 "AliRICH.h"
  17 #include "AliRICHParam.h"
  18 #include "AliRICHHelix.h" //ReadESD
  19 #include <TFile.h>
  20 #include <TObjArray.h>
  21 #include <TParticle.h>
  22 #include <AliStack.h>
  23 #include <AliRun.h>
  24 #include <AliMC.h>       //ctor
  25 #include <AliHeader.h>
  26 #include <AliGenEventHeader.h>
  27 #include <AliGenHijingEventHeader.h>
  28 #include <AliESD.h>
  29 #include <TH1F.h>        //HitQA()
  30 #include <TH2F.h>        //Display()
  31 #include <TBenchmark.h>
  32 #include <AliLog.h>
  33 #include <TLatex.h>      //Display()
  34 #include <TCanvas.h>     //Display()
  35 #include <TGraph.h>      //Display()
  36 #include <TStyle.h>      //Display()
  37 #include <TMarker.h>     //Display()
  38 ClassImp(AliRICH)
  39 //__________________________________________________________________________________________________
  40 // RICH manager class
  41 //BEGIN_HTML
  42 /*
  43   <img src="gif/alirich.gif">
  44 */
  45 //END_HTML
  46 //__________________________________________________________________________________________________
  47 AliRICH::AliRICH():AliDetector(),fSdig(0),fSdigCnt(0),fDig(0),fClu(0),fCounters(0),fNcham(0)
  48 {
  49 //Default ctor should not contain any new operators
  50 //AliDetector ctor deals with Hits and Digits
  51 }//AliRICH::AliRICH()
  52 //__________________________________________________________________________________________________
  53 AliRICH::AliRICH(const char *name, const char *title):AliDetector(name,title),fSdig(0),fSdigCnt(0),fDig(0),fClu(0),fCounters(0),fNcham(0)
  54 {
  55 //Named ctor
  56   AliDebug(1,"Start.");
  57 //AliDetector ctor deals with Hits and Digits (reset them to 0, does not create them)
  58   HitCreate();          gAlice->GetMCApp()->AddHitList(fHits);
  59   fNcham=7;
  60   fCounters.ResizeTo(40); fCounters.Zero();
  61   AliDebug(1,"Stop.");
  62 }//AliRICH::AliRICH(const char *name, const char *title)
  63 //__________________________________________________________________________________________________
  64 AliRICH::~AliRICH()
  65 {
  66 //dtor
  67   AliDebug(1,"Start.");
  68
  69
  70   if(fHits)      delete fHits;
  71   if(fSdig)      delete fSdig;
  72   if(fDigits)    delete fDigits;
  73   if(fDig)      {fDig->Delete();   delete fDig;}
  74   if(fClu)      {fClu->Delete();   delete fClu;}
  75   AliDebug(1,"Stop.");
  76 }//AliRICH::~AliRICH()
  77 //__________________________________________________________________________________________________
  78 void AliRICH::BuildGeometry()
  79 {
  80 //Builds a TNode geometry for event display
  81   AliDebug(1,"Start.");
  82
  83   AliDebug(1,"Stop.");
  84 }//void AliRICH::BuildGeometry()
  85 //__________________________________________________________________________________________________
  86 void AliRICH::MakeBranch(Option_t* option)
  87 {
  88 //Create Tree branches for the RICH.
  89   AliDebug(1,Form("Start with option= %s.",option));
  90
  91   const Int_t kBufferSize = 4000;
  92
  93   const char *cH = strstr(option,"H");
  94   const char *cD = strstr(option,"D");
  95   const char *cR = strstr(option,"R");
  96   const char *cS = strstr(option,"S");
  97
  98   if(cH&&TreeH()){//H
  99     HitCreate();      //branch will be created in AliDetector::MakeBranch
 100   }//H
 101   AliDetector::MakeBranch(option);//this is after cH because we need to guarantee that fHits array is created
 102
 103   if(cS&&fLoader->TreeS()){//S
 104     SDigCreate();   MakeBranchInTree(fLoader->TreeS(),"RICH",&fSdig,kBufferSize,0) ;
 105   }//S
 106
 107   if(cD&&fLoader->TreeD()){//D
 108     DigCreate();
 109     for(Int_t i=0;i<fNcham;i++){
 110       MakeBranchInTree(fLoader->TreeD(),Form("%s%d",GetName(),i+1),&((*fDig)[i]),kBufferSize,0);
 111     }
 112   }//D
 113
 114   if(cR&&fLoader->TreeR()){//R
 115     CluCreate();
 116     for(Int_t i=0;i<fNcham;i++)
 117       MakeBranchInTree(fLoader->TreeR(),Form("%sClusters%d",GetName(),i+1), &((*fClu)[i]), kBufferSize, 0);
 118   }//R
 119   AliDebug(1,"Stop.");
 120 }//void AliRICH::MakeBranch(Option_t* option)
 121 //__________________________________________________________________________________________________
 122 void AliRICH::SetTreeAddress()
 123 {
 124 //Set branch address for the Hits and Digits Tree.
 125   AliDebug(1,"Start.");
 126
 127   TBranch *branch;
 128
 129   if(fLoader->TreeH()){//H
 130     AliDebug(1,"tree H is requested.");
 131     HitCreate();//branch map will be in AliDetector::SetTreeAddress
 132   }//H
 133   AliDetector::SetTreeAddress();//this is after TreeH because we need to guarantee that fHits array is created
 134
 135   if(fLoader->TreeS()){//S
 136     AliDebug(1,"tree S is requested.");
 137     branch=fLoader->TreeS()->GetBranch(GetName());        if(branch){SDigCreate();   branch->SetAddress(&fSdig);}
 138   }//S
 139
 140   if(fLoader->TreeD()){//D
 141     AliDebug(1,"tree D is requested.");
 142     for(int i=0;i<fNcham;i++){   branch=fLoader->TreeD()->GetBranch(Form("%s%d",GetName(),i+1));  if(branch){DigCreate(); branch->SetAddress(&((*fDig)[i]));}}
 143   }//D
 144
 145   if(fLoader->TreeR()){//R
 146     AliDebug(1,"tree R is requested.");
 147     for(int i=0;i<fNcham;i++){   branch=fLoader->TreeR()->GetBranch(Form("%sClusters%d" ,GetName(),i+1));  if(branch){CluCreate(); branch->SetAddress(&((*fClu)[i]));}}
 148   }//R
 149   AliDebug(1,"Stop.");
 150 }//void AliRICH::SetTreeAddress()
 151 //__________________________________________________________________________________________________
 152 // AliRICHHit* AliRICH::Hit(Int_t tid)const
 153 // {
 154 // // Search for the first RICH hit belonging to the given tid
 155 //   GetLoader()->LoadHits();
 156 //   for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
 157 //     GetLoader()->TreeH()->GetEntry(iPrimN);
 158 //     for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){
 159 //       AliRICHHit *pHit=(AliRICHHit*)Hits()->At(iHitN);
 160 //       if(tid==pHit->Track()) {GetLoader()->UnloadHits();return pHit;}
 161 //     }//hits
 162 //   }//prims loop
 163 //   GetLoader()->UnloadHits();
 164 //   return 0;
 165 // }
 166 //__________________________________________________________________________________________________
 167 void AliRICH::HitPrint(Int_t iEvtN)const
 168 {
 169 //Prints a list of RICH hits for a given event. Default is event number 0.
 170   if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
 171   AliInfo(Form("List of RICH hits for event %i",iEvtN));
 172   if(GetLoader()->LoadHits()) return;
 173
 174   Int_t iTotalHits=0;
 175   for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
 176     GetLoader()->TreeH()->GetEntry(iPrimN);
 177     Hits()->Print();
 178     iTotalHits+=Hits()->GetEntries();
 179   }
 180   GetLoader()->UnloadHits();
 181   AliInfo(Form("totally %i hits",iTotalHits));
 182 }
 183 //__________________________________________________________________________________________________
 184 void AliRICH::SDigPrint(Int_t iEvtN)const
 185 {
 186 //prints a list of RICH sdigits  for a given event
 187   if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
 188   Info("PrintSDigits","List of RICH sdigits for event %i",iEvtN);
 189   if(GetLoader()->LoadSDigits()) return;
 190
 191   GetLoader()->TreeS()->GetEntry(0);
 192   fSdig->Print();
 193   GetLoader()->UnloadSDigits();
 194   Printf("totally %i sdigits",fSdig->GetEntries());
 195 }
 196 //__________________________________________________________________________________________________
 197 void AliRICH::DigPrint(Int_t iEvtN)const
 198 {
 199 //prints a list of RICH digits  for a given event
 200   if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
 201   Printf("List of RICH digits for event %i",iEvtN);
 202   if(GetLoader()->LoadDigits()) return;
 203
 204   Int_t iTotalDigits=0;
 205   GetLoader()->TreeD()->GetEntry(0);
 206   if(!fDig) return;
 207   for(Int_t iCham=0;iCham<fNcham;iCham++){
 208     TClonesArray *pDigs=(TClonesArray*)fDig->At(iCham);    iTotalDigits+=pDigs->GetEntries();    pDigs->Print();
 209   }
 210   GetLoader()->UnloadDigits();
 211   Printf("totally %i Digits",iTotalDigits);
 212 }
 213 //__________________________________________________________________________________________________
 214 void AliRICH::OccupancyPrint(Int_t iEvtNreq)const
 215 {
 216 //prints occupancy for each chamber in a given event
 217   Int_t iEvtNmin,iEvtNmax;
 218   if(iEvtNreq==-1){
 219     iEvtNmin=0;
 220     iEvtNmax=gAlice->GetEventsPerRun();
 221   } else {
 222     iEvtNmin=iEvtNreq;iEvtNmax=iEvtNreq+1;
 223   }
 224
 225   if(GetLoader()->GetRunLoader()->LoadHeader()) return;
 226   if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
 227
 228 //  Info("Occupancy","for event %i",iEvtN);
 229   if(GetLoader()->LoadHits()) return;
 230   if(GetLoader()->LoadDigits()) return;
 231
 232   Int_t totPadsPerChamber = AliRICHParam::NpadsX()*AliRICHParam::NpadsY();
 233
 234
 235   for(Int_t iEvtN=iEvtNmin;iEvtN<iEvtNmax;iEvtN++){
 236     Int_t nDigCh[kNchambers]={0,0,0,0,0,0,0};
 237     Int_t iChHits[kNchambers]={0,0,0,0,0,0,0};
 238     Int_t nPrim[kNchambers]={0,0,0,0,0,0,0};
 239     Int_t nSec[kNchambers]={0,0,0,0,0,0,0};
 240     AliInfo(Form("events processed %i",iEvtN));
 241     if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
 242     AliStack *pStack = GetLoader()->GetRunLoader()->Stack();
 243     for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
 244       GetLoader()->TreeH()->GetEntry(iPrimN);
 245       for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){
 246         AliRICHHit *pHit = (AliRICHHit *)Hits()->At(iHitN);
 247         if(pHit->Eloss()>0){
 248           iChHits[pHit->C()-1]++;
 249           if(pStack->Particle(pHit->GetTrack())->Rho()<0.01) nPrim[pHit->C()-1]++;else nSec[pHit->C()-1]++;
 250         }
 251       }
 252     }
 253     GetLoader()->TreeD()->GetEntry(0);
 254     for(Int_t iCh=0;iCh<fNcham;iCh++) {
 255       nDigCh[iCh]= ((TClonesArray*)fDig->At(iCh))->GetEntries();
 256       Double_t occupancy = (Double_t)nDigCh[iCh]/(Double_t)totPadsPerChamber;
 257       Info("Occupancy","for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i",
 258         iCh+1,occupancy*100.,nPrim[iCh],nSec[iCh],iChHits[iCh]);
 259     }
 260   }
 261   GetLoader()->UnloadHits();
 262   GetLoader()->UnloadDigits();
 263   GetLoader()->GetRunLoader()->UnloadHeader();
 264   GetLoader()->GetRunLoader()->UnloadKinematics();
 265 }
 266 //__________________________________________________________________________________________________
 267 void AliRICH::CluPrint(Int_t iEvtN)const
 268 {
 269 //prints a list of RICH clusters  for a given event
 270   Printf("List of RICH clusters for event %i",iEvtN);
 271   GetLoader()->GetRunLoader()->GetEvent(iEvtN);
 272   if(GetLoader()->LoadRecPoints()) return;
 273
 274   Int_t iCluCnt=0;
 275   GetLoader()->TreeR()->GetEntry(0);
 276   for(Int_t iCham=0;iCham<fNcham;iCham++){
 277     TClonesArray *pClus=(TClonesArray*)fClu->At(iCham);    iCluCnt+=pClus->GetEntries();    pClus->Print();
 278   }
 279   GetLoader()->UnloadRecPoints();
 280   Printf("totally %i clusters for event %i",iCluCnt,iEvtN);
 281 }
 282 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 283 void AliRICH::DisplayEvent(Int_t iEvtNmin,Int_t iEvtNmax)const
 284 {
 285 // Display digits, reconstructed tracks intersections and RICH rings if available
 286   TH2F *pH2[8];
 287
 288   GetLoader()->LoadDigits();
 289
 290   TLatex t;  t.SetTextSize(0.1);
 291   TCanvas *pC = new TCanvas("RICHDisplay","RICH Display",0,0,1226,900);  pC->Divide(3,3);  pC->cd(9); t.DrawText(0.2,0.4,"View to IP");
 292   gStyle->SetPalette(1);
 293
 294
 295   for(Int_t iCh=1;iCh<=fNcham;iCh++) {
 296     pH2[iCh] = new TH2F(Form("RichDigH2_%i",iCh),Form("Chamber %i;cm;cm",iCh),165,0,AliRICHParam::PcSizeX(),144,0,AliRICHParam::PcSizeY());
 297     pH2[iCh]->SetMarkerColor(kGreen);
 298     pH2[iCh]->SetMarkerStyle(29);
 299     pH2[iCh]->SetMarkerSize(0.4);
 300     pH2[iCh]->SetStats(kFALSE);
 301     pH2[iCh]->SetMaximum(300);
 302   }
 303
 304   if(iEvtNmax>gAlice->GetEventsPerRun()||iEvtNmax==0) iEvtNmax=gAlice->GetEventsPerRun()-1;
 305
 306   for(Int_t iEvt=iEvtNmin;iEvt<=iEvtNmax;iEvt++) {//events loop
 307     pC->cd(3);  t.DrawText(0.2,0.4,Form("Event %i",iEvt));
 308
 309     GetLoader()->GetRunLoader()->GetEvent(iEvt); //get event
 310     GetLoader()->TreeD()->GetEntry(0);           //get list of digits
 311     for(Int_t iCh=1;iCh<=fNcham;iCh++) {//chambers loop
 312       pH2[iCh]->Reset();
 313       for(Int_t iDig=0;iDig < Digs(iCh)->GetEntries();iDig++) {//digits loop
 314         AliRICHDigit *pDig = (AliRICHDigit*)Digs(iCh)->At(iDig);
 315         TVector2 x2=AliRICHParam::Pad2Loc(pDig->Pad());
 316         pH2[pDig->C()]->Fill(x2.X(),x2.Y(),pDig->Qdc());
 317       }//digits loop
 318       if(iCh==1) pC->cd(9);
 319       if(iCh==2) pC->cd(8);
 320       if(iCh==3) pC->cd(6);
 321       if(iCh==4) pC->cd(5);
 322       if(iCh==5) pC->cd(4);
 323       if(iCh==6) pC->cd(2);
 324       if(iCh==7) pC->cd(1);
 325       pH2[iCh]->Draw("col");
 326       ReadESD(iEvt,iCh);
 327       AliRICHParam::DrawSectors();
 328     }//chambers loop
 329     pC->Update();
 330     pC->Modified();
 331
 332     if(iEvt<iEvtNmax) {gPad->WaitPrimitive();pC->Clear();}
 333   }//events loop
 334 }//ShowEvent()
 335 //__________________________________________________________________________________________________
 336 void AliRICH::Display()const
 337 {
 338 //Provides fast event display
 339 //For RICH only, full display is .x Display.C
 340   Bool_t isHits    =!GetLoader()->LoadHits();
 341   Bool_t isDigits  =!GetLoader()->LoadDigits();
 342   Bool_t isClusters=!GetLoader()->LoadRecPoints();
 343
 344   if(!isHits && !isDigits && !isClusters){Error("Exec","No hits digits and clusters. Nothing to display.");return;}
 345
 346   TCanvas *pCanvas = new TCanvas("Display","RICH Display",0,0,600,600);
 347
 348   TH2F *pHitsH2=0,*pDigitsH2=0,*pClustersH2=0;
 349
 350   if(isHits)     pHitsH2     = new TH2F("pHitsH2"  ,  "Event Display;x,cm;y,cm",165,0,AliRICHParam::PcSizeX(),
 351                                                                                 144,0,AliRICHParam::PcSizeY());
 352   if(pHitsH2)    pHitsH2->SetStats(kFALSE);
 353
 354   if(isDigits)   pDigitsH2   = new TH2F("pDigitsH2"  ,"Event Display",165,0,AliRICHParam::PcSizeX(),
 355                                                                       144,0,AliRICHParam::PcSizeY());
 356   if(isClusters) pClustersH2 = new TH2F("pClustersH2","Event Display",165,0,AliRICHParam::PcSizeX(),
 357                                                                       144,0,AliRICHParam::PcSizeY());
 358
 359   for(Int_t iEvt=0;iEvt<GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvt++){//events Loop
 360     GetLoader()->GetRunLoader()->GetEvent(iEvt);
 361 //display all the staff on chamber by chamber basis
 362     for(Int_t iCh=1;iCh<=fNcham;iCh++){//chambers loop
 363       if(isHits)     pHitsH2    ->Reset();
 364       if(isDigits)   pDigitsH2  ->Reset();
 365       if(isClusters) pClustersH2->Reset();
 366 //deals with hits
 367       for(Int_t i=0;i<GetLoader()->TreeH()->GetEntries();i++){//TreeH loop
 368         GetLoader()->TreeH()->GetEntry(i);
 369         for(Int_t iHit=0;iHit<Hits()->GetEntries();iHit++){//hits loop
 370           AliRICHHit *pHit = (AliRICHHit*)Hits()->At(iHit);
 371           if(pHit->C()==iCh){
 372             TVector2 hitLocX2 = AliRICHParam::Instance()->Mars2Lors(iCh,pHit->OutX3());
 373             pHitsH2->Fill(hitLocX2.X(),hitLocX2.Y(),200);
 374           }//if
 375         }//hits loop
 376       }//TreeH loop
 377       pHitsH2->SetTitle(Form("event %i chamber %2i",iEvt,iCh));
 378       pHitsH2->SetMarkerColor(kRed); pHitsH2->SetMarkerStyle(29); pHitsH2->SetMarkerSize(0.4);
 379       pHitsH2->Draw();
 380       AliRICHParam::DrawSectors();
 381       TLatex l; l.SetNDC(); l.SetTextSize(0.02);
 382       if(!isHits)     {l.SetTextColor(kRed)  ;l.DrawLatex(0.1,0.01,"No Hits"    );}
 383       if(!isDigits)   {l.SetTextColor(kGreen);l.DrawLatex(0.4,0.01,"No DIGITS"  );}
 384       if(!isClusters) {l.SetTextColor(kBlue) ;l.DrawLatex(0.8,0.01,"No CLUSTERS");}
 385       pCanvas->Update();        pCanvas->Modified();       gPad->WaitPrimitive();
 386 //deals with digits
 387       if(isDigits){
 388         GetLoader()->TreeD()->GetEntry(0);
 389         for(Int_t iDig=0;iDig < Digs(iCh)->GetEntries();iDig++){//digits loop
 390           AliRICHDigit *pDig = (AliRICHDigit*)Digs(iCh)->At(iDig);
 391                 TVector2 x2=AliRICHParam::Pad2Loc(pDig->Pad());
 392                 pDigitsH2->Fill(x2.X(),x2.Y(),100);
 393         }//digits loop
 394         pDigitsH2->SetMarkerColor(kGreen); pDigitsH2->SetMarkerStyle(29); pDigitsH2->SetMarkerSize(0.4);
 395         pDigitsH2->Draw("same");
 396         pCanvas->Update();        pCanvas->Modified();       gPad->WaitPrimitive();
 397       }//if(isDigits)
 398 //deals with clusters
 399       if(isClusters){
 400         GetLoader()->TreeR()->GetEntry(0);
 401         for(Int_t iClu=0;iClu<Clus(iCh)->GetEntries();iClu++){//clusters loop
 402           AliRICHCluster *pClu = (AliRICHCluster*)Clus(iCh)->At(iClu);
 403           pClustersH2->Fill(pClu->X(),pClu->Y(),50);
 404         }//clusters loop
 405         pClustersH2->SetMarkerColor(kBlue); pClustersH2->SetMarkerStyle(29);  pClustersH2->SetMarkerSize(0.4);
 406         pClustersH2->Draw("same");
 407         pCanvas->Update();        pCanvas->Modified();       gPad->WaitPrimitive();
 408       }//if(isClusters)
 409     }//chambers loop
 410   }//events Loop
 411
 412   delete pCanvas;
 413   GetLoader()->UnloadHits();
 414   if(isDigits)   GetLoader()->UnloadDigits();
 415   if(isClusters) GetLoader()->UnloadRecPoints();
 416 }//Display()
 417 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 418 void AliRICH::ReadESD(Int_t iEventN, Int_t iChamber)const
 419 {
 420 //
 421   TFile *pFile=TFile::Open("AliESDs.root","read");
 422   if(!pFile || !pFile->IsOpen()) {AliInfo("ESD file not open.");return;}      //open AliESDs.root
 423   TTree *pTree = (TTree*) pFile->Get("esdTree");
 424   if(!pTree){AliInfo("ESD not found.");return;}                               //get ESD tree
 425
 426
 427   AliESD *pESD=new AliESD;  pTree->SetBranchAddress("ESD", &pESD);
 428
 429   pTree->GetEvent(iEventN);
 430
 431   Int_t iNtracks=pESD->GetNumberOfTracks();
 432
 433   for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
 434     AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
 435     Double_t mom[3], pos[3];
 436     pTrack->GetPxPyPz(mom); TVector3 mom3(mom[0],mom[1],mom[2]);
 437     pTrack->GetXYZ(pos); TVector3 pos3(pos[0],pos[1],pos[2]);
 438     AliRICHHelix helix(pos3,mom3,(Int_t)pTrack->GetSign(),-0.1*pESD->GetMagneticField());
 439     Int_t iChamberOnRICH=helix.RichIntersect(AliRICHParam::Instance());
 440     if(iChamberOnRICH==iChamber) {
 441       TMarker *trackImpact = new TMarker(helix.PosPc().X(),helix.PosPc().Y(),kStar);
 442       trackImpact->SetMarkerColor(kRed);
 443       trackImpact->Draw();
 444 //
 445       Int_t iChamberRecon = pTrack->GetRICHcluster()/1000000;
 446       if(iChamberRecon==iChamber) {
 447         Double_t thetaCer = pTrack->GetRICHsignal();
 448         if(thetaCer<0) continue;
 449         TVector3 entrance(helix.PosRad().X(),helix.PosRad().Y(),0);
 450         Float_t thetaTrack,phiTrack;
 451         pTrack->GetRICHthetaPhi(thetaTrack,phiTrack);
 452         TVector3 vectorTrack;
 453         vectorTrack.SetMagThetaPhi(pTrack->GetP(),thetaTrack,phiTrack);
 454         AliInfo(Form("Draw ring started for track %i on chamber %i",iTrackN,iChamber));
 455         AliInfo(Form("ThetaCer %f TrackTheta %f TrackPhi %f Momentum %f",thetaCer,thetaTrack,phiTrack,pTrack->GetP()));
 456         Float_t dx,dy;
 457         pTrack->GetRICHdxdy(dx,dy);
 458         DrawRing(entrance,vectorTrack,thetaCer);
 459       }
 460     }
 461   }
 462   delete pESD;  pFile->Close();//close AliESDs.root
 463 }
 464 //__________________________________________________________________________________________________
 465 void AliRICH::DrawRing(TVector3 entrance,TVector3 vectorTrack,Double_t thetaCer)const
 466 {
 467   Double_t xGraph[100],yGraph[100];
 468   Int_t nPointsToDraw = 0;
 469   for(Int_t i=0;i<100;i++) {
 470     Double_t phiCer = 2*TMath::Pi()*i/100;
 471     TVector3 pos = AliRICHParam::Instance()->ForwardTracing(entrance,vectorTrack,thetaCer,phiCer);
 472     if(pos.X()==-999) continue;
 473     xGraph[nPointsToDraw] = pos.X();yGraph[nPointsToDraw] = pos.Y();nPointsToDraw++;
 474   }
 475 //  AliInfo(Form("Npoints per ring %i",nPointsToDraw));
 476   TGraph *gra = new TGraph(nPointsToDraw,xGraph,yGraph);
 477   gra->Draw("C");
 478 }
 479 //__________________________________________________________________________________________________
 480 void AliRICH::SummaryOfEvent(Int_t iEvtN) const
 481 {
 482 //prints a summary for a given event
 483   AliInfo(Form("Summary of event %i",iEvtN));
 484   GetLoader()->GetRunLoader()->GetEvent(iEvtN);
 485   if(GetLoader()->GetRunLoader()->LoadHeader()) return;
 486   if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
 487   AliStack *pStack=GetLoader()->GetRunLoader()->Stack();
 488
 489   AliGenEventHeader* pGenHeader =  gAlice->GetHeader()->GenEventHeader();
 490   if(pGenHeader->InheritsFrom("AliGenHijingEventHeader")) {
 491     AliInfo(Form(" Hijing event with impact parameter b = %.2f (fm)",((AliGenHijingEventHeader*) pGenHeader)->ImpactParameter()));
 492   }
 493   Int_t nChargedPrimaries=0;
 494   for(Int_t i=0;i<pStack->GetNtrack();i++) {
 495     TParticle *pParticle = pStack->Particle(i);
 496     if(pParticle->IsPrimary()&&pParticle->GetPDG()->Charge()!=0) nChargedPrimaries++;
 497     }
 498   AliInfo(Form("Total number of         primaries %i",pStack->GetNprimary()));
 499   AliInfo(Form("Total number of charged primaries %i",nChargedPrimaries));
 500   AliInfo(Form("Total n. of tracks in stack(+sec) %i",pStack->GetNtrack()));
 501   GetLoader()->GetRunLoader()->UnloadHeader();
 502   GetLoader()->GetRunLoader()->UnloadKinematics();
 503 }
 504 //__________________________________________________________________________________________________
 505 void AliRICH::HitQA(Double_t cut,Double_t cutele,Double_t cutR)
 506 {
 507 // Provides a set of control plots intended primarily for charged particle flux analisys
 508 // Arguments: cut (GeV)    - cut on momentum of any charged particles but electrons,
 509 //            cetele (GeV) - the same for electrons-positrons
 510 //            cutR (cm)    - cut on production vertex radius (cylindrical system)
 511   gBenchmark->Start("HitsAna");
 512
 513   Double_t cutPantiproton    =cut;
 514   Double_t cutPkaonminus     =cut;
 515   Double_t cutPpionminus     =cut;
 516   Double_t cutPmuonminus     =cut;
 517   Double_t cutPpositron      =cutele;
 518
 519   Double_t cutPelectron      =cutele;
 520   Double_t cutPmuonplus      =cut;
 521   Double_t cutPpionplus      =cut;
 522   Double_t cutPkaonplus      =cut;
 523   Double_t cutPproton        =cut;
 524
 525
 526   TH2F *pEleHitRZ    =new TH2F("EleHitRZ"    ,Form("e^{+} e^{-} hit %s;z[cm];R[cm]" ,GetName())     , 400,-300,300 ,400,-500,500);   //R-z plot 0cm<R<550cm -300cm<z<300cm
 527   TH2F *pEleHitRP    =new TH2F("EleHitRP"    ,Form("e^{+} e^{-} hit %s;p[GeV];R[cm]",GetName())     ,1000,-1  ,1   ,400,   0,550);   //R-p plot 0cm<R<550cm -1GeV<p<1GeV
 528   TH1F *pEleAllP     =new TH1F("EleAllP"     ,     "e^{+} e^{-} all;p[GeV]"                         ,1000,-1  ,1                );
 529   TH1F *pEleHitP     =new TH1F("EleHitP"     ,Form("e^{+} e^{-} hit %s;p[GeV]"      ,GetName())     ,1000,-1  ,1                );
 530   TH1F *pMuoHitP     =new TH1F("MuoHitP"     ,Form("#mu^{-} #mu^{+} hit %s;p[GeV]"  ,GetName())     ,1000,-4  ,4                );
 531   TH1F *pPioHitP     =new TH1F("PioHitP"     ,Form("#pi^{-} #pi^{+} hit %s;p[GeV]"  ,GetName())     ,1000,-4  ,4                );
 532   TH1F *pKaoHitP     =new TH1F("KaoHitP"     ,Form("K^{-} K^{+} hit %s;p[GeV]"      ,GetName())     ,1000,-4  ,4                );
 533   TH1F *pProHitP     =new TH1F("ProHitP"     ,Form("p^{-} p^{+} hit %s;p[GeV]"      ,GetName())     ,1000,-4  ,4                );
 534   TH2F *pFlux        =new TH2F("flux"        ,Form("%s flux with Rvertex<%.1fcm"    ,GetName(),cutR),10  ,-5  ,5   , 10,0    ,10); //special text hist
 535   TH2F *pVertex      =new TH2F("vertex"      ,Form("%s 2D vertex of RICH hit;x;y"   ,GetName())     ,120 ,0   ,600 ,120,0    ,600); //special text hist
 536   TH1F *pRho         =new TH1F("rho"         ,Form("%s r of RICH hit"               ,GetName())     ,600 ,0   ,600); //special text hist
 537   pFlux->SetStats(0);
 538   pFlux->GetXaxis()->SetBinLabel(1 ,Form("p^{-}>%.3fGeV/c"   ,cutPantiproton));
 539   pFlux->GetXaxis()->SetBinLabel(2 ,Form("K^{-}>%.3fGeV/c"   ,cutPkaonminus ));
 540   pFlux->GetXaxis()->SetBinLabel(3 ,Form("#pi^{-}>%.3fGeV/c" ,cutPpionminus ));
 541   pFlux->GetXaxis()->SetBinLabel(4 ,Form("#mu^{-}>%.3fGeV/c" ,cutPmuonminus ));
 542   pFlux->GetXaxis()->SetBinLabel(5 ,Form("e^{+}>%.3fGeV/c"   ,cutPpositron  ));
 543
 544   pFlux->GetXaxis()->SetBinLabel(6 ,Form("e^{-}>%.3fGeV/c"   ,cutPelectron  ));
 545   pFlux->GetXaxis()->SetBinLabel(7 ,Form("#mu^{+}>%.3fGeV/c" ,cutPmuonplus  ));
 546   pFlux->GetXaxis()->SetBinLabel(8 ,Form("#pi^{+}>%.3fGeV/c" ,cutPpionplus  ));
 547   pFlux->GetXaxis()->SetBinLabel(9 ,Form("K^{+}>%.3fGeV/c"   ,cutPkaonplus  ));
 548   pFlux->GetXaxis()->SetBinLabel(10,Form("p^{+}>%.3fGeV/c"   ,cutPproton    ));
 549
 550   pFlux->GetYaxis()->SetBinLabel(1,"sum");
 551   pFlux->GetYaxis()->SetBinLabel(2,"ch1");
 552   pFlux->GetYaxis()->SetBinLabel(3,"ch2");
 553   pFlux->GetYaxis()->SetBinLabel(4,"ch3");
 554   pFlux->GetYaxis()->SetBinLabel(5,"ch4");
 555   pFlux->GetYaxis()->SetBinLabel(6,"ch5");
 556   pFlux->GetYaxis()->SetBinLabel(7,"ch6");
 557   pFlux->GetYaxis()->SetBinLabel(8,"ch7");
 558   pFlux->GetYaxis()->SetBinLabel(9,"prim");
 559   pFlux->GetYaxis()->SetBinLabel(10,"tot");
 560
 561 //end of hists definition
 562
 563   Int_t iNevents=fLoader->GetRunLoader()->GetAliRun()->GetEventsPerRun(),iCntPrimParts=0,iCntTotParts=0;
 564 //load all needed trees
 565   fLoader->LoadHits();
 566   fLoader->GetRunLoader()->LoadHeader();
 567   fLoader->GetRunLoader()->LoadKinematics();
 568
 569   for(Int_t iEvtN=0;iEvtN < iNevents;iEvtN++){//events loop
 570     fLoader->GetRunLoader()->GetEvent(iEvtN);
 571     AliInfo(Form(" %i event processes",fLoader->GetRunLoader()->GetEventNumber()));
 572     AliStack *pStack= fLoader->GetRunLoader()->Stack();
 573
 574     for(Int_t iParticleN=0;iParticleN<pStack->GetNtrack();iParticleN++){//stack loop
 575       TParticle *pPart=pStack->Particle(iParticleN);
 576
 577       if(iParticleN%10000==0) AliInfo(Form(" %i particles read",iParticleN));
 578
 579       switch(pPart->GetPdgCode()){
 580         case kProtonBar: pFlux->Fill(-4.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-4.5,8); break;
 581         case kKMinus:    pFlux->Fill(-3.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-3.5,8); break;
 582         case kPiMinus:   pFlux->Fill(-2.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-2.5,8); break;
 583         case kMuonMinus: pFlux->Fill(-1.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-1.5,8); break;
 584         case kPositron:  pFlux->Fill(-0.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-0.5,8); pEleAllP->Fill(-pPart->P()); break;
 585
 586         case kElectron:  pFlux->Fill( 0.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 0.5,8); pEleAllP->Fill( pPart->P()); break;
 587         case kMuonPlus:  pFlux->Fill( 1.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 1.5,8); break;
 588         case kPiPlus:    pFlux->Fill( 2.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 2.5,8); break;
 589         case kKPlus:     pFlux->Fill( 3.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 3.5,8); break;
 590         case kProton:    pFlux->Fill( 4.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 4.5,8); break;
 591       }//switch
 592     }//stack loop
 593 //now hits analiser
 594     for(Int_t iEntryN=0;iEntryN < fLoader->TreeH()->GetEntries();iEntryN++){//TreeH loop
 595       fLoader->TreeH()->GetEntry(iEntryN);                                  //get current entry (prim)
 596       for(Int_t iHitN=0;iHitN < Hits()->GetEntries();iHitN++){//hits loop
 597         AliRICHHit *pHit = (AliRICHHit*)Hits()->At(iHitN);            //get current hit
 598         TParticle  *pPart=pStack->Particle(pHit->GetTrack());      //get stack particle which produced the current hit
 599
 600         if(pPart->GetPDG()->Charge()!=0&&pPart->Rho()>0.1) pVertex->Fill(pPart->Vx(),pPart->Vy()); //safe margin for sec.
 601         if(pPart->GetPDG()->Charge()!=0) pRho->Fill(pPart->Rho()); //safe margin for sec.
 602         if(pPart->R()>cutR) continue;                                   //cut on production radius (cylindrical system)
 603
 604         switch(pPart->GetPdgCode()){
 605           case kProtonBar: if(pPart->P()>cutPantiproton) {pProHitP->Fill(-pPart->P()); pFlux->Fill(-4.5,pHit->C());}break;
 606           case kKMinus   : if(pPart->P()>cutPkaonminus)  {pKaoHitP->Fill(-pPart->P()); pFlux->Fill(-3.5,pHit->C());}break;
 607           case kPiMinus  : if(pPart->P()>cutPpionminus)  {pPioHitP->Fill(-pPart->P()); pFlux->Fill(-2.5,pHit->C());}break;
 608           case kMuonMinus: if(pPart->P()>cutPmuonminus)  {pMuoHitP->Fill(-pPart->P()); pFlux->Fill(-1.5,pHit->C());}break;
 609           case kPositron : if(pPart->P()>cutPpositron)   {pEleHitP->Fill(-pPart->P()); pFlux->Fill(-0.5,pHit->C());
 610                pEleHitRP->Fill(-pPart->P(),pPart->R());  pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
 611
 612           case kElectron : if(pPart->P()>cutPelectron)   {pEleHitP->Fill( pPart->P()); pFlux->Fill( 0.5,pHit->C());
 613                pEleHitRP->Fill( pPart->P(),pPart->R());  pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
 614           case kMuonPlus : if(pPart->P()>cutPmuonplus)   {pMuoHitP->Fill( pPart->P()); pFlux->Fill( 1.5,pHit->C());}break;
 615           case kPiPlus   : if(pPart->P()>cutPpionplus)   {pPioHitP->Fill( pPart->P()); pFlux->Fill( 2.5,pHit->C());}break;
 616           case kKPlus    : if(pPart->P()>cutPkaonplus)   {pKaoHitP->Fill( pPart->P()); pFlux->Fill( 3.5,pHit->C());}break;
 617           case kProton   : if(pPart->P()>cutPproton)     {pProHitP->Fill( pPart->P()); pFlux->Fill( 4.5,pHit->C());}break;
 618         }
 619       }//hits loop
 620     }//TreeH loop
 621     iCntPrimParts +=pStack->GetNprimary();
 622     iCntTotParts  +=pStack->GetNtrack();
 623   }//events loop
 624 //unload all loaded staff
 625   fLoader->UnloadHits();
 626   fLoader->GetRunLoader()->UnloadHeader();
 627   fLoader->GetRunLoader()->UnloadKinematics();
 628 //Calculater some sums
 629   Stat_t sum=0;
 630 //sum row, sum over rows
 631   for(Int_t i=1;i<=pFlux->GetNbinsX();i++){
 632     sum=0; for(Int_t j=2;j<=8;j++)    sum+=pFlux->GetBinContent(i,j);
 633     pFlux->SetBinContent(i,1,sum);
 634   }
 635
 636 //display everything
 637   new TCanvas("canvas1",Form("Events %i Nprims=%i Nparticles=%i",iNevents,iCntPrimParts,iCntTotParts),1000,900); pFlux->Draw("text");  gPad->SetGrid();
 638 //total prims and particles
 639   TLatex latex; latex.SetTextSize(0.02);
 640   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,10);    latex.DrawLatex(5.1,9.5,Form("%.0f",sum));
 641   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,9);     latex.DrawLatex(5.1,8.5,Form("%.0f",sum));
 642   for(Int_t iChN=1;iChN<=kNchambers;iChN++) {
 643     sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,iChN+1);latex.DrawLatex(5.1,iChN+0.5,Form("%.0f",sum));
 644   }
 645   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,1);    latex.DrawLatex(5.1,0.5,Form("%.0f",sum));
 646
 647   new TCanvas("cEleAllP"   ,"e" ,200,100); pEleAllP->Draw();
 648   new TCanvas("cEleHitRP"  ,"e" ,200,100); pEleHitRP->Draw();
 649   new TCanvas("cEleHitRZ"  ,"e" ,200,100); pEleHitRZ->Draw();
 650   new TCanvas("cEleHitP"   ,"e" ,200,100); pEleHitP->Draw();
 651   new TCanvas("cMuoHitP"   ,"mu",200,100); pMuoHitP->Draw();
 652   new TCanvas("cPioHitP"   ,"pi",200,100); pPioHitP->Draw();
 653   new TCanvas("cKaoHitP"   ,"K" ,200,100); pKaoHitP->Draw();
 654   new TCanvas("cProHitP"   ,"p" ,200,100); pProHitP->Draw();
 655   new TCanvas("cVertex"    ,"2d vertex" ,200,100); pVertex->Draw();
 656   new TCanvas("cRho"    ,"Rho of sec" ,200,100); pRho->Draw();
 657
 658   gBenchmark->Show("HitsPlots");
 659 }//HitsPlots()