Removing obsolete macros

[u/mrichter/AliRoot.git] / ITS / AliITSReadPlotData.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include<iostream.h>
#include<TROOT.h>
#include<TArrayI.h>
#include<TBranch.h>
#include<TCanvas.h>
#include<TClassTable.h>
#include<TClonesArray.h>
#include<TFile.h>
#include<TH1.h>
#include<TH2.h>
#include<TLegend.h>
#include<TObject.h>
#include<TObjArray.h>
#include<TTree.h>
#include <AliRun.h>
#include <AliITS.h>
#include <AliITSgeom.h>
#include <AliITSDetType.h>
#include <AliITSRecPoint.h>
#include <AliITSdigit.h>
#include <AliITShit.h>
#include <AliITSmodule.h> 
#include <AliITSsegmentation.h>
#include <AliITSsegmentationSPD.h> 
#include <AliITSsegmentationSDD.h>
#include <AliITSsegmentationSSD.h>
#endif

Int_t GetModuleHits (TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St);
Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z, Int_t isfastpoints);
Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z);

Int_t AliITSReadPlotData( Int_t evNum = 0, char *filename = "galice.root", TString FileDigits="galice.root", TString FileRec="galice.root", Int_t isfastpoints = 0) {

  /*********************************************************************
   *                                                                   *
   *  Macro used to read hits, digits and recpoints for a module       *
   *  It draws 3 TH2Fs where stores the 2-D coordinates of these       *
   *  data for a specified module (for which the macro asks when       *
   *  it starts, with some checks to avoid wrong detector coords.      *
   *                                                                   *
   *  Only a little 'experimental' warning:                            *
   *  with all the tests I have made, I wasn't able to plot the        *
   *  digits fo the 5th layer...                                       *
   *  I skipped this problem with an alert to beware th user, while    *
   *  in this situation the macro plots only recpoints and hits        *
   *                                                                   *
   *  Author: Alberto Pulvirenti                                       *
   *                                                                   *
   *********************************************************************/

	
  // First of all, here are put some variable declarations
  // that are useful in the following part:
  Int_t nparticles; // number of particles
  // ITS module coordinates [layer = 1, ladder = 2, det = 3] and absolute ID[0] of module [0]
  TArrayI ID(4);
  Int_t nmodules, dtype; // Total number of modules and module type (SSD, SPD, SDD)
  Float_t *x = 0, *y = 0, *z = 0; // Arrays where to store read coords
  Bool_t *St = 0; // Status of the track (hits only)
  /*
	// It's necessary to load the gAlice shared libs
	// if they aren't already stored in memory...
	if (gClassTable->GetID("AliRun") < 0) {
   	gROOT->LoadMacro("loadlibs.C");
    loadlibs();
	}
    // Anyway, this macro needs to read a gAlice file, so it
    // clears the gAlice object if there is already one in memory...
    else {
  */
  if(gAlice){
    delete gAlice;
    gAlice = 0;
  }
  //	}

  // Now is opened the Root input file containing Geometry, Kine and Hits
  // by default its name must be "galice.root".
  // When the file is opened, its contens are shown.
  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
  if (!file) file = new TFile(filename);
  file->ls();
	
  // Then, the macro gets the AliRun object from file.
  // If this object is not present, an error occurs
  // and the execution is stopped.
  // Anyway, this operation needs some time,
  // don't worry about an apparent absence of output and actions...
  cout << "\nSearching in '" << filename << "' for an AliRun object ... " << flush;
  gAlice = (AliRun*)file->Get("gAlice");
  if (gAlice)
    cout << "FOUND!" << endl;
  else {
    cout<<"NOT FOUND! The Macro can't continue!" << endl;
    return 0;
  }
  if(!(FileDigits.Data() == filename)){
    gAlice->SetTreeDFileName(FileDigits);
  }
  if(!(FileRec.Data() == filename)){
    gAlice->SetTreeRFileName(FileRec);
  }
	
  // Then, the macro selects the event number specified. Default is 0.
  nparticles = gAlice->GetEvent(evNum);
  cout << "\nNumber of particles   = " << nparticles << endl;
  if (!nparticles) {
    cout << "With no particles I can't do much... Goodbye!" << endl;
    return 0;
  }
	
  // The next step is filling the ITS from the AliRun object.
  AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
  ITS->InitModules(-1, nmodules);
  cout << "Number of ITS modules = " << nmodules << endl;
  cout << "\nFilling modules (it takes a while, now)..." << flush;
  ITS->FillModules(0, 0, nmodules, " ", " ");
  cout << "DONE!" << endl;
  AliITSgeom *gm = ITS->GetITSgeom();
  //	AliITSDetType *det = ITS->DetType(dtype);	
  //	AliITSsegmentation *seg = det->GetSegmentationModel();	
	
  for(;;) {

    // Input phase.
    // The macro asks if the user wants to put a detector ID[0]
    // or prefers to input layer, ladder and detector.
    for (Int_t i = 0; i < 4; i++) ID[i] = -1;
    Int_t answ;
    do {
      cout << "\nSelection modes:" << endl;
      cout << "1) by layer - ladder - detector numbers" << endl;
      cout << "2) by unique detector ID" << endl;
      cout << "0) exit macro" << endl;
      cout << "\nEnter your choice: ";
      cin >> answ;
    } while (answ < 0 || answ > 2);
    switch (answ) {
    case 0:
      // input = 0 ---> EXIT
      return 0;
      break;
    case 1:
      // input = 1 ---> SELECTION BY COORD
      do {
        cout << "\nLayer number [1-6, 0 to exit]: ";
        cin >> ID[1];
        if (!ID[1]) return 0;
      } while (ID[1] < 0 || ID[1] > 6);
				
      // Detector type: 0 --> SPD, 1 --> SDD, 2 --> SSD.
      // Layer 1,2 --> 0 / Layer 3,4 --> 1 / Layer 5,6 --> 2
      dtype = (ID[1] - 1) / 2;
				
      // Once fixed the layer number, the macro calculates the max number
      // for ladder and detector from geometry, and accepts only suitable values.
      do {
        ID[2] = gm->GetNladders(ID[1]);
        cout << "Ladder number [1-" << ID[2] << ", 0 to exit]: ";
        cin >> ID[2];
        if (!ID[2]) return 0;
      } while (ID[2] < 0 || ID[2] > gm->GetNladders(ID[1]));
      do {
        ID[3] = gm->GetNdetectors(ID[1]);
        cout << "Detector number [1-" << ID[3] << ", 0 to exit]: ";
        cin >> ID[3];
        if (!ID[3]) return 0;
      } while (ID[3] < 0 || ID[3] > gm->GetNdetectors(ID[1]));
      break;
    case 2:
      // input = 2 ---> SELECTION BY ID[0]
      do {
        ID[0] = gm->GetIndexMax();
        cout << "\n Detector ID number [0-" << ID[0]-1 << ", -1 to exit]: ";
        cin >> ID[0];
        ID[0]++;
        if (ID[0] == -1) return 0;
      } while (ID[0] < 0 || ID[0] > gm->GetIndexMax());
      break;
    };

    if (ID[0] == -1)
      // If ID[0] = -1 the chioce was by coords, so it's necessary to assign the ID:
      ID[0] = gm->GetModuleIndex(ID[1], ID[2], ID[3]);
    else {	
      // Else we must get the layer, ladder and detector by the ID;
      // ...but first we must remember that the ID goest from 0 to NModules - 1!
      ID[0]--;
      ID[1] = ITS->GetModule(ID[0])->GetLayer();
      ID[2] = ITS->GetModule(ID[0])->GetLadder();
      ID[3] = ITS->GetModule(ID[0])->GetDet();
    }
				
    // Defines the histograms inside the `for' cycle, so they are destroyed at the end
    // of every read sequqnce, in order to mek another withour segmentation faults
    Text_t msg[250];
    Float_t xm = 0.0, zm = 0.0;
    dtype = (ID[1] - 1) / 2;
    switch (dtype) {
    case 0: xm = 1.5; zm = 7.0; break;
    case 1: xm = 7.5; zm = 8.0; break;
    case 2: xm = 7.5; zm = 4.5; break;
    }
    sprintf(msg, "Module index=%d lay=%d lad=%d det=%d", ID[0], ID[1], ID[2], ID[3]);
    TH2F *hhits = new TH2F("hhits", msg, 500, -xm, xm, 500, -zm, zm);     // Histogram of hits
    TH2F *hrecs = new TH2F("hrecs", msg, 500, -xm, xm, 500, -zm, zm);     // Histogram of recpoints
    TH2F *hdigits = new TH2F("hdigits", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of digits
		
    cout << endl;
		
    // Reads hits...
    Int_t tmpint = ID[0];
    Int_t hits = GetModuleHits(ITS, tmpint, x, y, z, St);
    if (!hits) {
      cout << "No hits in module!" << endl;
      continue;
    }
    else
      cout << "Hits scanned..." << endl;
    for (Int_t i = 0; i < hits; i++) if (!St[i]) hhits->Fill(x[i], z[i]);
		
    // Reads recpoints...
    Int_t recs = GetModuleRecPoints(ITS, ID[0], x, z, isfastpoints);
    if (!recs) {
      cout << "No recpoints in module!" << endl;
      continue;
    }
    else
      cout << "Recpoints scanned..." << endl;
    for (Int_t i = 0; i < recs; i++) hrecs->Fill(x[i], z[i]);
		
    // Reads digits...
    Int_t digits = GetModuleDigits(ITS, ID[0], dtype, x, z);
    if (!digits) {
      cout << "No digits in module!" << endl;
      //continue;
    }
    else
      cout << "Digits scanned..." << endl;
    for (Int_t i = 0; i < digits; i++) hdigits->Fill(x[i], z[i]);

    // Draws read data...
    // HITS -------> red (2) crosses.
    // DIGITS -----> green (8) boxes.
    // REC-POINTS -> blue (4) St. Andrew's crosses.

    TCanvas *viewer = new TCanvas("viewer", "Module viewer canvas", 0, 0, 800, 800);
    viewer->cd();
		
    hdigits->SetMarkerStyle(25);
    hdigits->SetMarkerColor(8);
    hdigits->SetMarkerSize(2);
    hdigits->SetStats(kFALSE);
    hdigits->SetXTitle("Local X (cm)");
    hdigits->SetYTitle("Local Z (cm)");
    hdigits->Draw();
		
    hhits->SetMarkerStyle(5);
    hhits->SetMarkerColor(2);
    hhits->SetMarkerSize(3);
    hhits->SetStats(kFALSE);
    hhits->Draw("same");
	
    hrecs->SetMarkerStyle(2);
    hrecs->SetMarkerColor(4);
    hrecs->SetMarkerSize(3);
    hrecs->SetStats(kFALSE);
    hrecs->Draw("same");
		
    TLegend *legend = new TLegend(0.7, 0.8, 0.9, 0.9);
    legend->SetMargin(0.2);
    legend->AddEntry(hhits, "hits","P");
    legend->AddEntry(hrecs, "recpoints","P");
    legend->AddEntry(hdigits, "digits","P");
    legend->SetTextSizePixels(14);
    legend->Draw();
		
    viewer->Update();
		
		
    Text_t fname[250],ans;
    cout << "Do you want to save the current canvas on a file (y/n) ? ";
    cin >> ans;
    if(ans == 'y' || ans == 'Y') {
      cout << "Enter filename: ";
      cin >> fname;
      TString *control = new TString(fname);
      Bool_t ok=control->Contains(".C") || control->Contains(".root") || control->Contains(".ps") || control->Contains(".eps") || control->Contains(".gif");
      if(!ok){ 
        cout << "File extension is not recognized. The canvas will be saved as Postscript file";
        strcat(fname,".ps");
      }  
      viewer->SaveAs(fname);
    }
  }
	
  cout << "Done. Goodbye" << endl;
  return 0;
}

Int_t GetModuleHits(TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St) {	
  // First of all, the macro selects the specified module,
  // then gets the array of hits in it and their number.
  AliITS *ITS = (AliITS*) its;
  AliITSmodule *module = ITS->GetModule(ID);
  TObjArray *hits_array = module->GetHits();
  Int_t hits_num = hits_array->GetEntriesFast();
	
	// Now, if this count returns 0, there's nothing to do,
	// while, if it doesn't, the first thing to do is dimensioning
	// the coordinate arrays, and then the loop can start.
  if (!hits_num)
    return 0;
  else {
    if (X) delete [] X;	
    if (Y) delete [] Y;
    if (Z) delete [] Z;
    if (St) delete [] St;
    X = new Float_t[hits_num];
    Y = new Float_t[hits_num];
    Z = new Float_t[hits_num];
    St = new Bool_t[hits_num];
  }

  for (Int_t j = 0; j < hits_num; j++) {
    AliITShit *hit = (AliITShit*) hits_array->At(j);
    X[j]  = hit->GetXL();
    Y[j]  = hit->GetYL();
    Z[j]  = hit->GetZL();
    St[j] = hit->StatusEntering();
  }
  return hits_num;
}

Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z, Int_t isfastpoints) {
	
  // First of all, the macro selects the specified module,
  // then gets the array of recpoints in it and their number.
  AliITS *ITS = (AliITS*) its;
  TTree *TR = gAlice->TreeR();
  ITS->ResetRecPoints();
	
  TClonesArray *recs_array = ITS->RecPoints();

  TBranch *branch = 0;
  if(TR && recs_array){
    if(isfastpoints==1){
      branch = gAlice->TreeR()->GetBranch("ITSRecPointsF");
      cout<<"using fast points\n";
    }
    else {
      branch = gAlice->TreeR()->GetBranch("ITSRecPoints");
    }
    if(branch)branch->SetAddress(&recs_array);
  }
        
  branch->GetEvent(ID);

  Int_t recs_num = recs_array->GetEntries();


  // Now, if this count returns 0, there's nothing to do,
  // while, if it doesn't, the first thing to do is dimensioning
  // the coordinate and energy loss arrays, and then the loop can start.
  if (!recs_num)
    return 0;
  else {
    if (X) delete [] X;	
    if (Z) delete [] Z;
    X = new Float_t[recs_num];
    Z = new Float_t[recs_num];
  }
  for(Int_t j = 0; j < recs_num; j++) {
    AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j);
    X[j] = recp->GetX();
    Z[j] = recp->GetZ();
  }
  return recs_num;	
}

Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z) {

  // First of all, the macro selects the specified module,
  // then gets the array of recpoints in it and their number.
  AliITSdigit *digit;
  AliITS *ITS = (AliITS*)its;
  TTree *TD = gAlice->TreeD();
  ITS->ResetDigits();
  TD->GetEvent(ID);
  TClonesArray *digits_array = ITS->DigitsAddress(dtype);
  AliITSgeom *gm = ITS->GetITSgeom();	
  AliITSDetType *det = ITS->DetType(dtype);	
  AliITSsegmentation *seg = det->GetSegmentationModel();	
  TArrayI ssdone(5000);  // used to match p and n side digits of strips
  TArrayI pair(5000);    // as above 	
  Int_t digits_num = digits_array->GetEntries();
  // Now, if this count returns 0, there's nothing to do,
  // while, if it doesn't, the first thing to do is dimensioning
  // the coordinate and energy loss arrays, and then the loop can start.

  if(dtype==2){
    Int_t layer, ladder, detec;
    gm->GetModuleId(ID,layer,ladder,detec);
    seg->SetLayer(layer);
  }

  if (!digits_num)
    return 0;
  else {
    cout << "Digits to scan: " << digits_num << endl;
    if (X) delete [] X;			
    if (Z) delete [] Z;
    X = new Float_t[digits_num];		
    Z = new Float_t[digits_num];
  }
	
  // Get the coordinates of the module
  if (dtype == 2) {
    for (Int_t j = 0; j < digits_num; j++) {
      ssdone[j] = 0;			
      pair[j] = 0;
    }
  }
  //    AliITSdigit *digit;
  for (Int_t j = 0; j < digits_num; j++) {
    digit = (AliITSdigit*)digits_array->UncheckedAt(j);
    Int_t iz=digit->fCoord1;  // cell number z
    Int_t ix=digit->fCoord2;  // cell number x
    // Get local coordinates of the element (microns)
    // ******************************* PARTE CORRETTA ***************************************
    if(dtype < 2) {
      Float_t xx, zz; // aggiunta
      seg->DetToLocal(ix, iz, xx, zz);
      X[j] = xx; // aggiunta
      Z[j] = zz; // aggiunta
    }
    // ******************************* FINE PARTE CORRETTA ***************************************
    else {
      // SSD: if iz==0 ---> N side; if iz==1 P side
      if (ssdone[j] == 0) {
        ssdone[j]=1;
        pair[j]=-1;
        Bool_t pside = (iz == 1);
        Bool_t impaired = kTRUE;
        Int_t pstrip = 0;
        Int_t nstrip = 0;
        if (pside) pstrip = ix; else nstrip = ix;
        for (Int_t k = 0; k < digits_num; k++) {
          if (ssdone[k] == 0 && impaired) {
            AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k);
            if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
              ssdone[k]=2;
              pair[j]=k;
              if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
              impaired=kFALSE;
            }
          }
        }
        if (!impaired) seg->GetPadCxz(pstrip, nstrip, X[j], Z[j]);
        X[j] /= 10000.0;  // convert microns to cm
        Z[j] /= 10000.0;  // convert microns to cm
      }
    }
  }
  return digits_num;
}