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