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