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