MUON/MUONstraggling.C

   1    TH1F *mass1 = new TH1F("mass1","Invariant Mass",120,0,12);
   2    TH1F *mass2 = new TH1F("mass2","Invariant Mass",100,9.0,10.5);
   3    TH1F *mass3 = new TH1F("mass3","Invariant Mass",100,2.5,3.5);
   4 void MUONstraggling (Int_t evNumber1=0,Int_t evNumber2=0)
   5 {
   6 /////////////////////////////////////////////////////////////////////////
   7 //   This macro is a small example of a ROOT macro
   8 //   illustrating how to read the output of GALICE
   9 //   and do some analysis.
  10 //
  11 /////////////////////////////////////////////////////////////////////////
  12
  13 // Dynamically link some shared libs
  14     static Float_t xmuon, ymuon;
  15
  16    if (gClassTable->GetID("AliRun") < 0) {
  17       gROOT->LoadMacro("loadlibs.C");
  18       loadlibs();
  19    }
  20
  21 // Connect the Root Galice file containing Geometry, Kine and Hits
  22
  23     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  24
  25
  26     if (!file) {
  27         printf("\n Creating galice.root \n");
  28         file = new TFile("galice.root");
  29     } else {
  30         printf("\n galice.root found in file list");
  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) {
  37             printf("\n Create new gAlice object");
  38             gAlice = new AliRun("gAlice","Alice test program");
  39         }
  40     }
  41
  42
  43 //  Create some histograms
  44
  45
  46    AliMUONChamber*  iChamber;
  47 //
  48 //   Loop over events
  49 //
  50    Int_t Nh=0;
  51    Int_t Nh1=0;
  52    for (Int_t nev=0; nev<= evNumber2; nev++) {
  53        Int_t nparticles = gAlice->GetEvent(nev);
  54        //cout << "nparticles  " << nparticles <<endl;
  55        if (nev < evNumber1) continue;
  56        if (nparticles <= 0) return;
  57
  58        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  59
  60
  61        TTree *TH = gAlice->TreeH();
  62        Int_t ntracks = TH->GetEntries();
  63 //
  64 //   Loop over tracks
  65 //
  66
  67        Float_t pups[4]={0,0,0,0};
  68        for (Int_t track=0; track<ntracks;track++) {
  69            Int_t Nelt=0;
  70            gAlice->ResetHits();
  71            Int_t nbytes += TH->GetEvent(track);
  72
  73
  74            if (MUON)  {
  75                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
  76                    mHit;
  77                    mHit=(AliMUONHit*)MUON->NextHit())
  78                {
  79                    Int_t   nch   = mHit->Chamber();  // chamber number
  80                    Float_t x     = mHit->X();        // x-pos of hit
  81                    Float_t y     = mHit->Y();        // y-pos
  82                    Float_t z     = mHit->Z();        // y-pos
  83                    Float_t p  = mHit->Momentum();
  84                    Float_t px = mHit->Px();
  85                    Float_t py = mHit->Py();
  86                    Float_t pz = mHit->Pz();
  87
  88                    if (nch != 1) continue;
  89
  90                    Int_t ipart = mHit->Particle();
  91                    TParticle *Part;
  92                    Int_t ftrack = mHit->Track();
  93                    Part = gAlice->Particle(ftrack);
  94                    Int_t ipart = Part->GetPdgCode();
  95                    TParticle *Mother;
  96                    Float_t px0 = Part->Px();
  97                    Float_t py0 = Part->Py();
  98                    Float_t pz0 = Part->Pz();
  99                    Float_t p0  = Part->P();
 100
 101                    Float_t e0=Part->Energy();
 102
 103                    if (ipart == kMuonPlus || ipart == kMuonMinus) {
 104 //
 105 // Branson Correction
 106 //
 107                        Float_t zch = z;
 108                        Float_t r   = TMath::Sqrt(x*x+y*y);
 109                        Float_t zb;
 110
 111                        if (r<26.3611) {
 112                            zb = 466.;
 113                        } else {
 114                            zb = 441.;
 115                        }
 116
 117                        Float_t xb  = x-(zch-zb)*px/pz;
 118                        Float_t yb  = y-(zch-zb)*py/pz;
 119                        Float_t pzB = p * zb/TMath::Sqrt(zb*zb+yb*yb+xb*xb);
 120                        Float_t pxB = pzB * xb/zb;
 121                        Float_t pyB = pzB * yb/zb;
 122                        Float_t theta
 123                            = TMath::ATan2(TMath::Sqrt(pxB*pxB+pyB*pyB), pzB) *
 124                            180./TMath::Pi();
 125
 126
 127 //
 128 // Energy Correction
 129 //
 130 //
 131
 132                        Float_t corr=(p+CorrectP(p, theta))/p;
 133                        pups[0]+=p*corr;
 134                        pups[1]+=pxB*corr;
 135                        pups[2]+=pyB*corr;
 136                        pups[3]+=pzB*corr;
 137                    }
 138                } // hits
 139            } // if MUON
 140        } // tracks
 141        Float_t mass =
 142            TMath::Sqrt(pups[0]*pups[0]
 143                        -pups[1]*pups[1]-pups[2]*pups[2]-pups[3]*pups[3]);
 144        mass1->Fill(mass, 1.);
 145        mass2->Fill(mass, 1.);
 146        mass3->Fill(mass, 1.);
 147     } // event
 148 //Create a canvas, set the view range, show histograms
 149    TCanvas *c1 = new TCanvas("c1","Mass Spectra",400,10,600,700);
 150
 151    mass2->SetFillColor(42);
 152    mass2->SetXTitle("Mass (GeV)");
 153    mass2->Draw();
 154 }
 155
 156
 157
 158 Float_t CorrectP(Float_t p, Float_t theta)
 159 {
 160     Float_t c;
 161
 162     if (theta<3.) {
 163 //W
 164         if (p<15) {
 165             c = 2.916+0.0156*p-0.00000866*p*p;
 166         } else {
 167             c = 2.998+0.0137*p;
 168         }
 169     } else {
 170 //Cu+Fe+Cc+C
 171         if (p<15) {
 172             c = 2.464+0.00877*p-0.0000106*p*p;
 173         } else {
 174             c = 2.496+0.00817*p;
 175         }
 176     }
 177 //
 178 //
 179     c*=0.95;
 180 //
 181     return c;
 182 }
 183
 184
 185
 186
 187
 188