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       gROOT->LoadMacro("loadlibs.C");
  20       loadlibs();
  21    } else {
  22       delete gAlice;
  23       gAlice = 0;
  24    }
  25 //
  26 // Connect the Root Galice file containing Geometry, Kine and Hits
  27    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  28    if (!file) file = new TFile("galice.root");
  29
  30 // Get AliRun object from file or create it if not on file
  31    if (!gAlice) {
  32       gAlice = (AliRun*)file->Get("gAlice");
  33       if (gAlice) printf("AliRun object found on file\n");
  34       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  35    }
  36    TClonesArray *PartArray = gAlice->Particles();
  37    TParticle *Part;
  38
  39
  40 // Import the Kine and Hits Trees for the event evNumber in the file
  41    Int_t nparticles = gAlice->GetEvent(evNumber);
  42    if (nparticles <= 0) return;
  43 //
  44 //
  45    TH1F *dmass = new TH1F("dmass","Dimuon-Mass Distribution"
  46                          ,25,0.,5.);
  47    TH1F *dmassc = new TH1F("dmassc","Dimuon-Mass Distribution"
  48                          ,50,0.,10.);
  49    TH1F *dmassd = new TH1F("dmassd","Dimuon-Mass Distribution"
  50                          ,50,0.,10.);
  51
  52    TH1F *pt    = new TH1F("pt","pT-single"
  53                          ,50,0.,10.);
  54 //
  55 //  Generator Loop
  56 //
  57
  58    AliGenCocktailEntry *Entry, *e1, *e2;
  59    AliGenCocktail*  Cocktail = (AliGenCocktail*) gAlice->Generator();
  60 //                   Single Generator
  61    for (Entry=Cocktail->FirstGenerator();
  62         Entry;
  63         Entry=Cocktail->NextGenerator()
  64        ) {
  65        Entry->PrintInfo();
  66      }
  67 //                   Pairs of Generators
  68 //
  69 // Initialize Combinator
  70    AliDimuCombinator* Combinator = new AliDimuCombinator(PartArray);
  71    Combinator->SetEtaCut(-5, 5);
  72    Combinator->SetPtMin(0.0);
  73
  74    Int_t i=0;
  75 //
  76 // Single Muon  Loop
  77 //
  78    Combinator->ResetRange();
  79
  80
  81
  82    for(Muon=Combinator->FirstMuon();
  83        Muon;
  84        Muon=Combinator->NextMuon()) {
  85 //
  86        Int_t chfirst= Muon->GetFirstDaughter();
  87        Int_t chlast = Muon->GetLastDaughter();
  88        Int_t parent = Muon->GetFirstMother();
  89        Float_t ptm  = Muon->Pt();
  90        Float_t eta  = Muon->Eta();
  91
  92        printf("\n Particle %d Parent %d first child %d last child %d",
  93               i,parent, chfirst, chlast);
  94        printf("\n Particle pt, eta: %f , %f ",pt,eta);
  95        i++;
  96    }
  97 //
  98 // Di-Muon Loop
  99
 100    Float_t pt1,pt2;
 101    TParticle* Muon1, *Muon2;
 102
 103    Combinator->ResetRange();
 104 /*
 105    for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
 106         (Muon1 && Muon2);
 107         Combinator->NextMuonPairSelected(Muon1,Muon2))
 108    {
 109        pt1=Muon1->Pt();
 110        pt2=Muon2->Pt();
 111        Float_t mass=Combinator->Mass(Muon1, Muon2);
 112        Float_t wgt =Combinator->Weight(Muon1, Muon2);
 113        pt->Fill(pt1, wgt);
 114        pt->Fill(pt2, wgt);
 115        Float_t smeared_mass=mass;
 116        Combinator->SmearGauss(0.05*mass, smeared_mass);
 117        if (Combinator->Correlated(Muon1, Muon2)) {
 118            dmassc->Fill(mass, wgt);
 119        } else {
 120            dmass->Fill(mass, wgt);
 121        }
 122     }
 123 */
 124 //
 125 //   Dimuon Loop controlled by Generator Loop
 126 //
 127    for (Cocktail->FirstGeneratorPair(e1,e2);
 128         (e1&&e2);
 129         Cocktail->NextGeneratorPair(e1,e2)
 130        ){
 131        printf("\n ----------------------------------------------------");
 132        e1->PrintInfo();
 133        e2->PrintInfo();
 134        printf("\n ----------------------------------------------------");
 135        Combinator->SetFirstRange (e1->GetFirst(), e1->GetLast());
 136        Combinator->SetSecondRange(e2->GetFirst(), e2->GetLast());
 137        Combinator->SetRate(e1->Rate(), e2->Rate());
 138        for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
 139             (Muon1 && Muon2);
 140             Combinator->NextMuonPairSelected(Muon1,Muon2))
 141        {
 142            pt1=Muon1->Pt();
 143            pt2=Muon2->Pt();
 144            Float_t mass=Combinator->Mass(Muon1, Muon2);
 145            Float_t wgt =Combinator->Weight(Muon1, Muon2);
 146            pt->Fill(pt1, wgt);
 147            pt->Fill(pt2, wgt);
 148            Float_t smeared_mass=mass;
 149            Combinator->SmearGauss(0.05*mass, smeared_mass);
 150            Float_t DecayWeight=
 151                Combinator->Decay_Prob(Muon2)*
 152                Combinator->Decay_Prob(Muon2);
 153 //         if (TMath::Min(pt1,pt2) > -0.5*TMath::Max(pt1,pt2)+2.) {
 154                if (Combinator->Correlated(Muon1, Muon2)) {
 155                    dmassc->Fill(smeared_mass, wgt);
 156                }
 157                    dmass->Fill(smeared_mass, wgt);
 158                    dmassd->Fill(smeared_mass, wgt*DecayWeight);
 159 //         }
 160        } // Dimuon Loop
 161    }// Generator Loop
 162 //
 163 //Create a canvas, set the view range, show histograms
 164 //
 165    TCanvas *c1 = new TCanvas("c1","Dimuon Plots",400,10,600,700);
 166    TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 167    TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 168    TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 169    TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 170    pad11->SetFillColor(11);
 171    pad12->SetFillColor(11);
 172    pad13->SetFillColor(11);
 173    pad14->SetFillColor(11);
 174    pad11->Draw();
 175    pad12->Draw();
 176    pad13->Draw();
 177    pad14->Draw();
 178
 179    pad11->cd();
 180    pt->SetFillColor(42);
 181    pt->SetXTitle("pT (GeV)");
 182    pt->Draw();
 183
 184    pad12->cd();
 185    pad12->SetLogy(0);
 186    dmass->SetFillColor(42);
 187    dmass->SetXTitle("m (GeV)");
 188    dmass->Draw();
 189
 190    pad13->cd();
 191    pad13->SetLogy(0);
 192    dmassc->SetFillColor(42);
 193    dmassc->SetXTitle("m (GeV)");
 194    dmassc->Draw();
 195
 196    pad14->cd();
 197    pad14->SetLogy(0);
 198    dmassd->SetFillColor(42);
 199    dmassd->SetXTitle("m (GeV)");
 200    dmassd->Draw();
 201
 202
 203 }
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213