ITS/SPDclusterTest.C

   1 #include "iostream.h"
   2
   3 void SPDclusterTest (Int_t evNumber1=0,Int_t evNumber2=0)
   4 {
   5 /////////////////////////////////////////////////////////////////////////
   6 //   This macro is a small example of a ROOT macro
   7 //   illustrating how to read the output of GALICE
   8 //   and do some analysis.
   9 //
  10 /////////////////////////////////////////////////////////////////////////
  11
  12 // Dynamically link some shared libs
  13
  14    if (gClassTable->GetID("AliRun") < 0) {
  15       gROOT->LoadMacro("loadlibs.C");
  16       loadlibs();
  17    }
  18
  19 // Connect the Root Galice file containing Geometry, Kine and Hits
  20
  21    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  22    if (!file) file = new TFile("galice.root");
  23    file->ls();
  24
  25 // Get AliRun object from file or create it if not on file
  26
  27    if (!gAlice) {
  28       gAlice = (AliRun*)file->Get("gAlice");
  29       if (gAlice) printf("AliRun object found on file\n");
  30       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  31    }
  32
  33           // ------------ Cluster and point analysis histogramms ------------
  34
  35           TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.);
  36           TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.);
  37           TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.);
  38           TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.);
  39           TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.);
  40           TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.);
  41           TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.);
  42           TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.);
  43
  44           TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.);
  45           TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.);
  46           TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.);
  47           TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.);
  48
  49           // -------------- Create ntuples --------------------
  50           //  ntuple structures:
  51
  52           struct {
  53             Int_t lay;
  54             Int_t nx;
  55             Int_t nz;
  56             Int_t hitprim;
  57             Int_t partcode;
  58             Int_t ntrover;
  59             Float_t dx;
  60             Float_t dz;
  61             Float_t pmod;
  62           } ntuple_st;
  63
  64           struct {
  65             Int_t lay;
  66             Int_t lad;
  67             Int_t det;
  68             Int_t nx;
  69             Int_t nz;
  70             Int_t ntrover;
  71             Int_t noverlaps;
  72             Int_t noverprim;
  73             Float_t qcl;
  74             Float_t dx;
  75             Float_t dz;
  76           } ntuple1_st;
  77
  78           struct {
  79             //      Int_t lay;
  80             Int_t nx;
  81             Int_t nz;
  82           } ntuple2_st;
  83
  84           ntuple = new TTree("ntuple","Demo ntuple");
  85           ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
  86           ntuple->Branch("nx",&ntuple_st.nx,"nx/I");
  87           ntuple->Branch("nz",&ntuple_st.nz,"nz/I");
  88           ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
  89           ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
  90           ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I");
  91           ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
  92           ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
  93           ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
  94
  95           ntuple1 = new TTree("ntuple1","Demo ntuple1");
  96           ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
  97           ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
  98           ntuple1->Branch("det",&ntuple1_st.det,"det/I");
  99           ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I");
 100           ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I");
 101           ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F");
 102           ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");
 103           ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
 104           ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
 105           ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
 106           ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
 107
 108
 109           ntuple2 = new TTree("ntuple2","Demo ntuple2");
 110           //      ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
 111           ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I");
 112           ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I");
 113
 114 // ------------------------------------------------------------------------
 115 //
 116 //   Loop over events
 117 //
 118    for (int nev=0; nev<= evNumber2; nev++) {
 119      Int_t nparticles = gAlice->GetEvent(nev);
 120      cout << "nev         " << nev <<endl;
 121      cout << "nparticles  " << nparticles <<endl;
 122      if (nev < evNumber1) continue;
 123      if (nparticles <= 0) return;
 124
 125      TTree *TH = gAlice->TreeH();
 126      Int_t ntracks = TH->GetEntries();
 127      cout<<"ntracks "<<ntracks<<endl;
 128
 129 // Get pointers to Alice detectors and Digits containers
 130    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 131    TClonesArray *Particles = gAlice->Particles();
 132
 133    if (ITS) {
 134      // fill modules with sorted by module hits
 135      Int_t nmodules;
 136      ITS->InitModules(-1,nmodules);
 137      //    ITS->FillModules(nev,-1,evNumber2,nmodules," "," ");
 138      ITS->FillModules(nev,evNumber2,nmodules," "," ");
 139      //get pointer to modules array
 140      TObjArray *ITSmodules = ITS->GetModules();
 141      AliITShit *itsHit;
 142
 143      // get the Tree for clusters
 144      ITS->GetTreeC(nev);
 145      TTree *TC=ITS->TreeC();
 146      Int_t nent=TC->GetEntries();
 147      printf("Found %d entries in the tree (must be one per module per event!)\n",nent);
 148
 149      for (Int_t idettype=0;idettype<3;idettype++) {
 150
 151        TClonesArray *ITSclusters  = ITS->ClustersAddress(idettype);
 152        //printf ("ITSclusters %p \n",ITSclusters);
 153
 154           if (idettype != 0) continue;
 155
 156           // Module loop
 157           for (Int_t mod=0; mod<nent; mod++) {
 158               AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod);
 159
 160               Int_t nhits = itsModule->GetNhits();
 161               //if(nhits) printf("module nhits %d %d\n",mod,nhits);
 162               if(!nhits) continue;
 163
 164               ITS->ResetClusters();
 165               TC->GetEvent(mod);
 166               Int_t nclust = ITSclusters->GetEntries();
 167               if (!nclust) continue;
 168
 169               // cluster/hit loops
 170
 171         for (Int_t clu=0;clu<nclust;clu++) {
 172                 itsclu   = (AliITSRawClusterSPD*)ITSclusters->UncheckedAt(clu);
 173
 174                 Int_t noverlaps = 0;
 175                 Int_t noverprim = 0;
 176
 177                 Int_t clustersizex = itsclu->NclX();
 178                 Int_t clustersizez = itsclu->NclZ();
 179                 //       Int_t xstart = itsclu->XStart();
 180                 //       Int_t xstop = itsclu->XStop();
 181                 Int_t xstart = itsclu->XStartf();
 182                 Int_t xstop = itsclu->XStopf();
 183                 Float_t fxstart = xstart*50;
 184                 Float_t fxstop = (xstop+1)*50;
 185                 Float_t zstart = itsclu->ZStart();
 186                 Float_t zstop = itsclu->ZStop();
 187                 Int_t zend = itsclu->Zend();
 188                 Int_t ntrover = itsclu->NTracks();
 189                 Float_t clusterx = itsclu->X();
 190                 Float_t clusterz = itsclu->Z();
 191                 Float_t clusterQ = itsclu->Q();
 192
 193                 ntuple2_st.nx = clustersizex;
 194                 ntuple2_st.nz = clustersizez;
 195
 196                 ntuple2->Fill();
 197
 198                 Int_t icl = 0;
 199                 Float_t dxprimlast = 10.e+6;
 200                 Float_t dzprimlast = 10.e+6;
 201
 202                 Float_t SPDlength = 83600;
 203                 Float_t SPDwidth = 12800;
 204                 Float_t xhit0 = 1e+5;
 205                 Float_t zhit0 = 1e+5;
 206
 207
 208        for (Int_t hit=0;hit<nhits;hit++) {
 209
 210                   // Find coordinate differences between the hit and cluster positions
 211                   // for the resolution determination.
 212
 213                   itsHit   = (AliITShit*)itsModule->GetHit(hit);
 214
 215                   Int_t hitlayer = itsHit->GetLayer();
 216                   Int_t hitladder= itsHit->GetLadder();
 217                   Int_t hitdet= itsHit->GetDetector();
 218
 219                   Int_t clusterlayer = hitlayer;
 220                   Int_t clusterladder= hitladder;
 221                   Int_t clusterdetector = hitdet;
 222
 223                   Int_t track = itsHit->GetTrack();
 224                   Int_t dray = 0;
 225                   Int_t hitstat = itsHit->GetTrackStatus();
 226
 227
 228                   Float_t zhit = 10000*itsHit->GetZL();
 229                   Float_t xhit = 10000*itsHit->GetXL();
 230
 231                   if(abs(zhit) > SPDlength/2) {
 232                     if(hitstat == 66) zhit0 = 1e+5;
 233                     continue;
 234                   }
 235
 236                   if(abs(xhit) > SPDwidth/2) {
 237                     if(hitstat == 66) xhit0 = 1e+5;
 238                     continue;
 239                   }
 240
 241                   zhit += SPDlength/2;
 242                   xhit += SPDwidth/2;
 243                   Float_t yhit = 10000*itsHit->GetYL();
 244
 245                   if(hitlayer == 1 && hitstat == 66 && yhit > 71) {
 246                     xhit0 = xhit;
 247                     zhit0 = zhit;
 248                   }
 249                   if(hitlayer == 2 && hitstat == 66 && yhit < -71) {
 250                     xhit0 = xhit;
 251                     zhit0 = zhit;
 252                   }
 253
 254                   if(hitstat != 68) continue; // Take only the hit if the last
 255                   // track point went out from
 256                   // the detector.
 257
 258                   if(xhit0 > 9e+4 || zhit0 > 9e+4) continue;
 259
 260
 261                   Float_t xmed = (xhit + xhit0)/2;
 262                   Float_t zmed = (zhit + zhit0)/2;
 263
 264                   Float_t xdif = xmed - clusterx;
 265                   Float_t zdif = zmed - clusterz;
 266
 267
 268         // Consider the hits inside of cluster region only
 269
 270   if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) {
 271
 272         icl = 1;
 273
 274                 //        part = (TParticle *)particles.UncheckedAt(track);
 275                 //        Int_t partcode = part->GetPdgCode();
 276                 //              Int_t primery = gAlice->GetPrimary(track);
 277
 278         Int_t parent = itsHit->GetParticle()->GetFirstMother();
 279         Int_t partcode = itsHit->GetParticle()->GetPdgCode();
 280
 281 //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
 282 //                      310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 283
 284         Float_t pmod = itsHit->GetParticle()->P(); // the momentum at the
 285                                                    // vertex
 286         pmod *= 1.0e+3;
 287
 288         /*
 289         Float_t px = itsHit->GetPXL();  // the momenta at this GEANT point
 290         Float_t py = itsHit->GetPYL();
 291         Float_t pz = itsHit->GetPZL();
 292         */
 293
 294         Int_t hitprim = 0;
 295
 296         if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 297                                                  // at p < 6 MeV/c
 298
 299         if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but
 300                                                  // not for delta ray which
 301                                                  // also went out from the
 302                                                  // detector and returned
 303                                                  // again
 304
 305         if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
 306                                               // particles
 307
 308         if(hitprim > 0) noverprim = noverprim + 1;
 309
 310         if(hitprim > 0) {
 311          dxprimlast = xdif;
 312          dzprimlast = zdif;
 313         }
 314
 315         // fill ntuple
 316
 317          ntuple_st.lay = hitlayer;
 318          ntuple_st.nx = clustersizex;
 319          ntuple_st.nz = clustersizez;
 320          ntuple_st.hitprim = hitprim;
 321          ntuple_st.partcode = partcode;
 322          ntuple_st.ntrover = ntrover;
 323          ntuple_st.dx = xdif;
 324          ntuple_st.dz = zdif;
 325          ntuple_st.pmod = pmod;
 326
 327          ntuple->Fill();
 328
 329
 330
 331       if(hitprim > 0) {   // for primary particles
 332         if(hitlayer == 1) {
 333            Xres1->Fill(xdif);
 334            Zres1->Fill(zdif);
 335         }
 336         if(hitlayer == 2) {
 337            Xres2->Fill(xdif);
 338            Zres2->Fill(zdif);
 339         }
 340       } // primery particles
 341
 342      } // end of cluster region
 343    } // end of hit loop
 344
 345
 346       if(icl == 1) {
 347
 348         // fill ntuple1
 349
 350         //      ntuple1->Fill(clusterlayer,clustersizex,clustersizez,noverlaps,\
 351 noverprim,dx,dz);
 352
 353         if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
 354                                           // delta rays only
 355
 356          ntuple1_st.lay = clusterlayer;
 357          ntuple1_st.lad = clusterladder;
 358          ntuple1_st.det = clusterdetector;
 359          ntuple1_st.nx = clustersizex;
 360          ntuple1_st.nz = clustersizez;
 361          ntuple1_st.qcl = clusterQ;
 362          ntuple1_st.ntrover = ntrover;
 363          ntuple1_st.noverlaps = noverlaps;
 364          ntuple1_st.noverprim = noverprim;
 365          ntuple1_st.dx = dxprimlast;
 366          ntuple1_st.dz = dzprimlast;
 367
 368          ntuple1->Fill();
 369
 370      } // icl = 1
 371     } // cluster loop
 372    } // module loop
 373   } // idettype loop
 374  } // end if ITS
 375 } // event loop
 376
 377
 378    //  Write and Draw Histogramms and ntuples
 379
 380
 381
 382    TFile fhistos("SPD_his.root","RECREATE");
 383
 384    ntuple->Write();
 385    ntuple1->Write();
 386    ntuple2->Write();
 387
 388    Nxpix1->Write();
 389    Nzpix1->Write();
 390    Nxpix2->Write();
 391    Nzpix2->Write();
 392
 393    Xpix1->Write();
 394    Zpix1->Write();
 395    Xpix2->Write();
 396    Zpix2->Write();
 397
 398    Xres1->Write();
 399    Zres1->Write();
 400    Xres2->Write();
 401    Zres2->Write();
 402
 403    fhistos.Close();
 404    cout<<"!!! Histogramms and ntuples were written"<<endl;
 405
 406    TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
 407    c1->Divide(2,2);
 408    c1->cd(1);
 409    gPad->SetFillColor(33);
 410          Xres1->SetFillColor(42);
 411          Xres1->Draw();
 412    c1->cd(2);
 413    gPad->SetFillColor(33);
 414          Zres1->SetFillColor(46);
 415          Zres1->Draw();
 416    c1->cd(3);
 417    gPad->SetFillColor(33);
 418          Xres2->SetFillColor(42);
 419          Xres2->Draw();
 420    c1->cd(4);
 421    gPad->SetFillColor(33);
 422          Zres2->SetFillColor(46);
 423          Zres2->Draw();
 424
 425      cout<<"END  test for clusters and hits "<<endl;
 426
 427 //     file->Close();
 428 }
 429
 430
 431