MUON/MUONcombi.C

   1 void MUONcombi (Int_t evNumber=0)
   2 {
   3 /////////////////////////////////////////////////////////////////////////
   4 //   This macro is a small example of a ROOT macro
   5 //   illustrating how to read the output of GALICE
   6 //   and fill some histograms.
   7 //
   8 //     Root > .L anal.C   //this loads the macro in memory
   9 //     Root > anal();     //by default process first event
  10 //     Root > anal(2);    //process third event
  11 //Begin_Html
  12 /*
  13 <img src="gif/anal.gif">
  14 */
  15 //End_Html
  16 /////////////////////////////////////////////////////////////////////////
  17 // Dynamically link some shared libs
  18     if (gClassTable->GetID("AliRun") < 0) {
  19         gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
  20         gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
  21         gSystem->Load("$ROOTSYS/lib/libEG.so");
  22         gSystem->Load("$ROOTSYS/lib/libEGPythia.so");
  23         gSystem->Load("libGeant3Dummy.so");    //a dummy version of Geant3
  24         gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes
  25         gSystem->Load("libgalice.so");         // the standard Alice classes
  26     }
  27 //
  28 // Connect the Root Galice file containing Geometry, Kine and Hits
  29    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  30    if (!file) file = new TFile("galice.root");
  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    }
  38 //   TClonesArray &Partarray = *(gAlice->Particles());
  39    TClonesArray *PartArray = gAlice->Particles();
  40    GParticle *Part;
  41
  42
  43 // Import the Kine and Hits Trees for the event evNumber in the file
  44    Int_t nparticles = gAlice->GetEvent(evNumber);
  45    if (nparticles <= 0) return;
  46 //
  47 //
  48    TH1F *dmass = new TH1F("dmass","Dimuon-Mass Distribution"
  49                          ,50,0.,10.);
  50    TH1F *pt    = new TH1F("pt","pT-single"
  51                          ,50,0.,10.);
  52 //
  53 //  Generator Loop
  54 //
  55
  56
  57    AliGenCocktailEntry *Entry, *e1, *e2;
  58    AliGenCocktail*  Cocktail = (AliGenCocktail*) gAlice->Generator();
  59 //                   Single Generator
  60    for (Entry=Cocktail->FirstGenerator();
  61         Entry;
  62         Entry=Cocktail->NextGenerator()
  63        ) {
  64        Entry->PrintInfo();
  65      }
  66 //                   Pairs of Generators
  67
  68 //
  69 // Initialize Combinator
  70    DimuonCombinator Combinator = DimuonCombinator(PartArray);
  71    Combinator->SetEtaCut(-4, 4);
  72    Combinator->SetPtMin(0.5);
  73
  74    Int_t i=0;
  75 //
  76 // Single Muon  Loop
  77 //
  78 /*
  79
  80    for(Muon=Combinator.FirstMuon();
  81        Muon;
  82        Muon=Combinator.NextMuon()) {
  83 //
  84        Int_t chfirst= Muon->GetFirstChild();
  85        Int_t chlast = Muon->GetLastChild();
  86        Int_t parent = Muon->GetParent();
  87        Float_t ptm  = Muon->GetPT();
  88        Float_t eta  = Muon->GetEta();
  89
  90        printf("\n Particle %d Parent %d first child %d last child %d",
  91               i,parent, chfirst, chlast);
  92        printf("\n Particle pt, eta: %f , %f ",pt,eta);
  93        i++;
  94    }
  95 */
  96 //
  97 // Di-Muon Loop
  98
  99    Float_t pt1,pt2;
 100    GParticle* Muon1, *Muon2;
 101 /*
 102    Combinator->ResetRange();
 103    for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
 104         (Muon1 && Muon2);
 105         Combinator->NextMuonPairSelected(Muon1,Muon2))
 106    {
 107        pt1=Muon1->GetPT();
 108        pt2=Muon2->GetPT();
 109        Float_t mass=Combinator->Mass(Muon1, Muon2);
 110        Float_t wgt =Combinator->Weight(Muon1, Muon2);
 111        pt->Fill(pt1, wgt);
 112        pt->Fill(pt2, wgt);
 113        Float_t smeared_mass=mass;
 114        Combinator->SmearGauss(0.05*mass, smeared_mass);
 115        dmass->Fill(smeared_mass , wgt);
 116     }
 117 */
 118 //
 119 //   Dimuon Loop controlled by Generator Loop
 120 //
 121    for (Cocktail->FirstGeneratorPair(e1,e2);
 122         (e1&&e2);
 123         Cocktail->NextGeneratorPair(e1,e2)
 124        ){
 125        printf("\n ----------------------------------------------------");
 126        e1->PrintInfo();
 127        e2->PrintInfo();
 128        printf("\n ----------------------------------------------------");
 129        Combinator->SetFirstRange (e1->GetFirst(), e1->GetLast());
 130        Combinator->SetSecondRange(e2->GetFirst(), e2->GetLast());
 131        Combinator->SetRate(e1->Rate(), e2->Rate());
 132        for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
 133             (Muon1 && Muon2);
 134             Combinator->NextMuonPairSelected(Muon1,Muon2))
 135        {
 136            pt1=Muon1->GetPT();
 137            pt2=Muon2->GetPT();
 138            Float_t mass=Combinator->Mass(Muon1, Muon2);
 139            Float_t wgt =Combinator->Weight(Muon1, Muon2);
 140            pt->Fill(pt1, wgt);
 141            pt->Fill(pt2, wgt);
 142            Float_t smeared_mass=mass;
 143            Combinator->SmearGauss(0.05*mass, smeared_mass);
 144            dmass->Fill(smeared_mass , wgt);
 145        } // Dimuon Loop
 146    }// Generator Loop
 147
 148 //
 149 //Create a canvas, set the view range, show histograms
 150 //
 151    TCanvas *c1 = new TCanvas("c1","Dimuon Plots",400,10,600,700);
 152    TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 153    TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 154    TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 155    TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 156    pad11->SetFillColor(11);
 157    pad12->SetFillColor(11);
 158    pad13->SetFillColor(11);
 159    pad14->SetFillColor(11);
 160    pad11->Draw();
 161    pad12->Draw();
 162    pad13->Draw();
 163    pad14->Draw();
 164
 165    pad11->cd();
 166    pt->SetFillColor(42);
 167    pt->SetXTitle("pT (GeV)");
 168    pt->Draw();
 169
 170    pad12->cd();
 171    dmass->SetFillColor(42);
 172    dmass->SetXTitle("m (GeV)");
 173    dmass->Draw();
 174 }
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184