MUON/MUONcombi.C

   1 void MUONcombi (Int_t evNumber=0)
   2 {
   3     const Float_t runWeight = 4.e8;
   4
   5 // Dynamically link some shared libs
   6    if (gClassTable->GetID("AliRun") < 0) {
   7       gROOT->LoadMacro("loadlibs.C");
   8       loadlibs();
   9    } else {
  10       delete gAlice;
  11       gAlice = 0;
  12    }
  13
  14 // Connect the Root Galice file containing Geometry, Kine and Hits
  15    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  16    if (!file) file = new TFile("galice.root");
  17
  18 // Get AliRun object from file or create it if not on file
  19    if (!gAlice) {
  20       gAlice = (AliRun*)file->Get("gAlice");
  21       if (gAlice) printf("AliRun object found on file\n");
  22       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  23    }
  24
  25 // Import the Kine and Hits Trees for the event evNumber in the file
  26    Int_t nparticles = gAlice->GetEvent(evNumber);
  27    if (nparticles <= 0) return;
  28 //
  29 //
  30    TH1F *dmass  = new TH1F("dmass","Dimuon-Mass Distribution"
  31                           ,100,0.,5.);
  32
  33    TH1F *dmassc = new TH1F("dmassc","Dimuon-Mass Distribution"
  34                            ,200,0.,10.);
  35    TH1F *dmassd = new TH1F("dmassd","Dimuon-Mass Distribution"
  36                            ,200,0.,10.);
  37
  38    TH1F *pt     = new TH1F("pt","pT-single"
  39                            ,50,0.,10.);
  40    TH1F *hCont[30];
  41
  42 //
  43 //  Generator Loop
  44 //
  45
  46    Int_t i=0;
  47
  48    AliGenCocktailEntry *Entry, *e1, *e2;
  49    AliGenCocktail*  Cocktail = (AliGenCocktail*) gAlice->Generator();
  50 //                   Single Generator
  51    for (Entry=Cocktail->FirstGenerator();
  52         Entry;
  53         Entry=Cocktail->NextGenerator()
  54        ) {
  55        Entry->PrintInfo();
  56      }
  57 //
  58 // Initialize Combinator
  59 //
  60    AliDimuCombinator* Combinator = new AliDimuCombinator();
  61    Combinator->SetEtaCut(2.5, 4.);
  62    Combinator->SetPtMin(1.0);
  63
  64    i=0;
  65
  66 //
  67 // Single Muon  Loop
  68 //
  69    Combinator->ResetRange();
  70
  71
  72
  73    for(Muon=Combinator->FirstMuon();
  74        Muon;
  75        Muon=Combinator->NextMuon()) {
  76 //
  77        Int_t chfirst = Muon->GetFirstDaughter();
  78        Int_t chlast  = Muon->GetLastDaughter();
  79        Int_t parent  = Muon->GetFirstMother();
  80        Float_t ptm   = Muon->Pt();
  81        Float_t eta   = Muon->Eta();
  82
  83 //       printf("\n Particle %d Parent %d first child %d last child %d",
  84 //            i,parent, chfirst, chlast);
  85 //       printf("\n Particle pt, eta: %f , %f ", ptm, eta);
  86        i++;
  87    }
  88 //
  89 // Di-Muon Loop
  90
  91    Float_t pt1,pt2;
  92    TParticle* Muon1;
  93    TParticle* Muon2;
  94
  95    Combinator->ResetRange();
  96
  97 //
  98 //   Dimuon Loop controlled by Generator Loop
  99 //
 100    Float_t  sig = 0;
 101    Int_t  icont = 0;
 102    char name[30];
 103    for (Cocktail->FirstGeneratorPair(e1,e2);
 104         (e1&&e2);
 105         Cocktail->NextGeneratorPair(e1,e2)
 106        ){
 107
 108        sprintf(name, "%s-%s", e1->GetName(), e2->GetName());
 109        hCont[icont] = new TH1F(name,"Dimuon-Mass Distribution",100,0.,5.);
 110        printf("\n ----------------------------------------------------");
 111        e1->PrintInfo();
 112        e2->PrintInfo();
 113        printf("\n ----------------------------------------------------");
 114        Combinator->SetFirstRange (e1->GetFirst(), e1->GetLast());
 115        Combinator->SetSecondRange(e2->GetFirst(), e2->GetLast());
 116        Combinator->SetRate(e1->Rate(), e2->Rate());
 117        for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
 118             (Muon1 && Muon2);
 119             Combinator->NextMuonPairSelected(Muon1,Muon2))
 120        {
 121            pt1 = Muon1->Pt();
 122            pt2 = Muon2->Pt();
 123            Float_t mass = Combinator->Mass(Muon1, Muon2);
 124            Float_t wgt  = runWeight*Combinator->Weight(Muon1, Muon2)*
 125                Combinator->DecayProbability(Muon1)*
 126                Combinator->DecayProbability(Muon2);
 127            pt->Fill(pt1, wgt);
 128            pt->Fill(pt2, wgt);
 129            Float_t smeared_mass = mass;
 130            Combinator->SmearGauss(0.02*mass, smeared_mass);
 131
 132            if (TMath::Min(pt1,pt2) > -0.5*TMath::Max(pt1,pt2)+2.) {
 133                if (Combinator->Correlated(Muon1, Muon2)) {
 134                    dmassc->Fill(smeared_mass, wgt);
 135                    sig += wgt;
 136                } else {
 137 // account for the fact that we sum like-sign and unlike-sign
 138                    wgt *= 0.5;
 139                }
 140                dmass->Fill(smeared_mass, wgt);
 141                hCont[icont]->Fill(smeared_mass, wgt);
 142            }
 143        } // Dimuon Loop
 144        icont++;
 145    }// Generator Loop
 146    printf("\n Signal %e \n \n", sig);
 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    pad12->SetLogy(0);
 172    dmass->SetFillColor(42);
 173    dmass->SetXTitle("m (GeV)");
 174    dmass->Draw();
 175
 176    pad13->cd();
 177    pad13->SetLogy(0);
 178    dmassc->SetFillColor(42);
 179    dmassc->SetXTitle("m (GeV)");
 180    dmassc->Draw();
 181
 182    pad14->cd();
 183    pad14->SetLogy(0);
 184    dmassd->SetFillColor(42);
 185    dmassd->SetXTitle("m (GeV)");
 186    dmassd->Draw();
 187
 188
 189 }
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200