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     file->ls();
  33
  34 // Get AliRun object from file or create it if not on file
  35     if (!gAlice) {
  36         gAlice = (AliRun*)(file->Get("gAlice"));
  37         if (gAlice) printf("AliRun object found on file\n");
  38         if (!gAlice) {
  39             printf("\n Create new gAlice object");
  40             gAlice = new AliRun("gAlice","Alice test program");
  41         }
  42     }
  43
  44
  45 //  Create some histograms
  46
  47
  48    AliMUONChamber*  iChamber;
  49 //
  50 //   Loop over events
  51 //
  52    Int_t Nh=0;
  53    Int_t Nh1=0;
  54    for (Int_t nev=0; nev<= evNumber2; nev++) {
  55        Int_t nparticles = gAlice->GetEvent(nev);
  56        //cout << "nparticles  " << nparticles <<endl;
  57        if (nev < evNumber1) continue;
  58        if (nparticles <= 0) return;
  59
  60        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  61
  62
  63        TTree *TH = gAlice->TreeH();
  64        Int_t ntracks = TH->GetEntries();
  65 //
  66 //   Loop over tracks
  67 //
  68
  69        Float_t pups[4]={0,0,0,0};
  70        for (Int_t track=0; track<ntracks;track++) {
  71            Int_t Nelt=0;
  72            gAlice->ResetHits();
  73            Int_t nbytes += TH->GetEvent(track);
  74
  75
  76            if (MUON)  {
  77                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
  78                    mHit;
  79                    mHit=(AliMUONHit*)MUON->NextHit())
  80                {
  81                    Int_t   nch   = mHit->fChamber;  // chamber number
  82                    Float_t x     = mHit->fX;        // x-pos of hit
  83                    Float_t y     = mHit->fY;        // y-pos
  84                    Float_t z     = mHit->fZ;        // y-pos
  85                    Float_t p=mHit->fPTot;
  86                    Float_t px=mHit->fCxHit;
  87                    Float_t py=mHit->fCyHit;
  88                    Float_t pz=mHit->fCzHit;
  89
  90                    if (nch != 1) continue;
  91
  92                    Int_t ipart = mHit->fParticle;
  93                    TClonesArray *fPartArray = gAlice->Particles();
  94                    TParticle *Part;
  95                    Int_t ftrack = mHit->fTrack;
  96                    Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
  97                    Int_t ipart = Part->GetPdgCode();
  98                    TParticle *Mother;
  99                    Float_t px0=Part->Px();
 100                    Float_t py0=Part->Py();
 101                    Float_t pz0=Part->Pz();
 102                    Float_t e0=Part->Energy();
 103
 104                    if (ipart == kMuonPlus || ipart == kMuonMinus) {
 105 //
 106 // Branson Correction
 107 //
 108                        Float_t zch=505.;
 109                        Float_t r=TMath::Sqrt(x*x+y*y);
 110                        Float_t zb;
 111
 112                        if (r<26.3611) {
 113                            zb=466.;
 114                        } else {
 115                            zb=441.;
 116                        }
 117
 118                        Float_t xb=x-(zch-zb)*px/pz;
 119                        Float_t yb=y-(zch-zb)*py/pz;
 120                        Float_t tx=xb/zb;
 121                        Float_t ty=yb/zb;
 122                        pz=zb*p/TMath::Sqrt(zb*zb+xb*xb+yb*yb);
 123                        px=pz*tx;
 124                        py=pz*ty;
 125
 126 //
 127 // Energy Correction
 128 //
 129 //
 130                        Float_t corr=(p+CorrectP(p,mHit->fTheta))/p;
 131                        pups[0]+=p*corr;
 132                        pups[1]+=px*corr;
 133                        pups[2]+=py*corr;
 134                        pups[3]+=pz*corr;
 135                    }
 136                } // hits
 137            } // if MUON
 138        } // tracks
 139        Float_t mass=TMath::Sqrt(pups[0]*pups[0]-pups[1]*pups[1]-pups[2]*pups[2]-pups[3]*pups[3]);
 140        mass1->Fill(mass, 1.);
 141        mass2->Fill(mass, 1.);
 142        mass3->Fill(mass, 1.);
 143     } // event
 144 //Create a canvas, set the view range, show histograms
 145    TCanvas *c1 = new TCanvas("c1","Vertices from electrons and positrons",400,10,600,700);
 146    c1->Divide(2,2);
 147
 148    c1->cd(1);
 149    mass1->SetFillColor(42);
 150    mass1->SetXTitle("Mass (GeV)");
 151    mass1->Draw();
 152
 153    c1->cd(2);
 154    mass2->SetFillColor(42);
 155    mass2->SetXTitle("Mass (GeV)");
 156    mass2->Draw();
 157
 158    c1->cd(3);
 159    mass3->SetFillColor(42);
 160    mass3->SetXTitle("Mass (GeV)");
 161    mass3->Draw();
 162
 163 }
 164
 165
 166
 167 Float_t CorrectP(Float_t p, Float_t theta)
 168 {
 169     if (theta<3.) {
 170 //W
 171         if (p<15) {
 172             return 2.737+0.0494*p-0.001123*p*p;
 173         } else {
 174             return 3.0643+0.01346*p;
 175         }
 176     } else {
 177 //Pb
 178         if (p<15) {
 179             return 2.1380+0.0351*p-0.000853*p*p;
 180         } else {
 181             return 2.407+0.00702*p;
 182         }
 183     }
 184 }
 185
 186
 187
 188
 189
 190