macros/analHits.C

   1 void analHits (const char *filename="galice.root",Int_t evNumber=0, char *opt="Liny"){
   2 /////////////////////////////////////////////////////////////////////////
   3 //   This macro is a small example of a ROOT macro
   4 //   illustrating how to read the output of GALICE
   5 //   and fill some histograms.
   6 //
   7 //     Root > .L analHits.C               //this loads the macro in memory
   8 //     Root > analHits();                 //by default process first event
   9 //     Root > analHits("galice2.root",2); //process third event from
  10 //                                          galice2.root file.
  11 //Begin_Html
  12 /*
  13 <img src="picts/analHits.gif">
  14 */
  15 //End_Html
  16 /////////////////////////////////////////////////////////////////////////
  17   if(gAlice){
  18     delete gAlice;
  19     gAlice=0;
  20   }
  21   else{
  22     // Dynamically link some shared libs
  23     if(gClassTable->GetID("AliRun") < 0) {
  24       gROOT->LoadMacro("loadlibs.C");
  25       loadlibs();
  26     } // end if
  27   }
  28 // Connect the Root Galice file containing Geometry, Kine and Hits
  29     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
  30     if(!file) file = new TFile(filename);
  31
  32 // Get AliRun object from file or create it if not on file
  33     if(!gAlice) {
  34         gAlice = (AliRun*)file->Get("gAlice");
  35         if(gAlice) printf("AliRun object found on file\n");
  36         if(!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  37     } // end if !gAlice
  38
  39 // Set event pointer to this event
  40     Int_t nparticles = gAlice->GetEvent(evNumber);
  41     if (nparticles <= 0){
  42         cout << "No particles found for event " << evNumber;
  43         cout << " in file " << filename << endl;
  44         return;
  45     } // end if nparticles <=0
  46
  47 // Pointer to specific detector hits.
  48     AliFMDhit    *fmdHit;
  49     AliITShit    *itsHit;
  50     AliMUONHit   *muonHit;
  51     AliPHOSHit   *phosHit;
  52     AliPMDhit    *pmdHit;
  53     AliRICHHit   *richHit;
  54     AliSTARThit  *startHit;
  55     AliTOFhit    *tofHit;
  56     AliTPChit    *tpcHit;
  57     AliTPCTrackHits *tpc2Hit;
  58     AliTRDhit    *trdHit;
  59     AliZDCHit    *zdcHit;
  60     AliEMCALHit  *emcalHit;
  61
  62 // Get pointers to ALL Alice detectors and Hits containers
  63     AliFMD    *FMD    = (AliFMD*)    gAlice->GetDetector("FMD");
  64     AliITS    *ITS    = (AliITS*)    gAlice->GetDetector("ITS");
  65     AliMUON   *MUON   = (AliMUON*)   gAlice->GetDetector("MUON");
  66     AliPHOS   *PHOS   = (AliPHOS*)   gAlice->GetDetector("PHOS");
  67     AliPMD    *PMD    = (AliPMD*)    gAlice->GetDetector("PMD");
  68     AliRICH   *RICH   = (AliRICH*)   gAlice->GetDetector("RICH");
  69     AliSTART  *START  = (AliSTART*)  gAlice->GetDetector("START");
  70     AliTOF    *TOF    = (AliTOF*)    gAlice->GetDetector("TOF");
  71     AliTPC    *TPC    = (AliTPC*)    gAlice->GetDetector("TPC");
  72     AliTPC    *TPC    = (AliTPC*)    gAlice->GetDetector("TPC");
  73     AliTRD    *TRD    = (AliTRD*)    gAlice->GetDetector("TRD");
  74     AliZDC    *ZDC    = (AliZDC*)    gAlice->GetDetector("ZDC");
  75     AliEMCAL  *EMCAL  = (AliEMCAL*)  gAlice->GetDetector("EMCAL");
  76
  77 // Get pointer to the particles
  78 //    TClonesArray *Particles = gAlice->Particles();
  79 //    TParticle    *part;
  80
  81     // Create histograms
  82     if(FMD)   TH1F *hFMD   = new TH1F("hFMD"   ,"Hit Radius",100,0.,100.);
  83     if(ITS)   TH1F *hITS   = new TH1F("hITS"   ,"Ionization",100,0.,3.e-3);
  84     if(MUON)  TH1F *hMUON  = new TH1F("hMUON"  ,"Hit Radius",100,0.,500.);
  85     if(PHOS)   TH1F *hPHOS  = new TH1F("hPHOS"  ,"Energy Dep.",100,0.,0.5);
  86     if(PMD)     TH1F *hPMD   = new TH1F("hPMD"   ,"Energy Dep.",100,0.,1.e+5);
  87     if(RICH)  TH1F *hRICH  = new TH1F("hRICH"  ,"Energy loss",100,0.,1.e-5);
  88     if(START) TH1F *hSTART = new TH1F("hSTART" ,"Time of Flight",100,0.,10.);
  89     if(TOF)   TH1F *hTOF   = new TH1F("hTOF"   ,"Time of Flight",100,0.,1.e-5);
  90     if(TPC)   TH1F *hTPC   = new TH1F("hTPC"   ,"Charge",100,0.,70.2);
  91     if(TRD)   TH1F *hTRD   = new TH1F("hTRD"   ,"Charge",100,0.,10.);
  92     if(ZDC)   TH1F *hZDC   = new TH1F("hZDC"   ,"Energy",100,0.,5.);
  93     if(EMCAL) TH1F *hEMCAL = new TH1F("hEMCAL" ,"Energy",100,0.,2.);
  94 //    TH1F *hTPAR  = new TH1F("hTPAR" ,"?",6,1,7);
  95     Int_t track,ntracks = gAlice->TreeH()->GetEntries();
  96 // Start loop on tracks in the hits containers
  97     for(track=0; track<ntracks;track++){
  98       //MI change
  99       gAlice->ResetHits();
 100       gAlice->TreeH()->GetEvent(track);
 101       if(FMD){
 102             for(fmdHit=(AliFMDhit*)FMD->FirstHit(-1);fmdHit;
 103                 fmdHit=(AliFMDhit*)FMD->NextHit()){
 104                 hFMD->Fill(TMath::Hypot(fmdHit->X(),fmdHit->Y()));
 105             } // end for fmdHit
 106         } // end if FMD
 107         if(ITS){
 108             for(itsHit=(AliITShit*)ITS->FirstHit(-1);itsHit;
 109                 itsHit=(AliITShit*)ITS->NextHit()){
 110                 if(itsHit->GetIonization()>0.0){//only after a step in the ITS
 111                     hITS->Fill(itsHit->GetIonization());
 112                 } // end if
 113             } // end for itsHit
 114         } // end if ITS
 115         if(MUON){
 116             for(muonHit=(AliMUONHit*)MUON->FirstHit(-1);muonHit;
 117                 muonHit=(AliMUONHit*)MUON->NextHit()){
 118                 hMUON->Fill(TMath::Hypot(muonHit->X(),muonHit->Y()));
 119             } // end for muonHit
 120         } // end if MUON
 121         if(PHOS){
 122             for(phosHit=(AliPHOSHit*)PHOS->FirstHit(-1);phosHit;
 123                 phosHit=(AliPHOSHit*)PHOS->NextHit()){
 124                 hPHOS->Fill(phosHit->GetEnergy());
 125             } // end for phosHit
 126         } // end if PHOS
 127         if(PMD){
 128             for(pmdHit=(AliPMDhit*)PMD->FirstHit(-1);pmdHit;
 129                 pmdHit=(AliPMDhit*)PMD->NextHit()){
 130                 hPMD->Fill(pmdHit->GetEnergy());
 131             } // end for pmdHit
 132         } // end if PMD
 133         if(RICH){
 134             for(richHit=(AliRICHHit*)RICH->FirstHit(-1);richHit;
 135                 richHit=(AliRICHHit*)RICH->NextHit()){
 136                 hRICH->Fill(richHit->Eloss());
 137             } // end for richHit
 138         } // end if RICH
 139         if(START){
 140             for(startHit=(AliSTARThit*)START->FirstHit(-1);startHit;
 141                 startHit=(AliSTARThit*)START->NextHit()){
 142                 hSTART->Fill(startHit->fTime);
 143             } // end for startHit
 144         } // end if START
 145         if(TOF){
 146             for(tofHit=(AliTOFhit*)TOF->FirstHit(-1);tofHit;
 147                 tofHit=(AliTOFhit*)TOF->NextHit()){
 148                 hTOF->Fill(tofHit->GetTof());
 149             } // end for tofHit
 150         } // end if TOF
 151
 152         if(TRD) {
 153             for(trdHit=(AliTRDhit*)TRD->FirstHit(-1);trdHit;
 154                 trdHit=(AliTRDhit*)TRD->NextHit()) {
 155                 hTRD->Fill((Float_t)(trdHit->GetCharge()));
 156             } // end for
 157         } // end if TRD
 158         if(ZDC){
 159             for(zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);zdcHit;
 160                 zdcHit=(AliZDCHit*)ZDC->NextHit()){
 161                 hZDC->Fill(zdcHit->GetEnergy());
 162             } // end for zdcdHit
 163         } // end if ZDC
 164         if(TPC) {
 165           for(tpcHit=(AliTPChit*)TPC->FirstHit(-1);tpcHit;
 166               tpcHit=(AliTPChit*)TPC->NextHit()) {
 167             hTPC->Fill((Float_t)(tpcHit->fQ));
 168           } // end for tpcHit
 169         } // end if TPC
 170         if(EMCAL) {
 171           for(emcalHit=(AliEMCALHit*)EMCAL->FirstHit(-1);emcalHit;
 172               emcalHit=(AliEMCALHit*)EMCAL->NextHit()) {
 173             hEMCAL->Fill((Float_t)(emcalHit->GetEnergy()));
 174           } // end for tpcHit
 175         } // end if TPC
 176
 177     } // end for track
 178
 179 //Create a canvas, set the view range, show histograms
 180     TCanvas *c0 = new TCanvas("c0","Alice Detectors",400,10,600,700);
 181     if(opt=="Logy")c0->SetLogy();
 182     if(FMD){
 183         hFMD->SetFillColor(42);
 184         hFMD->Draw();
 185         c0->Print("analHitsFMD.ps");
 186     } // end if FMD
 187     if(ITS){
 188         hITS->SetFillColor(42);
 189         hITS->Draw();
 190         c0->SaveAs("analHitsITS.ps");
 191     } // end if ITS
 192     if(MUON){
 193         hMUON->SetFillColor(42);
 194         hMUON->Draw();
 195         c0->SaveAs("analHitsMUON.ps");
 196     } // end if MUON
 197     if(PHOS){
 198         hPHOS->SetFillColor(42);
 199         hPHOS->Draw();
 200         c0->SaveAs("analHitsPHOS.ps");
 201     } // end if PHOS
 202     if(PMD){
 203         hPMD->SetFillColor(42);
 204         hPMD->Draw();
 205         c0->SaveAs("analHitsPMD.ps");
 206     } // end if PMD
 207     if(RICH){
 208         hRICH->SetFillColor(42);
 209         hRICH->Draw();
 210         c0->SaveAs("analHitsRICH.ps");
 211     } // end if RICH
 212     if(START){
 213         hSTART->SetFillColor(42);
 214         hSTART->Draw();
 215         c0->SaveAs("analHitsSTART.ps");
 216     } // end if START
 217     if(TOF){
 218         hTOF->SetFillColor(42);
 219         hTOF->Draw();
 220         c0->SaveAs("analHitsTOF.ps");
 221     } // end if TOF
 222     if(TPC){
 223         hTPC->SetFillColor(42);
 224         hTPC->Draw();
 225         c0->SaveAs("analHitsTPC.ps");
 226     } // end if TPC
 227     if(TRD){
 228         hTRD->SetFillColor(42);
 229         hTRD->Draw();
 230         c0->SaveAs("analHitsTRD.ps");
 231     } // end if TRD
 232     if(ZDC){
 233         hZDC->SetFillColor(42);
 234         hZDC->Draw();
 235         c0->SaveAs("analHitsZDC.ps");
 236     } // end if ZDC
 237     if(EMCAL){
 238         hEMCAL->SetFillColor(42);
 239         hEMCAL->Draw();
 240         c0->SaveAs("analHitsEMCAL.ps");
 241     } // end if ZDC
 242 }