ACORDE/macros/TestACORDEHits.C

   1 /////////////////////////////////////////////////////////////////////////
   2 //   This macro is a small example of a ROOT macro
   3 //   illustrating how to read the hits in ACORDE,
   4 //   fill some histograms and store it on a file
   5 //
   6 //   Author:
   7 //
   8 //      Eleazar Cuautle <ecuautle@nucleares.unam.mx>
   9 //
  10 //   Also:
  11 //
  12 //      Mario Rodriguez Cahuantzi <mrodrigu@mail.cern.ch>
  13 //
  14 /////////////////////////////////////////////////////////////////////////
  15
  16 void TestACORDEHits (const char *filename="galice.root",Int_t numberOfEvents=0){
  17    // delete existing gAlice
  18   // Dynamically link some shared libs
  19   if (gClassTable->GetID("AliRun") < 0)
  20     {
  21       gROOT->LoadMacro("loadlibs.C");
  22       loadlibs();
  23     }
  24   if (gAlice)
  25     {
  26       //delete gAlice->GetRunLoader();
  27       delete gAlice;
  28       gAlice = 0x0;
  29     }
  30
  31   AliRunLoader *rl = AliRunLoader::Open("galice.root",
  32                                         AliConfig::GetDefaultEventFolderName(),"read");
  33   if (!rl)
  34     {
  35       cerr<<"Can't load RunLoader from file! \n";      return 0x0;
  36     }
  37
  38   rl->LoadgAlice();
  39   gAlice = rl->GetAliRun();
  40   if (!gAlice)
  41     {
  42       cerr << " AliRun object not found on file \n";     return 0x0;
  43     }
  44   rl->LoadHeader();
  45
  46   // Get the pointer to the ACORDE detector
  47   AliLoader *acordel = rl->GetLoader("ACORDELoader");
  48   AliACORDE *ACORDE = (AliACORDE *) gAlice->GetDetector("ACORDE");
  49   if (ACORDE == 0x0 || acordel == 0x0) {
  50     cerr << " Can not find ACORDE or ACORDELoader \n";    return 0x0;
  51   }
  52   //Determine how many events would you like to analyze
  53   if (numberOfEvents==0) numberOfEvents=(Int_t)(rl->GetNumberOfEvents());
  54
  55   //create histograms
  56
  57   TH1F  *fACORDEEnergy = new TH1F("fACORDEEnergy","Energy distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Energy (GeV)",1200,0,360);
  58   TH1F  *fACORDEEnergyLoss = new TH1F("fACORDEEnergyLoss","Loss energy distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);Energy loss (Gev)",100,0,0.030);
  59   TH1F  *fACORDEAzimuth = new TH1F("fACORDEAzimuth","Azimuth angle distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Azimuth angle",360,-180,180);
  60   TH1F  *fACORDEPolar = new TH1F("fACORDEPolar","Polar angle distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Polar angle",180,0,180);
  61   TH2F  *fACORDExz = new TH2F("fACORDExz" ,"Distribution XZ  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); X (cms); Z (cms) ",900,-850.,850.,1200,-500.,500.);
  62   TH2F  *fACORDEAzimPol =  new TH2F("fACORDEAzimPol" ,"Azimuth vs Polar angle distribution ",360,-180.,180.,180,0.,180.);
  63   TH1F  *fACORDEPx = new TH1F("fACORDEPx" ,"P_{x} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{x}",60,-120.,120.);
  64   TH1F  *fACORDEPy = new TH1F("fACORDEPy" ,"P_{y} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{y}",60,-120.,120.);
  65   TH1F  *fACORDEPz = new TH1F("fACORDEPz" ,"P_{z} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{z}",60,-120.,120.);
  66   TH1F  *fACORDEPt = new TH1F("fACORDEPt" ,"P_{t} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{t}",1200,0.,360.);
  67   TH2F  *fELossAzimuth = new TH2F("fELossAzimuth","Dist. Azimuth Angle VS Energy Loss; Azimuth angle; Energy (GeV)",360,0,350,100,0,0.030);
  68   TH2F  *fEnergyAzimuth = new TH2F("fEnergyAzimuth","Dist. Azimuth Angle VS Energy; Azimuth angle; Energy (GeV)",360,0,360,300,0,350);
  69   TH2F  *fELossPolar = new TH2F("fELossPolar","Dist. Polar Angle VS Energy Loss; Polar angle; Energy (GeV)",180,0,180,100,0,0.030);
  70   TH2F  *fEnergyPolar = new TH2F("fEnergyPolar","Dist. Polar Angle VS Energy; Polar angle; Energy (GeV)",180,0,180,300,0,350);
  71
  72  //------------------------------------------------------------------//
  73
  74   for (Int_t ievent=0; ievent<numberOfEvents; ievent++) {
  75     if ((ievent%10) == 0)  printf ("Processing event %d \n", ievent);
  76     rl->GetEvent(ievent);
  77
  78     // Get the pointer Hit tree
  79     acordel->LoadHits();
  80     TTree *hitTree = acordel->TreeH();
  81     ACORDE->SetTreeAddress();
  82     if (!hitTree) {
  83       cout << " No TreeH found" << endl;      return 0x0; //rc;
  84     }
  85
  86     rl->LoadKinematics();
  87     Int_t nTrack = (Int_t) hitTree->GetEntries();
  88     //    cout <<"<AliACORDEanalyzeHits> Found "<< nTrack <<" primary particles with hits \n";
  89
  90     // Start loop on tracks in the hits containers
  91     for(Int_t iTrack=0; iTrack<nTrack;iTrack++){
  92       ACORDE->ResetHits();
  93       hitTree->GetEvent(iTrack);
  94       if(ACORDE) {
  95         for(acordeHit=(AliACORDEhit*)ACORDE->FirstHit(-1);acordeHit;
  96             acordeHit=(AliACORDEhit*)ACORDE->NextHit()) {
  97                 fACORDEEnergy->Fill(acordeHit->Energy());
  98                 fACORDEEnergyLoss->Fill(acordeHit->Eloss());
  99                 fACORDEAzimuth->Fill(acordeHit->AzimuthAngle());
 100                 fACORDEPolar->Fill(acordeHit->PolarAngle());
 101                 fACORDExz->Fill(acordeHit->X(),acordeHit->Z());
 102                 fELossAzimuth->Fill((Float_t)(acordeHit->AzimuthAngle()),(Float_t)(acordeHit->Eloss()));
 103                 fACORDEPx->Fill( (Float_t)(acordeHit->Px()));
 104                 fACORDEPy->Fill( (Float_t)(acordeHit->Py()));
 105                 fACORDEPz->Fill( (Float_t)(acordeHit->Pz()));
 106                 fACORDEPt->Fill( TMath::Sqrt( (acordeHit->Px())*(acordeHit->Px())+
 107                                              (acordeHit->Py())*(acordeHit->Py())) );
 108                 fACORDEAzimPol->Fill( (Float_t)(acordeHit->AzimuthAngle()),
 109                                 (Float_t)(acordeHit->PolarAngle()));
 110                 fEnergyAzimuth->Fill(acordeHit->Energy(),acordeHit->AzimuthAngle());
 111                 fEnergyPolar->Fill(acordeHit->Energy(),acordeHit->PolarAngle());
 112                 fELossPolar->Fill(acordeHit->PolarAngle(),acordeHit->Eloss());
 113         //cout << "Muon triggered by ACORDE "<<endl;
 114         //cout << "Energy: "<<acordeHit->Energy()<<" EnergyLoss: "<<acordeHit->Eloss()<<" Polar angle: "<<acordeHit->PolarAngle()<<" Azimuth angle: "<<acordeHit->AzimuthAngle()<<endl;
 115
 116         } // end for acordeHits
 117       } // end if ACORDE
 118     } // end for iTrack
 119   }
 120
 121  // Put the histograms in a TCanvas
 122
 123         TCanvas *c1 = new TCanvas("c1","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
 124         c1->Divide(2,1);
 125
 126         c1->cd(1)->SetLogy();
 127         fACORDEEnergy->SetMarkerStyle(kFullCircle);
 128         fACORDEEnergy->SetMarkerSize(0.7);
 129         fACORDEEnergy->Draw("E");
 130
 131         c1->cd(2)->SetLogy();
 132         fACORDEEnergyLoss->SetMarkerStyle(kFullCircle);
 133         fACORDEEnergyLoss->SetMarkerSize(0.7);
 134         fACORDEEnergyLoss->Draw("E");
 135
 136         c1->Draw();
 137
 138         TCanvas *c2 = new TCanvas("c2","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
 139         c2->Divide(2,2);
 140         c2->cd(1);
 141         fACORDEAzimuth->Draw();
 142
 143         c2->cd(2);
 144         fACORDEPolar->Draw();
 145
 146         c2->cd(3);
 147         gStyle->SetPalette(1);
 148         fACORDExz->DrawCopy("colz");
 149
 150         c2->cd(4);
 151         gStyle->SetPalette(1);
 152         fACORDEAzimPol->DrawCopy("colz");
 153
 154         c2->Draw();
 155
 156         TCanvas *c3 = new TCanvas("c3","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
 157         c3->Divide(2,2);
 158
 159         c3->cd(1)->SetLogy();
 160         fACORDEPx->SetMarkerStyle(kFullCircle);
 161         fACORDEPx->SetMarkerSize(0.7);
 162         fACORDEPx->Draw("E");
 163
 164         c3->cd(2)->SetLogy();
 165         fACORDEPy->SetMarkerStyle(kFullCircle);
 166         fACORDEPy->SetMarkerSize(0.7);
 167         fACORDEPy->Draw("E");
 168
 169         c3->cd(3)->SetLogy();
 170         fACORDEPz->SetMarkerStyle(kFullCircle);
 171         fACORDEPz->SetMarkerSize(0.7);
 172         fACORDEPz->Draw("E");
 173
 174         c3->cd(4)->SetLogy();
 175         fACORDEPt->SetMarkerStyle(kFullCircle);
 176         fACORDEPt->SetMarkerSize(0.7);
 177         fACORDEPt->Draw("E");
 178
 179         c3->Draw();
 180
 181         TCanvas *c4 = new TCanvas("c4","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
 182         c4->Divide(2,2);
 183
 184         c4->cd(1);
 185         gStyle->SetPalette(1);
 186         fEnergyPolar->DrawCopy("colz");
 187
 188         c4->cd(2);
 189         gStyle->SetPalette(1);
 190         fEnergyAzimuth->DrawCopy("colz");
 191
 192         c4->cd(3);
 193         gStyle->SetPalette(1);
 194         fELossPolar->DrawCopy("colz");
 195
 196         c4->cd(4);
 197         gStyle->SetPalette(1);
 198         fELossAzimuth->DrawCopy("colz");
 199
 200         c4->Draw();
 201
 202
 203  // save histos in a root file
 204   TFile *fout = new TFile("ACORDE_hits.root","RECREATE");
 205         fACORDEEnergy->Write();
 206         fACORDEEnergyLoss->Write();
 207         fACORDEAzimuth->Write();
 208         fACORDEPolar->Write();
 209         fACORDExz->Write();
 210         fELossAzimuth->Write();
 211         fACORDEPx->Write();
 212         fACORDEPy->Write();
 213         fACORDEPz->Write();
 214         fACORDEPt->Write();
 215         fACORDEAzimPol->Write();
 216         c1->Write();
 217         c2->Write();
 218         c3->Write();
 219         c4->Write();
 220 }