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     AliHMPIDHit   *richHit;
  54     AliT0hit  *startHit;
  55     AliTOFhit    *tofHit;
  56     AliTPChit    *tpcHit;
  57     AliTPCTrackHits *tpc2Hit;
  58     AliTRDhit    *trdHit;
  59     AliCASTORhit *castorHit;
  60     AliZDCHit    *zdcHit;
  61     AliEMCALHit  *emcalHit;
  62
  63 // Get pointers to ALL Alice detectors and Hits containers
  64     AliFMD    *FMD    = (AliFMD*)    gAlice->GetDetector("FMD");
  65     AliITS    *ITS    = (AliITS*)    gAlice->GetDetector("ITS");
  66     AliMUON   *MUON   = (AliMUON*)   gAlice->GetDetector("MUON");
  67     AliPHOS   *PHOS   = (AliPHOS*)   gAlice->GetDetector("PHOS");
  68     AliPMD    *PMD    = (AliPMD*)    gAlice->GetDetector("PMD");
  69     AliHMPID   *HMPID   = (AliHMPID*)   gAlice->GetDetector("HMPID");
  70     AliT0  *T0  = (AliT0*)  gAlice->GetDetector("T0");
  71     AliTOF    *TOF    = (AliTOF*)    gAlice->GetDetector("TOF");
  72     AliTPC    *TPC    = (AliTPC*)    gAlice->GetDetector("TPC");
  73     AliTPC    *TPC    = (AliTPC*)    gAlice->GetDetector("TPC");
  74     AliTRD    *TRD    = (AliTRD*)    gAlice->GetDetector("TRD");
  75     AliCASTOR *CASTOR = (AliCASTOR*) gAlice->GetDetector("CASTOR");
  76     AliZDC    *ZDC    = (AliZDC*)    gAlice->GetDetector("ZDC");
  77     AliEMCAL  *EMCAL  = (AliEMCAL*)  gAlice->GetDetector("EMCAL");
  78
  79 // Get pointer to the particles
  80 //    TClonesArray *Particles = gAlice->Particles();
  81 //    TParticle    *part;
  82
  83     // Create histograms
  84     if(FMD)   TH1F *hFMD   = new TH1F("hFMD"   ,"Hit Radius",100,0.,100.);
  85     if(ITS)   TH1F *hITS   = new TH1F("hITS"   ,"Ionization",100,0.,3.e-3);
  86     if(MUON)  TH1F *hMUON  = new TH1F("hMUON"  ,"Hit Radius",100,0.,500.);
  87     if(PHOS)   TH1F *hPHOS  = new TH1F("hPHOS"  ,"Energy Dep.",100,0.,0.5);
  88     if(PMD)     TH1F *hPMD   = new TH1F("hPMD"   ,"Energy Dep.",100,0.,1.e+5);
  89     if(HMPID)  TH1F *hHMPID  = new TH1F("hHMPID"  ,"Energy loss",100,0.,1.e-5);
  90     if(T0) TH1F *hT0 = new TH1F("hT0" ,"Time of Flight",100,0.,10.);
  91     if(TOF)   TH1F *hTOF   = new TH1F("hTOF"   ,"Time of Flight",100,0.,1.e-5);
  92     if(TPC)   TH1F *hTPC   = new TH1F("hTPC"   ,"Charge",100,0.,70.2);
  93     if(TRD)   TH1F *hTRD   = new TH1F("hTRD"   ,"Charge",100,0.,10.);
  94     if(CASTOR)TH1F *hCASTOR= new TH1F("hCASTOR","Hit Radius",100,0.,10.);
  95     if(ZDC)   TH1F *hZDC   = new TH1F("hZDC"   ,"Energy",100,0.,5.);
  96     if(EMCAL) TH1F *hEMCAL = new TH1F("hEMCAL" ,"Energy",100,0.,2.);
  97 //    TH1F *hTPAR  = new TH1F("hTPAR" ,"?",6,1,7);
  98     Int_t track,ntracks = gAlice->TreeH()->GetEntries();
  99 // Start loop on tracks in the hits containers
 100     for(track=0; track<ntracks;track++){
 101       //MI change
 102       gAlice->ResetHits();
 103       gAlice->TreeH()->GetEvent(track);
 104       if(FMD){
 105             for(fmdHit=(AliFMDhit*)FMD->FirstHit(-1);fmdHit;
 106                 fmdHit=(AliFMDhit*)FMD->NextHit()){
 107                 hFMD->Fill(TMath::Hypot(fmdHit->X(),fmdHit->Y()));
 108             } // end for fmdHit
 109         } // end if FMD
 110         if(ITS){
 111             for(itsHit=(AliITShit*)ITS->FirstHit(-1);itsHit;
 112                 itsHit=(AliITShit*)ITS->NextHit()){
 113                 if(itsHit->GetIonization()>0.0){//only after a step in the ITS
 114                     hITS->Fill(itsHit->GetIonization());
 115                 } // end if
 116             } // end for itsHit
 117         } // end if ITS
 118         if(MUON){
 119             for(muonHit=(AliMUONHit*)MUON->FirstHit(-1);muonHit;
 120                 muonHit=(AliMUONHit*)MUON->NextHit()){
 121                 hMUON->Fill(TMath::Hypot(muonHit->X(),muonHit->Y()));
 122             } // end for muonHit
 123         } // end if MUON
 124         if(PHOS){
 125             for(phosHit=(AliPHOSHit*)PHOS->FirstHit(-1);phosHit;
 126                 phosHit=(AliPHOSHit*)PHOS->NextHit()){
 127                 hPHOS->Fill(phosHit->GetEnergy());
 128             } // end for phosHit
 129         } // end if PHOS
 130         if(PMD){
 131             for(pmdHit=(AliPMDhit*)PMD->FirstHit(-1);pmdHit;
 132                 pmdHit=(AliPMDhit*)PMD->NextHit()){
 133                 hPMD->Fill(pmdHit->GetEnergy());
 134             } // end for pmdHit
 135         } // end if PMD
 136         if(HMPID){
 137             for(richHit=(AliHMPIDHit*)HMPID->FirstHit(-1);richHit;
 138                 richHit=(AliHMPIDHit*)HMPID->NextHit()){
 139                 hHMPID->Fill(richHit->fEloss);
 140             } // end for richHit
 141         } // end if HMPID
 142         if(T0){
 143             for(startHit=(AliT0hit*)T0->FirstHit(-1);startHit;
 144                 startHit=(AliT0hit*)T0->NextHit()){
 145                 hT0->Fill(startHit->fTime);
 146             } // end for startHit
 147         } // end if T0
 148         if(TOF){
 149             for(tofHit=(AliTOFhit*)TOF->FirstHit(-1);tofHit;
 150                 tofHit=(AliTOFhit*)TOF->NextHit()){
 151                 hTOF->Fill(tofHit->GetTof());
 152             } // end for tofHit
 153         } // end if TOF
 154
 155         if(TRD) {
 156             for(trdHit=(AliTRDhit*)TRD->FirstHit(-1);trdHit;
 157                 trdHit=(AliTRDhit*)TRD->NextHit()) {
 158                 hTRD->Fill((Float_t)(trdHit->GetCharge()));
 159             } // end for
 160         } // end if TRD
 161         if(CASTOR) {
 162             for(castorHit=(AliCASTORhit*)CASTOR->FirstHit(-1);castorHit;
 163                 castorHit=(AliCASTORhit*)CASTOR->NextHit()) {
 164                 hCASTOR->Fill(TMath::Hypot(castorHit->X(),castorHit->Y()));
 165             } // end for
 166         } // end if CASTOR
 167         if(ZDC){
 168             for(zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);zdcHit;
 169                 zdcHit=(AliZDCHit*)ZDC->NextHit()){
 170                 hZDC->Fill(zdcHit->GetEnergy());
 171             } // end for zdcdHit
 172         } // end if ZDC
 173         if(TPC) {
 174           for(tpcHit=(AliTPChit*)TPC->FirstHit(-1);tpcHit;
 175               tpcHit=(AliTPChit*)TPC->NextHit()) {
 176             hTPC->Fill((Float_t)(tpcHit->fQ));
 177           } // end for tpcHit
 178         } // end if TPC
 179         if(EMCAL) {
 180           for(emcalHit=(AliEMCALHit*)EMCAL->FirstHit(-1);emcalHit;
 181               emcalHit=(AliEMCALHit*)EMCAL->NextHit()) {
 182             hEMCAL->Fill((Float_t)(emcalHit->GetEnergy()));
 183           } // end for tpcHit
 184         } // end if TPC
 185
 186     } // end for track
 187
 188 //Create a canvas, set the view range, show histograms
 189     TCanvas *c0 = new TCanvas("c0","Alice Detectors",400,10,600,700);
 190     if(opt=="Logy")c0->SetLogy();
 191     if(FMD){
 192         hFMD->SetFillColor(42);
 193         hFMD->Draw();
 194         c0->Print("analHitsFMD.ps");
 195     } // end if FMD
 196     if(ITS){
 197         hITS->SetFillColor(42);
 198         hITS->Draw();
 199         c0->SaveAs("analHitsITS.ps");
 200     } // end if ITS
 201     if(MUON){
 202         hMUON->SetFillColor(42);
 203         hMUON->Draw();
 204         c0->SaveAs("analHitsMUON.ps");
 205     } // end if MUON
 206     if(PHOS){
 207         hPHOS->SetFillColor(42);
 208         hPHOS->Draw();
 209         c0->SaveAs("analHitsPHOS.ps");
 210     } // end if PHOS
 211     if(PMD){
 212         hPMD->SetFillColor(42);
 213         hPMD->Draw();
 214         c0->SaveAs("analHitsPMD.ps");
 215     } // end if PMD
 216     if(HMPID){
 217         hHMPID->SetFillColor(42);
 218         hHMPID->Draw();
 219         c0->SaveAs("analHitsHMPID.ps");
 220     } // end if HMPID
 221     if(T0){
 222         hT0->SetFillColor(42);
 223         hT0->Draw();
 224         c0->SaveAs("analHitsT0.ps");
 225     } // end if T0
 226     if(TOF){
 227         hTOF->SetFillColor(42);
 228         hTOF->Draw();
 229         c0->SaveAs("analHitsTOF.ps");
 230     } // end if TOF
 231     if(TPC){
 232         hTPC->SetFillColor(42);
 233         hTPC->Draw();
 234         c0->SaveAs("analHitsTPC.ps");
 235     } // end if TPC
 236     if(TRD){
 237         hTRD->SetFillColor(42);
 238         hTRD->Draw();
 239         c0->SaveAs("analHitsTRD.ps");
 240     } // end if TRD
 241     if(CASTOR){
 242         hCASTOR->SetFillColor(42);
 243         hCASTOR->Draw();
 244         c0->SaveAs("analHitsCASTOR.ps");
 245     } // end if TRD
 246     if(ZDC){
 247         hZDC->SetFillColor(42);
 248         hZDC->Draw();
 249         c0->SaveAs("analHitsZDC.ps");
 250     } // end if ZDC
 251     if(EMCAL){
 252         hEMCAL->SetFillColor(42);
 253         hEMCAL->Draw();
 254         c0->SaveAs("analHitsEMCAL.ps");
 255     } // end if ZDC
 256
 257 // Clean Up
 258     /*
 259     if(FMD)    delete hFMD;
 260     if(ITS)    delete hITS;
 261     if(MUON)   delete hMUON;
 262     if(PHOS)   delete hPHOS;
 263     if(PMD)    delete hPMD;
 264     if(HMPID)   delete hHMPID;
 265     if(T0)  delete hT0;
 266     if(TOF)    delete hTOF;
 267     if(TPC)    delete hTPC;
 268     if(TRD)    delete hTRD;
 269     if(CASTOR) delete hCASTOR;
 270     if(ZDC)    delete hZDC;
 271     */
 272 }