MUON/MUONtestabso.C

   1 void MUONtestabso (Int_t evNumber1=0,Int_t evNumber2=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 do some analysis.
   7 //
   8 /////////////////////////////////////////////////////////////////////////
   9
  10 // Dynamically link some shared libs
  11     static Float_t xmuon, ymuon;
  12
  13    if (gClassTable->GetID("AliRun") < 0) {
  14       gROOT->LoadMacro("loadlibs.C");
  15       loadlibs();
  16    }
  17
  18 // Connect the Root Galice file containing Geometry, Kine and Hits
  19
  20     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  21
  22
  23     if (!file) {
  24         printf("\n Creating galice.root \n");
  25         file = new TFile("galice.root");
  26     } else {
  27         printf("\n galice.root found in file list");
  28     }
  29     file->ls();
  30
  31 // Get AliRun object from file or create it if not on file
  32     if (!gAlice) {
  33         gAlice = (AliRun*)(file->Get("gAlice"));
  34         if (gAlice) printf("AliRun object found on file\n");
  35         if (!gAlice) {
  36             printf("\n Create new gAlice object");
  37             gAlice = new AliRun("gAlice","Alice test program");
  38         }
  39     }
  40
  41
  42 //  Create some histograms
  43
  44
  45    TH1F *zv = new TH1F("zv","z vertex" ,140,470.,505.);
  46    TH1F *xv = new TH1F("xv","x vertex" ,140,-70.,70.);
  47    TH1F *yv = new TH1F("yv","y vertex", 140,-70.,70.);
  48
  49    TH1F *ip  = new TH1F("ip","geant part",50,0.,10.);
  50    TH2F *rzv  = new TH2F("rzv","R-z vert",100,502,504.,100,0.,100.);
  51
  52    TH1F *ptUps  = new TH1F("ptUps","pT Upsilon",50,0.,25.);
  53    TH1F *hde  = new TH1F("hde","dE",200,0.,50);
  54    TH1F *hde2  = new TH1F("hde2","dE",100,1.5,5.5);
  55    TH2F *hdevsn  = new TH2F("hdevsn","dE vs N electron",100,0.,15., 20, 0.,20.);
  56
  57    TH1F *ekine  = new TH1F("ekine","E_kin electrons",70,-5,2);
  58    TH1F *etheta = new TH1F("etheta","Theta electrons",90,0,90);
  59    TH1F *edr = new TH1F("edr","Distance to muon",100,0,10);
  60
  61    TH1F *de  = new TH1F("de","correction",100,-1,1);
  62
  63    TH1F *dtheta = new TH1F("dtheta","Delta Theta" ,200,-5.,5.);
  64    AliMUONChamber*  iChamber;
  65 //
  66 //   Loop over events
  67 //
  68    Int_t Nh=0;
  69    Int_t Nh1=0;
  70        Int_t Nel1=0;
  71        Int_t Nel2=0;
  72        Int_t Nel3=0;
  73        Int_t Nel4=0;
  74
  75    for (Int_t nev=0; nev<= evNumber2; nev++) {
  76        //cout << "nev         " << nev <<endl;
  77        Int_t nparticles = gAlice->GetEvent(nev);
  78        //cout << "nparticles  " << nparticles <<endl;
  79        if (nev < evNumber1) continue;
  80        if (nparticles <= 0) return;
  81
  82        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  83
  84
  85        TTree *TH = gAlice->TreeH();
  86        Int_t ntracks = TH->GetEntries();
  87 //
  88 //   Loop over tracks
  89 //
  90
  91        Float_t dE;
  92
  93        for (Int_t track=0; track<ntracks;track++) {
  94            Int_t Nelt=0;
  95            printf("\n nev %d %d %d %d %d\n ", nev, Nel1, Nel2, Nel3, Nel4);
  96            gAlice->ResetHits();
  97            Int_t nbytes += TH->GetEvent(track);
  98            if (MUON)  {
  99                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
 100                    mHit;
 101                    mHit=(AliMUONHit*)MUON->NextHit())
 102                {
 103                    Int_t   nch   = mHit->fChamber;  // chamber number
 104                    Float_t x     = mHit->fX;        // x-pos of hit
 105                    Float_t y     = mHit->fY;        // y-pos
 106                    Float_t z     = mHit->fZ;        // y-pos
 107
 108                    if (nch != 1) continue;
 109
 110                    Int_t ipart = mHit->fParticle;
 111                    TClonesArray *fPartArray = gAlice->Particles();
 112                    TParticle *Part;
 113                    Int_t ftrack = mHit->fTrack;
 114                    Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
 115                    Int_t ipart = Part->GetPdgCode();
 116                    ip->Fill((float)ipart);
 117                    TParticle *Mother;
 118
 119                    Float_t px0=Part->Px();
 120                    Float_t py0=Part->Py();
 121                    Float_t pz0=Part->Pz();
 122                    Float_t thetax0=TMath::ATan2(px0,pz0);
 123                    Float_t thetay0=TMath::ATan2(py0,pz0);
 124
 125                    if (ipart == kMuonPlus || ipart == kMuonMinus) {
 126 //                     Int_t imo = Part->GetFirstMother();
 127 //                     Mother = (TParticle*) fPartArray->UncheckedAt(imo);
 128 //                     Float_t pt = Mother->Pt();
 129 //                     ptUps->Fill(pt, (float) 1);
 130                        Float_t E=Part->Energy();
 131                        Float_t Eloc=mHit->fPTot;
 132                        Float_t corr=Eloc+CorrectP(Eloc,mHit->fTheta);
 133                        printf("\n %f %f %f", E, Eloc, corr);
 134                        de->Fill(E-corr,1.);
 135                        dE = E-Eloc;
 136                        if (dE<50)  hde->Fill(dE, (float) 1);
 137                        if (dE<5.5) hde2->Fill(dE, (float) 1);
 138                        xmuon=mHit->fX;
 139                        ymuon=mHit->fY;
 140                        Float_t thetax=TMath::ATan2(mHit->Px(), mHit->Momentum());
 141                        Float_t thetay=TMath::ATan2(mHit->Py(), mHit->Momentum());
 142                        dtheta->Fill((thetax-thetax0)*1000., 1.);
 143                        dtheta->Fill((thetay-thetay0)*1000., 1.);
 144                    }
 145
 146                    if (ipart == kElectron || ipart == kPositron) {
 147
 148
 149                        Float_t xvert = Part->Vx();       //  vertex
 150                        Float_t yvert = Part->Vy();       //  vertex
 151                        Float_t zvert = Part->Vz();       // z vertex
 152                        if (zvert < 503 && mHit->fTheta<90) {
 153                            Nelt++;
 154                            Float_t px = Part->Px();
 155                            Float_t py = Part->Py();
 156                            Float_t pz = Part->Pz();
 157                            Float_t Ek = Part->Energy()-Part->GetMass();
 158
 159                            Int_t imo = Part->GetFirstMother();
 160                            Mother = (TParticle*) fPartArray->UncheckedAt(imo);
 161                            Int_t imot = Mother->GetPdgCode();
 162                            xv->Fill(xvert);
 163                            yv->Fill(yvert);
 164                            zv->Fill(zvert);
 165                            Float_t rvert=TMath::Sqrt(xvert*xvert+yvert*yvert);
 166                            rzv->Fill(zvert,rvert);
 167                            ekine->Fill(TMath::Log10(Ek),1.);
 168                            etheta->Fill(mHit->fTheta,1.);
 169                            Float_t ex=mHit->fX;
 170                            Float_t ey=mHit->fY;
 171                            dr=TMath::Sqrt((ex-xmuon)*(ex-xmuon)+(ey-ymuon)*(ey-ymuon));
 172                            edr->Fill(dr,1.);
 173
 174                        }
 175                    }
 176
 177                } // hits
 178            } // if MUON
 179            if (Nelt == 1) Nel1++;
 180            if (Nelt == 2) Nel2++;
 181            if (Nelt == 3) Nel3++;
 182            if (Nelt  > 3) Nel4++;
 183            hdevsn->Fill(dE, (float) Nelt, (float) 1);
 184
 185        } // tracks
 186
 187     } // event
 188
 189
 190 //Create a canvas, set the view range, show histograms
 191    TCanvas *c1 = new TCanvas("c1","Vetices from electrons and positrons",400,10,600,700);
 192    c1->Divide(2,2);
 193
 194    c1->cd(1);
 195    ip->SetFillColor(42);
 196    ip->SetXTitle("ipart");
 197    ip->Draw();
 198
 199    c1->cd(2);
 200    xv->SetFillColor(42);
 201    xv->SetXTitle("xvert");
 202    xv->Draw();
 203
 204    c1->cd(3);
 205    yv->SetFillColor(42);
 206    yv->SetXTitle("yvert");
 207    yv->Draw();
 208
 209    c1->cd(4);
 210    zv->SetFillColor(42);
 211    zv->SetXTitle("zvert");
 212    zv->Draw();
 213
 214    TCanvas *c2 = new TCanvas("c2"," ",400,10,600,700);
 215    c2->Divide(2,2);
 216    c2->cd(1);
 217    rzv->SetXTitle("zvert");
 218    rzv->SetYTitle("rvert");
 219    rzv->Draw();
 220
 221    c2->cd(2);
 222    ptUps->SetXTitle("pt");
 223    ptUps->Draw();
 224
 225    c2->cd(3);
 226    hde->SetXTitle("dE");
 227    hde->Draw();
 228
 229
 230    c2->cd(4);
 231    hde2->SetXTitle("dE");
 232    hde2->Draw();
 233    TCanvas *c3 = new TCanvas("c3"," ",400,10,600,700);
 234    c3->Divide(2,2);
 235    c3->cd(1);
 236    ekine->SetXTitle("E_kin");
 237    ekine->Draw();
 238
 239    c3->cd(2);
 240    etheta->SetXTitle("Theta");
 241    etheta->Draw();
 242
 243    c3->cd(3);
 244    edr->SetXTitle("Distance to muon");
 245    edr->Draw();
 246
 247    c3->cd(4);
 248    dtheta->SetXTitle(" ");
 249    dtheta->Draw();
 250
 251 }
 252
 253 Float_t CorrectP(Float_t p, Float_t theta)
 254 {
 255     printf("\n %f%", theta);
 256
 257     if (theta<3.) {
 258 //W
 259         if (p<15) {
 260             return 2.737+0.0494*p-0.001123*p*p;
 261         } else {
 262             return 3.0643+0.01346*p;
 263         }
 264     } else {
 265 //Pb
 266         if (p<15) {
 267             return 2.1380+0.0351*p-0.000853*p*p;
 268         } else {
 269             return 2.407+0.00702*p;
 270         }
 271     }
 272 }
 273