ITS/SSDrecpointTest.C

   1 void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
   2   //void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=999)
   3 {
   4 /////////////////////////////////////////////////////////////////////////
   5 //   This macro is a small example of a ROOT macro
   6 //   illustrating how to read the output of GALICE
   7 //   and fill some histograms.
   8 //
   9 //     Root > .L anal.C   //this loads the macro in memory
  10 //     Root > anal();     //by default process first event
  11 //     Root > anal(2);    //process third event
  12 //Begin_Html
  13 /*
  14 <img src="gif/anal.gif">
  15 */
  16 //End_Html
  17 /////////////////////////////////////////////////////////////////////////
  18
  19 // Dynamically link some shared libs
  20
  21    if (gClassTable->GetID("AliRun") < 0) {
  22       gROOT->LoadMacro("loadlibs.C");
  23       loadlibs();
  24    } else {
  25       delete gAlice;
  26       gAlice=0;
  27    }
  28
  29 // Connect the Root Galice file containing Geometry, Kine and Hits
  30    //TString *str = new TString("galice.root");
  31    //TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(str->Data());
  32    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  33    //if (!file) file = new TFile(str->Data(),"UPDATE");
  34    if (!file) file = new TFile("galice.root");
  35    file->ls();
  36
  37 // Get AliRun object from file or create it if not on file
  38    if (!gAlice) {
  39      gAlice = (AliRun*)file->Get("gAlice");
  40      if (gAlice) printf("AliRun object found on file\n");
  41      if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  42    }
  43
  44
  45      // -------------- Create ntuples --------------------
  46
  47      //  ntuple structures:
  48
  49    /*
  50           struct {
  51             Int_t lay;
  52             Int_t nxP;
  53             Int_t nxN;
  54             Int_t hitprim;
  55             Int_t partcode;
  56             Int_t ntrover;
  57             Float_t x;
  58             Float_t z;
  59             Float_t dx;
  60             Float_t dz;
  61             Float_t pmod;
  62           } ntuple_st;
  63    */
  64
  65           struct {
  66             Int_t lay;
  67             Int_t lad;
  68             Int_t det;
  69             Int_t nxP;
  70             Int_t nxN;
  71             Int_t noverlaps;
  72             Int_t noverprim;
  73             Int_t ntrover;
  74             Float_t qclP;
  75             Float_t qclN;
  76             Float_t qrec;
  77             Float_t qcut;
  78             Float_t dx;
  79             Float_t dz;
  80           } ntuple1_st;
  81
  82
  83           struct {
  84             Int_t nxP;
  85             Int_t nxN;
  86             Float_t x;
  87             Float_t z;
  88           } ntuple2_st;
  89
  90           /*
  91           ntuple = new TTree("ntuple","Demo ntuple");
  92           ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
  93           ntuple->Branch("nxP",&ntuple_st.nxP,"nxP/I");
  94           ntuple->Branch("nxN",&ntuple_st.nxN,"nxN/I");
  95           ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
  96           ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
  97           ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I");
  98           ntuple->Branch("x",&ntuple_st.x,"x/F");
  99           ntuple->Branch("z",&ntuple_st.z,"z/F");
 100           ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
 101           ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
 102           ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
 103           */
 104
 105           ntuple1 = new TTree("ntuple1","Demo ntuple1");
 106           ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
 107           ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
 108           ntuple1->Branch("det",&ntuple1_st.det,"det/I");
 109           ntuple1->Branch("nxP",&ntuple1_st.nxP,"nxP/I");
 110           ntuple1->Branch("nxN",&ntuple1_st.nxN,"nxN/I");
 111           ntuple1->Branch("qclP",&ntuple1_st.qclP,"qclP/F");
 112           ntuple1->Branch("qclN",&ntuple1_st.qclN,"qclN/F");
 113           ntuple1->Branch("qrec",&ntuple1_st.qrec,"qrec/F");
 114           ntuple1->Branch("qcut",&ntuple1_st.qcut,"qcut/F");
 115           ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
 116           ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
 117           ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
 118           ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
 119           ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");
 120
 121
 122           ntuple2 = new TTree("ntuple2","Demo ntuple2");
 123           ntuple2->Branch("nxP",&ntuple2_st.nxP,"nxP/I");
 124           ntuple2->Branch("nxN",&ntuple2_st.nxN,"nxN/I");
 125           ntuple2->Branch("x",&ntuple2_st.x,"x/F");
 126           ntuple2->Branch("z",&ntuple2_st.z,"z/F");
 127
 128
 129           // Create Histogramms
 130
 131           TH1F *NxP5 = new TH1F("NxP5","P cluster size for layer 5",20,0.,20.);
 132           TH1F *NxN5 = new TH1F("NxN5","N cluster size for layer 5",20,0.,20.);
 133           TH1F *NxP6 = new TH1F("NxP6","P cluster size for layer 6",20,0.,20.);
 134           TH1F *NxN6 = new TH1F("NxN6","N cluster size for layer 6",20,0.,20.);
 135
 136           TH1F *Xres5 = new TH1F("Xres5","Xrec and Xgen difference (micr) for layers 5",100,-200.,200.);
 137           TH1F *Xres6 = new TH1F("Xres6","Xrec and Xgen difference (micr) for layers 6",100,-200.,200.);
 138           TH1F *Zres5 = new TH1F("Zres5","Zrec and Zgen difference (micr) for layers 5",100,-8000.,8000.);
 139           TH1F *Zres6 = new TH1F("Zres6","Zrec and Zgen difference (micr) for layers 6",100,-8000.,8000.);
 140           TH1F *Path5 = new TH1F("Path5","Path length in Si",100,0.,600.);
 141           TH1F *Path6 = new TH1F("Path6","Path length in Si",100,0.,600.);
 142           TH1F *dEdX = new TH1F("dEdX","dEdX  (KeV)",100,0.,500.);
 143           TH2F *adcPadcN5all = new TH2F("adcPadcN5all","adcP/N correlation for lay5",100,0.,200.,100,0.,200.);
 144           TH2F *adcPadcN6all = new TH2F("adcPadcN6all","adcP/N correlation for lay6",100,0.,200.,100,0.,200.);
 145           TH2F *adcPadcN5cut = new TH2F("adcPadcN5cut","adcP/N correlation for lay5 and cut of P-N signas",100,0.,200.,100,0.,200.);
 146           TH2F *adcPadcN6cut = new TH2F("adcPadcN6cut","adcP/N correlation for lay6 and cut of P-N signals",100,0.,200.,100,0.,200.);
 147           //-----------------------------------------------------------
 148
 149 //
 150 //   Loop over events
 151 //
 152   for (int nev=0; nev<= evNumber2; nev++) {
 153      Int_t nparticles = gAlice->GetEvent(nev);
 154      cout << "nev         " << nev <<endl;
 155      //cout << "nparticles  " << nparticles <<endl;
 156      if (nev < evNumber1) continue;
 157      if (nparticles <= 0) return;
 158
 159      TTree *TH = gAlice->TreeH();
 160      Int_t ntracks = TH->GetEntries();
 161      cout<<"all entries to GEANT(charged and neutral) "<<ntracks<<endl;
 162
 163 // Get pointers to Alice detectors and Digits containers
 164      AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 165      TClonesArray *Particles = gAlice->Particles();
 166
 167    if (ITS) {
 168
 169      // fill modules with sorted by module hits
 170      Int_t nmodules;
 171      ITS->InitModules(-1,nmodules);
 172      ITS->FillModules(nev,evNumber2,nmodules," "," ");
 173
 174      //get pointer to modules array
 175      TObjArray *ITSmodules = ITS->GetModules();
 176      AliITShit *itsHit;
 177
 178      // get the Tree for clusters
 179      ITS->GetTreeC(nev);
 180      TTree *TC=ITS->TreeC();
 181      Int_t nent=TC->GetEntries();
 182      printf("Found %d entries in the TreeC (full number of the modules)\n",nent);
 183      TTree *TR = gAlice->TreeR();
 184      Int_t lay, lad, det;
 185      AliITSgeom *geom = ITS->GetITSgeom();
 186      Int_t mod;
 187
 188      Int_t first0 = geom->GetStartDet(0);  // SPD
 189      Int_t last0 = geom->GetLastDet(0);    // SPD
 190      Int_t first1 = geom->GetStartDet(1);  // SDD
 191      Int_t last1 = geom->GetLastDet(1);    // SDD
 192      Int_t first2 = geom->GetStartDet(2);  // SSD
 193      Int_t last2 = geom->GetLastDet(2);    // SSD
 194
 195      //  For the SPD: first0 = 0, last0 = 239     (240 modules);
 196      //  for the SDD: first1 = 240, last1 = 499   (260 modules);
 197      //  for the SSD: first2 = 500, last2 = 2197  (1698 modules).
 198
 199      printf("det type %d first0, last0 %d %d \n",0,first0,last0);
 200      printf("det type %d first1, last1 %d %d \n",1,first1,last1);
 201      printf("det type %d first2, last2 %d %d \n",2,first2,last2);
 202
 203
 204    AliITSDetType *iDetType=ITS->DetType(2);
 205    AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
 206    //AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
 207
 208    printf("SSD dimensions %f %f %f \n",seg2->Dx(),seg2->Dz(),seg2->Dy());
 209    printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
 210    Float_t ylim = seg2->Dy()/2 - 12;
 211
 212
 213
 214     for (Int_t idettype=0;idettype<3;idettype++) {
 215
 216        TClonesArray *ITSclu  = ITS->ClustersAddress(idettype);
 217        TClonesArray *ITSrec  = ITS->RecPoints();
 218        //printf ("ITSrec %p \n",ITSrec);
 219        //printf ("ITSclu %p \n",ITSclu);
 220
 221           if (idettype != 2) continue;
 222
 223           // Module loop
 224
 225      for (mod=first2; mod<last2+1; mod++) {  // for the "ALL" option
 226        //for (mod=0; mod<last2-first2+1; mod++) { //for the "SSD" option
 227
 228        //AliITSmodule *Mod = (AliITSmodule *)ITSmodules->At(mod+first2);
 229        // for the "SSD" option
 230
 231        AliITSmodule *Mod = (AliITSmodule *)ITSmodules->At(mod);
 232        // for the "ALL" option
 233
 234               geom->GetModuleId(mod,lay,lad,det);
 235
 236               Int_t nhits = Mod->GetNhits();
 237               //if(nhits) printf("module nhits %d %d\n",mod,nhits);
 238               if(!nhits) continue;
 239
 240               ITS->ResetClusters();
 241               TC->GetEvent(mod);
 242               Int_t nclust = ITSclu->GetEntries();
 243               if (!nclust) continue;
 244
 245               ITS->ResetRecPoints();
 246               TR->GetEvent(mod);
 247               Int_t nrecp = ITSrec->GetEntries();
 248               //if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
 249               if (!nrecp) continue;
 250
 251
 252
 253        Float_t epart = 0;
 254        cout <<" module,nrecp,nclust,nhits ="<<mod<<","<<nrecp<<","<<nclust<<","<<nhits<< endl;
 255
 256        // ---------------- cluster/recpoint/hit analysis ---------------------
 257
 258
 259        Float_t pathInSSD = 300.;
 260
 261        // ---- Recpoint loop
 262       for (Int_t pnt=0;pnt<nrecp;pnt++) {
 263
 264          itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
 265          if(!itsPnt) continue;
 266          itsClu  = (AliITSRawClusterSSD*)ITSclu->At(pnt);
 267          if(!itsClu) continue;
 268
 269          Int_t nxP = itsClu->fMultiplicity;
 270          Int_t nxN = itsClu->fMultiplicityN;
 271          Int_t ntrover = itsClu->fNtracks;
 272          Float_t qclP = itsClu->fSignalP;     // in ADC
 273          Float_t qclN = itsClu->fSignalN;     // in ADC
 274          Float_t dq = TMath::Abs(qclP - qclN);
 275          Float_t xrec = 10000*itsPnt->GetX();
 276          Float_t zrec = 10000*itsPnt->GetZ();
 277          Float_t qrec = itsPnt->GetQ();      // in ADC, maximum from fSignalP/N
 278          Float_t qcut = dq/qrec;
 279          //Float_t dedx = itsPnt->GetdEdX();   // in KeV (ADC * 2.16)
 280          Float_t dedx = itsPnt->fdEdX;   // in KeV (ADC * 2.16)
 281          Int_t ii = 0;
 282          Int_t tr1 = itsPnt->GetLabel(ii);
 283          Int_t ii = 1;
 284          Int_t tr2 = itsPnt->GetLabel(ii);
 285          Int_t ii = 2;
 286          Int_t tr3 = itsPnt->GetLabel(ii);
 287
 288          // fill ntuple2
 289          ntuple2_st.nxP = nxP;
 290          ntuple2_st.nxN = nxN;
 291          ntuple2_st.x = xrec/1000;
 292          ntuple2_st.z = zrec/1000;
 293
 294          ntuple2->Fill();
 295
 296          Int_t noverlaps = 0;
 297          Int_t noverprim = 0;
 298          Int_t flaghit = 0;
 299          Float_t xhit0 = 1e+7;
 300          Float_t yhit0 = 1e+7;
 301          Float_t zhit0 = 1e+7;
 302          Float_t dxprimbest = 1e+7;
 303          Float_t dzprimbest = 1e+7;
 304
 305          // Hit loop
 306        for (Int_t hit=0;hit<nhits;hit++) {
 307
 308            itsHit   = (AliITShit*)Mod->GetHit(hit);
 309
 310            Int_t flagtrack = 0;
 311            Int_t hitlayer = itsHit->GetLayer();
 312            Int_t hitladder= itsHit->GetLadder();
 313            Int_t hitdet= itsHit->GetDetector();
 314
 315            Int_t track = itsHit->GetTrack();
 316            Int_t dray = 0;
 317            Int_t hitstat = itsHit->GetTrackStatus();
 318
 319            Float_t zhit = 10000*itsHit->GetZL();
 320            Float_t xhit = 10000*itsHit->GetXL();
 321            Float_t yhit = 10000*itsHit->GetYL();
 322            Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
 323
 324            Int_t parent = itsHit->GetParticle()->GetFirstMother();
 325            Int_t partcode = itsHit->GetParticle()->GetPdgCode();
 326
 327    //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
 328    //  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 329
 330            Float_t pmod = itsHit->GetParticle()->P(); // the momentum at the
 331                                                       // vertex
 332            pmod *= 1.0e+3;
 333
 334           if(hitstat == 66 && yhit < -ylim) {
 335             xhit0 = xhit;
 336             yhit0 = yhit;
 337             zhit0 = zhit;
 338           }
 339
 340           if(hitstat == 66) continue; // Take the not entering hits only
 341
 342           if(xhit0 > 9e+6 || zhit0 > 9e+6 || yhit0 > 9e+6) {
 343             continue;
 344           }
 345
 346           // Consider the hits only with the track number equaled to one
 347           // of the recpoint
 348           if((track == tr1) || (track == tr2) || (track == tr3)) flagtrack = 1;
 349
 350           if(flagtrack == 1) {     // the hit corresponds to the recpoint
 351
 352             flaghit = 1;
 353
 354             //Float_t px = itsHit->GetPXL(); // the momenta at this GEANT point
 355             //Float_t py = itsHit->GetPYL();
 356             //Float_t pz = itsHit->GetPZL();
 357
 358             Int_t hitprim = 0;
 359
 360             if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 361             // at p < 6 MeV/c
 362
 363             if((hitstat == 68 || hitstat == 33) && dray == 0)  noverlaps=noverlaps + 1;
 364             // overlapps for all hits but
 365             // not for delta ray which
 366             // also went out from the
 367             // detector and returned
 368             // again
 369
 370
 371             // x,z resolution colculation
 372             if((hitstat == 68 || hitstat == 33) && dray == 0) {
 373               Float_t xmed = (xhit + xhit0)/2;
 374               Float_t zmed = (zhit + zhit0)/2;
 375               Float_t xdif = xmed - xrec;
 376               Float_t zdif = zmed - zrec;
 377
 378               if(parent < 0)  {
 379                 hitprim = 1; // hitprim=1 for the primery particles
 380                 noverprim += 1;
 381               }
 382               pathInSSD = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
 383
 384               // Find the best xdif and zdif from any ones for the primery
 385               // particles (to remove the wronge xdif and zdif if the hit
 386               // belongs to the other package containing the same P/N cluster
 387
 388               if(hitprim > 0) {
 389                 if(TMath::Abs(dxprimbest)>TMath::Abs(xdif)) dxprimbest = xdif;
 390                 if(TMath::Abs(dzprimbest)>TMath::Abs(zdif)) dzprimbest = zdif;
 391               }
 392
 393
 394               /*
 395               // fill ntuple
 396               ntuple_st.lay = hitlayer;
 397               ntuple_st.nxP = nxP;
 398               ntuple_st.nxN = nxN;
 399               ntuple_st.hitprim = hitprim;
 400               ntuple_st.partcode = partcode;
 401               ntuple_st.ntrover = ntrover;
 402               ntuple_st.x = xrec/1000;
 403               ntuple_st.z = zrec/1000;
 404               ntuple_st.dx = xdif;
 405               ntuple_st.dz = zdif;
 406               ntuple_st.pmod = pmod;
 407
 408               //if(qcut < 0.18) ntuple->Fill();
 409               ntuple->Fill();
 410               */
 411
 412               //if(hitlayer == 5 && qcut < 0.18) {
 413
 414               if(hitlayer == 5 ) {
 415                 Xres5->Fill(xdif);
 416                 Zres5->Fill(zdif);
 417                 Path5->Fill(pathInSSD);
 418               }
 419               //if(hitlayer == 6 && qcut < 0.18) {
 420               if(hitlayer == 6) {
 421                 Xres6->Fill(xdif);
 422                 Zres6->Fill(zdif);
 423                 Path6->Fill(pathInSSD);
 424               }
 425             } // hitstat 68/33
 426           } else {       // non correspondent hit
 427             xhit0 = 1e+7;
 428             zhit0 = 1e+7;
 429           } // end of hit-recpoint correspondence
 430        } // hit loop
 431
 432        if(flaghit == 1) {
 433
 434            if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
 435            // delta rays only
 436
 437
 438            // fill ntuple1
 439            ntuple1_st.lay = hitlayer;
 440            ntuple1_st.lad = hitladder;
 441            ntuple1_st.det = hitdet;
 442            ntuple1_st.nxP = nxP;
 443            ntuple1_st.nxN = nxN;
 444            ntuple1_st.qclP = qclP*300/pathInSSD;
 445            ntuple1_st.qclN = qclN*300/pathInSSD;
 446            ntuple1_st.qrec = qrec*300/pathInSSD;
 447            ntuple1_st.qcut = qcut;
 448            ntuple1_st.dx = dxprimbest;
 449            ntuple1_st.dz = dzprimbest;
 450            noverlaps -= 1;
 451            noverprim -= 1;
 452            ntuple1_st.noverlaps = noverlaps;
 453            ntuple1_st.noverprim = noverprim;
 454            ntuple1_st.ntrover = ntrover;
 455
 456            //if(qcut < 0.18) ntuple1->Fill();
 457            ntuple1->Fill();
 458
 459
 460            Float_t de = dedx*300./pathInSSD;
 461            dEdX->Fill(de);
 462            if(noverprim >=0) {
 463              if(hitlayer == 5 ) {
 464                adcPadcN5all->Fill(qclP,qclN);
 465              }
 466              if(hitlayer == 6 ) {
 467                adcPadcN6all->Fill(qclP,qclN);
 468              }
 469              if(hitlayer == 5 && qcut < 0.18) {
 470                //if(hitlayer == 5 && noverlaps == 0) {
 471                adcPadcN5cut->Fill(qclP,qclN);
 472                NxP5->Fill(nxP);
 473                NxN5->Fill(nxN);
 474              }
 475              if(hitlayer == 6 && qcut < 0.18) {
 476                //if(hitlayer == 6 && noverlaps == 0) {
 477                adcPadcN6cut->Fill(qclP,qclN);
 478                NxP6->Fill(nxP);
 479                NxN6->Fill(nxN);
 480              }
 481            }
 482        } // flaghit = 1
 483       } //b.b. recpoint loop
 484      } // module loop
 485     } // detector loop (iDetType)
 486    } // if ITS
 487   } //b.b. evnt loop
 488
 489    TFile fhistos("SSD_his.root","RECREATE");
 490
 491    //ntuple->Write();
 492    ntuple1->Write();
 493    ntuple2->Write();
 494    NxP5->Write();
 495    NxN5->Write();
 496    NxP6->Write();
 497    NxN6->Write();
 498    Xres5->Write();
 499    Zres5->Write();
 500    Xres6->Write();
 501    Zres6->Write();
 502    Path5->Write();
 503    Path6->Write();
 504    adcPadcN5all->Write();
 505    adcPadcN6all->Write();
 506    adcPadcN5cut->Write();
 507    adcPadcN6cut->Write();
 508    dEdX->Write();
 509
 510    fhistos.Close();
 511
 512    cout<<"!!! Histogramms and ntuples were written"<<endl;
 513
 514    TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
 515    c1->Divide(2,2);
 516    c1->cd(1);
 517    gPad->SetFillColor(33);
 518          Xres5->SetFillColor(42);
 519          Xres5->Draw();
 520    c1->cd(2);
 521    gPad->SetFillColor(33);
 522          Zres5->SetFillColor(46);
 523          Zres5->Draw();
 524    c1->cd(3);
 525    gPad->SetFillColor(33);
 526          Xres6->SetFillColor(42);
 527          Xres6->Draw();
 528    c1->cd(4);
 529    gPad->SetFillColor(33);
 530          Zres6->SetFillColor(46);
 531          Zres6->Draw();
 532
 533    cout<<"END  test for clusters and hits "<<endl;
 534
 535 }
 536
 537
 538
 539
 540
 541