[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDanalyseRecPoints.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TRD/AliTRDcluster.h>
#include <AliRunLoader.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TLine.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TPad.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TSystem.h>
#include <TTree.h>
#endif

//#define EVENT_BY_EVENT      // Show amplitude event-by-event
//#define EVENT_BY_EVENT2     // Show tracks event-by-event
//#define PRINT_OTHER         // Display various additional histograms

// Number of noisy detectors
const Int_t numNoisyDets = 10;

// Holds the detector number of detectors that produce too much noise and
// are therefore excluded from the calculation
const Int_t exDets[numNoisyDets] = {25, 511, 515, 538, 523, 249, 294, 27, 24, 247};   

Bool_t isBadCluster(Int_t det, Int_t Q, Int_t PadCol, Int_t PadRow)
{  
  if (Q < 20) return 1;

  // Exclude bad cols
  if (PadCol == 142 || PadCol == 136 || PadCol == 66 || PadCol == 65 || PadCol <= 2) return 1;
  //if (PadCol == 142 || PadCol == 66 || PadCol == 65 || PadCol <= 2) return 1;
  //if (PadCol == 142 || PadCol <= 2) return 1;
  
  // Exclude noisy detectors
  for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
  {
    if (det == exDets[indDet])
    {
      // Exclude several regions of bad detectors
      if ((PadCol >= 134 && PadRow <= 5) || (PadCol == 141 && PadRow >= 10)) return 1;
      //if (PadRow == 15 || (PadCol >= 134 && PadRow <= 11) || (PadCol >= 119 && PadCol <= 123 && PadRow == 11)) return 1; 

      break;
    }
  }
  
  // Cluster seems to be clean
  return 0;
}


//=====================================

void AliTRDanalyseRecPoints (const char *filename="galice.root") 
{
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);

  // open run loader and load gAlice, kinematics and header
  AliRunLoader* runLoader = AliRunLoader::Open(filename);
  if (!runLoader) 
  {
    printf("Error: readKine:\nGetting run loader from file \"%s\" failed", filename);
    return;
  }

  runLoader->LoadHeader();
  runLoader->LoadRecPoints("TRD");
  TObjArray *module = new TObjArray();

  runLoader->CdGAFile();
  

  // Flag indicating whether the signal comes from a noisy detector
  Bool_t fromNoisyDet = 0;

  // Counts the number of clusters per event
  Int_t clusters = 0;

  // Threshold for the amplitude
  const Int_t kQThreshold = 30;

  // Selection criterions for clusters
  const Int_t kThresholdLower = 30;
  const Int_t kThresholdUpper = 500;

  // Number of events passing the criterion
  Int_t numPassed = 0;

  // hist
#ifdef PRINT_OTHER
  TH1D *mNCls = new TH1D("NCls", ";number of clusters", 400, 0, 400);
  TH1D *mResY = new TH1D("resY", ";sigma Y (cm)", 100, 0, 1);
  TH1D *mResZ = new TH1D("resZ", ";sigma Z (cm)", 100, 1, 3);

  TH1D *mNPads = new TH1D("NPads", ";number of pads", 10, -0.5, 9.5);
  TH1D *mCenter = new TH1D("center", ";center", 100, -1, 1);
  TH1D *mTime = new TH1D("timeBin", ";time bin", 25, -0.5, 24.5);
  TH1D *mPlane = new TH1D("plane", ";plane", 6, -0.5, 5.5);

  TH1D *mX = new TH1D("x", "; X (cm)", 100, 0, 3.5);
  TH1D *mY = new TH1D("y", "; y / chamber width (%)", 100, -1, 1);
  TH1D *mZ = new TH1D("z", "; Z (cm)", 200, -400, 400);

  TH2D *mPlaneY = new TH2D("planeY", ";plane;Y (cm);#Clusters", 6, -0.5, 5.5, 100, -60, 60);
  TH2D *mYX = new TH2D("yx", ";Y (cm);X (cm);#Clusters", 100, -60, 60, 120, 280, 400);
#endif

  TH1D *mChrg = new TH1D("chrg", ";Charge(no pad-filter)", 200, 0, 200);
  TH1D *mChrg2 = new TH1D("chrg2", ";Charge(pads<9)", 200, 0, 200);
  TH1D *mChrg3 = new TH1D("chrg3", ";Charge(2<pads<9)", 200, 0, 200);
  
  TH2D *mQvsPads = new TH2D("QvsPads", ";Number of pads;Q;#Clusters", 15, -0.5, 14.5, 100, 0, 1000);

  // There are 540 detectors
  TH2D *mQvsDetector = new TH2D("QvsDet", ";Number of detector;Q;#Clusters", 540, -0.5, 539.5, 100, 0, 1000);

  // Supermoduleid=int(detectorNum / 30)
  TH2D *mQvsSModId = new TH2D("QvsSModId", ";Supermodule id;Q;#Clusters", 18, -0.5, 17.5, 100, 0, 1000);
  // Histograms for the already installed supermodules
  TH1D *mSMod0 = new TH1D("SMod0", ";Charge", 200, 0, 200);
  TH1D *mSMod8 = new TH1D("SMod8", ";Charge", 200, 0, 200);
  TH1D *mSMod9 = new TH1D("SMod9", ";Charge", 200, 0, 200);
  TH1D *mSMod17 = new TH1D("SMod17", ";Charge", 200, 0, 200);

  // Index + 2 corresponds to the number of activated pads
  TH1D *mPads[14];
  char name[8];
  for (Int_t ind = 0; ind < 14; ind++)
  {
    sprintf(name,"%dPads",ind + 2);
    mPads[ind] = new TH1D(name, ";Charge", 200, 0, 200);
  }

  // Amplitude distribution will be displayed event by event for the first 100 events
  TH1D *mChrgEBE = new TH1D("chrgEBE", ";Charge", 200, 0, 200);

  // Keep track of the number of clusters per event
  TH1D *mClusters = new TH1D("Clusters", ";Number of clusters per event", 300, -0.5, 299.5);

  // Keep track of the number of clusters per chamber with Q < kQThreshold
  char clusTitle[55];
  sprintf(clusTitle,";Number of clusters with Q < %d vs detector", kQThreshold);
  TH1D *mClusCham = new TH1D("ClusCham", clusTitle, 540, -0.5, 539.5);
  // Same but this time noisy detectors excluded
  sprintf(clusTitle,";Number of clusters with Q < %d vs detector", kQThreshold);
  TH1D *mClusChamEx = new TH1D("ClusChamEx", clusTitle, 540, -0.5, 539.5);

  TH2D *mColRow = new TH2D("ColRow", ";Column;Row;#Clusters", 144, -0.5, 143.5, 18, -0.5, 17.5);
  TH2D *mColRowFiltered = new TH2D("ColRowFiltered", ";Column;Row;#Clusters", 144, -0.5, 143.5, 18, -0.5, 17.5);

  TH2D *mColTime = new TH2D("ColTime", ";TimeBin;Column", 18, -0.5, 17.5, 144, -0.5, 143.5);
  TH2D *mColTimeFiltered = new TH2D("ColTimeFiltered", ";TimeBin;Column;#Clusters", 18, -0.5, 17.5, 144, -0.5, 143.5);

  // YX histograms for the already installed supermodules
  TH2D *mTrSMod0 = new TH2D("TrSMod0", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
  TH2D *mTrSMod8 = new TH2D("TrSMod8", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
  TH2D *mTrSMod9 = new TH2D("TrSMod9", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
  TH2D *mTrSMod17 = new TH2D("TrSMod17", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);


#ifdef EVENT_BY_EVENT  
  TCanvas* eventc = new TCanvas("eventc","Amplitude distribution for event 0",700,500);
  char title[40];
#endif

#ifdef EVENT_BY_EVENT2
  TCanvas* sTrc = new TCanvas("sTrc","Track -  Supermodules",800,800);
  sTrc->SetLogz();
  sTrc->Divide(2,2);
  char title2[40];
#endif

 
  int nEvents = runLoader->GetNumberOfEvents();
  
  AliTRDcluster *cls = 0;

  for(Int_t ev = 0; ev < nEvents - 1; ev++) 
  {
    clusters = 0;    

    TTree *tree = runLoader->GetTreeR("TRD", 0);
    tree->SetBranchAddress("TRDcluster", &module);

    int N = tree->GetEntries();

    // Check number of clusters
    for(Int_t ind = 0; ind < N; ind++)
    {
      tree->GetEntry(ind);
      Int_t m = module->GetEntries();
  
      for (Int_t j = 0; j < m; j++)
      {
        if (cls != 0) delete cls;
        cls = (AliTRDcluster*)module->At(j);
        if (!isBadCluster(cls->GetDetector(), (Int_t)cls->GetQ(), cls->GetPadCol(), cls->GetPadRow())) clusters++;
        else
	      {
          // Uncomment the following lines AND comment the next if-clause to display the histograms for the 
          // EXCLUDED clusters

	        //mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
	        //mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
	        //mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());	
	        //mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
          //mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());    
	        //mClusCham->Fill(cls->GetDetector());
          //mYX->Fill(cls->GetY(), cls->GetX());
        }
      }  
    }

    // Process event, if number of clusters is in the correct range
    if (clusters >= kThresholdLower && clusters <= kThresholdUpper)
    {
      clusters = 0;
      numPassed++;

      for (Int_t i = 0; i < N; i++) 
      {
	      tree->GetEntry(i);
	      int m = module->GetEntries();

#ifdef PRINT_OTHER
	      mNCls->Fill(m);
#endif

	      for(Int_t j = 0; j < m; j++) 
	      {
          fromNoisyDet = 0;
	        clusters++;

          if (cls != 0) delete cls;
	        cls = (AliTRDcluster*)module->At(j);
	      
	        // Filter the clusters of the events, too
          if (cls->GetQ() < 20) continue;
           
	        mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());

          if (cls->GetQ() <= kQThreshold) 
          {
            mClusCham->Fill(cls->GetDetector());
	        }
          
          // Exclude noisy detectors
          for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
	        {
	          if (cls->GetDetector() == exDets[indDet])
	          {
          	  fromNoisyDet = 1;
	            break;
	          }
	        }

          // Exclude I
          if (cls->GetPadCol() == 142 || cls->GetPadCol() == 136 || cls->GetPadCol() == 66 || cls->GetPadCol() == 65 || cls->GetPadCol() <= 2) continue;
          //if (cls->GetPadCol() == 142 || cls->GetPadCol() == 66 || cls->GetPadCol() == 65 || cls->GetPadCol() <= 2) continue;
          //if (cls->GetPadCol() == 142 || cls->GetPadCol() <= 2) continue;

	        if (fromNoisyDet) 
	        {
	          mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());    
 
	          // Exclude II
	          //if (cls->GetPadRow() == 15 || (cls->GetPadCol() >= 134 && cls->GetPadRow() <= 11) || (cls->GetPadCol() >= 119 && cls->GetPadCol() <= 123 && cls->GetPadRow() == 11)) continue; 
            
            if ((cls->GetPadCol() >= 134 && cls->GetPadRow() <= 5) || (cls->GetPadCol() == 141 && cls->GetPadRow() >= 10)) continue; 
            
            mColRowFiltered->Fill(cls->GetPadCol(), cls->GetPadRow());
	        }

          mColTimeFiltered->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());

          if (cls->GetQ() <= kQThreshold) 
          {
            mClusChamEx->Fill(cls->GetDetector());
	        }
	        
#ifdef PRINT_OTHER
	        mX->Fill(cls->GetX());
	        mZ->Fill(cls->GetZ());
          mYX->Fill(cls->GetY(), cls->GetX());
	        mResY->Fill(TMath::Sqrt(cls->GetSigmaY2()));
	        mResZ->Fill(TMath::Sqrt(cls->GetSigmaZ2()));
          mNPads->Fill(cls->GetNPads());
          mCenter->Fill(cls->GetCenter());
	        mTime->Fill(cls->GetLocalTimeBin());
#endif

	        if (cls->GetNPads() > 2 && cls->GetNPads() < 9) mChrg3->Fill(cls->GetQ());
          if (cls->GetNPads() < 9) mChrg2->Fill(cls->GetQ());
          mChrg->Fill(cls->GetQ());	         
	        
	        mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
	        mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
	        mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());	

          Int_t superModule = (Int_t)(cls->GetDetector() / 30.);

          // Position in supermodule = cls->GetDetector() - superModule * 30
          // There are 6 layers per supermodule:
          Int_t layer = (cls->GetDetector() - superModule * 30) % 6;
	        switch (superModule)
	        {
	          case 0:
	            mSMod0->Fill(cls->GetQ());
                    mTrSMod0->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
	            break;
	          case 8:
	            mSMod8->Fill(cls->GetQ());
                    mTrSMod8->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
	            break;
	          case 9:
	            mSMod9->Fill(cls->GetQ());
                    mTrSMod9->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
	            break;
	          case 17:
	            mSMod17->Fill(cls->GetQ());
                    mTrSMod17->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
	            break;
	          default:
	            break;
	        }

	        // Only consider events with maximal 15 pads here
	        if (cls->GetNPads() - 2 < 14 && cls->GetNPads() - 2 >= 0)
	        {
	          mPads[cls->GetNPads() - 2]->Fill(cls->GetQ());
	        }

	        mChrgEBE->Fill(cls->GetQ());
	      }
      }
      

#ifdef EVENT_BY_EVENT  
      // Event by event display for the first 100 events
      if (numPassed < 100)
      {
        sprintf(title,"Charge distribution for event %d", ev);
        eventc->SetTitle(title);
        mChrgEBE->Draw();
        eventc->Update();
        gSystem->Sleep(400);
        mChrgEBE->Reset();
      }
#endif

#ifdef EVENT_BY_EVENT2
      // Event by event display 
      sprintf(title2,"Supermodule tracks for event %d", ev);
      sTrc->SetTitle(title2);
      sTrc->cd(1);
      mTrSMod0->Draw("colz");
      sTrc->cd(2);
      mTrSMod8->Draw("colz");
      sTrc->cd(3);
      mTrSMod9->Draw("colz");
      sTrc->cd(4);
      mTrSMod17->Draw("colz");
      sTrc->Update();
      gSystem->Sleep(400);
      mTrSMod0->Reset();
      mTrSMod8->Reset();
      mTrSMod9->Reset();
      mTrSMod17->Reset();
#endif


      mClusters->Fill(clusters);    
    }
     
    runLoader->GetNextEvent();
  }

#ifdef PRINT_OTHER
  new TCanvas();
  gPad->SetLogy();
  mNCls->Draw();

  new TCanvas();
  gPad->SetLogy();
  mResY->Draw();

  new TCanvas();
  gPad->SetLogy();
  mResZ->Draw();
#endif

  new TCanvas("chrgc", "Charge(no pad-filter)", 700, 500);
  mChrg->Draw();
  // Try landau-fit
  TF1* lFit = new TF1("lFit", "landau", 20, 200);
  cout << "\n\nFit (no pad-filter):\n";
  mChrg->Fit(lFit, "", "", 20, 200);

  new TCanvas("chrg2c", "Charge(pads<9)", 700, 500);
  mChrg2->Draw();
  // Try landau-fit
  TF1* lFit2 = new TF1("lFit2", "landau", 20, 200);
  cout << "\n\nFit (pads<9):\n";
  mChrg2->Fit(lFit2, "", "", 20, 200);

  new TCanvas("chrg3c", "Charge(2<pads<9)", 700, 500);
  mChrg3->Draw();
  // Try landau-fit
  TF1* lFit3 = new TF1("lFit3", "landau", 20, 200);
  cout << "\n\nFit (2<pads<9):\n";
  mChrg3->Fit(lFit3, "", "", 20, 200);

#ifdef PRINT_OTHER
  new TCanvas();
  mNPads->Draw();

  new TCanvas();
  mTime->Draw();

  new TCanvas();
  mX->Draw();

  new TCanvas();
  mY->Draw();

  new TCanvas();
  mZ->Draw();

  new TCanvas();
  mCenter->Draw();

  new TCanvas();
  mYX->Draw("colz");

  new TCanvas();
  mPlane->Draw();

  new TCanvas();
  mPlaneY->Draw("colz");
#endif

  new TCanvas();
  gPad->SetLogz();
  mQvsPads->Draw("colz");  

  new TCanvas();
  gPad->SetLogz();
  mQvsDetector->Draw("colz");
  TLine* tl = 0;
  for (Int_t num = 30; num < 540; num += 30)
  {
    // Draw vertical lines every 30 bins to sketch the layer modules
    tl = new TLine(num - 0.5, 0, num - 0.5, 1000);
    tl->Draw();
  }

  new TCanvas();
  gPad->SetLogz();
  mQvsSModId->Draw("colz");

  TCanvas* sc = new TCanvas("sc","Supermodules", 800, 800);
  sc->SetLogz();
  sc->Divide(2,2);
  sc->cd(1);
  mSMod0->Draw();
  sc->cd(2);
  mSMod8->Draw();
  sc->cd(3);
  mSMod9->Draw();
  sc->cd(4);
  mSMod17->Draw();
 

  TCanvas* padc = new TCanvas("padc","Pads", 1000, 1000);
  padc->Divide(4,4);
  for (Int_t ind = 0; ind < 14; ind++)
  {
    padc->cd(ind + 1);
    mPads[ind]->Draw();
  }

  new TCanvas();
  mClusters->Draw();

  TCanvas* clusChamc = new TCanvas("clusChamc", "Clusters against detector with threshold", 800, 800);
  clusChamc->Divide(1,2);
  clusChamc->cd(1);
  mClusCham->Draw();
  clusChamc->cd(2);
  mClusChamEx->Draw();

  TCanvas* noisyc = new TCanvas("noisyc", "Col/Row of noisy detectors");
  noisyc->Divide(1,2);
  noisyc->cd(1);
  gPad->SetLogz();
  mColRow->Draw("colz");
  noisyc->cd(2);
  gPad->SetLogz();
  mColRowFiltered->Draw("colz");

  TCanvas* ColTimec = new TCanvas("ColTimec", "Col/Time");
  ColTimec->Divide(1,2);
  ColTimec->cd(1);
  gPad->SetLogz();
  mColTime->Draw("colz");
  ColTimec->cd(2);
  gPad->SetLogz();
  mColTimeFiltered->Draw("colz");

  cout << numPassed << " out of " << nEvents << " events passed the criterion (number of 'good' clusters >= ";
  cout << kThresholdLower << " && clusters <= " << kThresholdUpper << ").\n";
}
Commit	Line	Data
68252c5c	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <TRD/AliTRDcluster.h>
	3	#include <AliRunLoader.h>
	4	#include <TCanvas.h>
	5	#include <TF1.h>
	6	#include <TH1.h>
	7	#include <TH2.h>
	8	#include <TLine.h>
	9	#include <TMath.h>
	10	#include <TObjArray.h>
	11	#include <TPad.h>
	12	#include <TROOT.h>
	13	#include <TStyle.h>
	14	#include <TSystem.h>
	15	#include <TTree.h>
	16	#endif
	17
	18	//#define EVENT_BY_EVENT // Show amplitude event-by-event
	19	//#define EVENT_BY_EVENT2 // Show tracks event-by-event
	20	//#define PRINT_OTHER // Display various additional histograms
	21
	22	// Number of noisy detectors
	23	const Int_t numNoisyDets = 10;
	24
	25	// Holds the detector number of detectors that produce too much noise and
	26	// are therefore excluded from the calculation
	27	const Int_t exDets[numNoisyDets] = {25, 511, 515, 538, 523, 249, 294, 27, 24, 247};
	28
	29	Bool_t isBadCluster(Int_t det, Int_t Q, Int_t PadCol, Int_t PadRow)
	30	{
	31	if (Q < 20) return 1;
	32
	33	// Exclude bad cols
	34	if (PadCol == 142 \|\| PadCol == 136 \|\| PadCol == 66 \|\| PadCol == 65 \|\| PadCol <= 2) return 1;
	35	//if (PadCol == 142 \|\| PadCol == 66 \|\| PadCol == 65 \|\| PadCol <= 2) return 1;
	36	//if (PadCol == 142 \|\| PadCol <= 2) return 1;
	37
	38	// Exclude noisy detectors
	39	for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
	40	{
	41	if (det == exDets[indDet])
	42	{
	43	// Exclude several regions of bad detectors
	44	if ((PadCol >= 134 && PadRow <= 5) \|\| (PadCol == 141 && PadRow >= 10)) return 1;
	45	//if (PadRow == 15 \|\| (PadCol >= 134 && PadRow <= 11) \|\| (PadCol >= 119 && PadCol <= 123 && PadRow == 11)) return 1;
	46
	47	break;
	48	}
	49	}
	50
	51	// Cluster seems to be clean
	52	return 0;
	53	}
	54
	55
	56	//=====================================
	57
	58	void AliTRDanalyseRecPoints (const char *filename="galice.root")
	59	{
	60	gROOT->SetStyle("Plain");
	61	gStyle->SetPalette(1);
	62
	63	// open run loader and load gAlice, kinematics and header
	64	AliRunLoader* runLoader = AliRunLoader::Open(filename);
65	if (!runLoader)
66	{
67	printf("Error: readKine:\nGetting run loader from file \"%s\" failed", filename);
68	return;
69	}
70
71	runLoader->LoadHeader();
72	runLoader->LoadRecPoints("TRD");
73	TObjArray *module = new TObjArray();
74
75	runLoader->CdGAFile();
76
77
78	// Flag indicating whether the signal comes from a noisy detector
79	Bool_t fromNoisyDet = 0;
80
81	// Counts the number of clusters per event
82	Int_t clusters = 0;
83
84	// Threshold for the amplitude
85	const Int_t kQThreshold = 30;
86
87	// Selection criterions for clusters
88	const Int_t kThresholdLower = 30;
89	const Int_t kThresholdUpper = 500;
90
91	// Number of events passing the criterion
92	Int_t numPassed = 0;
93
94	// hist
95	#ifdef PRINT_OTHER
96	TH1D *mNCls = new TH1D("NCls", ";number of clusters", 400, 0, 400);
97	TH1D *mResY = new TH1D("resY", ";sigma Y (cm)", 100, 0, 1);
98	TH1D *mResZ = new TH1D("resZ", ";sigma Z (cm)", 100, 1, 3);
99
100	TH1D *mNPads = new TH1D("NPads", ";number of pads", 10, -0.5, 9.5);
101	TH1D *mCenter = new TH1D("center", ";center", 100, -1, 1);
102	TH1D *mTime = new TH1D("timeBin", ";time bin", 25, -0.5, 24.5);
103	TH1D *mPlane = new TH1D("plane", ";plane", 6, -0.5, 5.5);
104
105	TH1D *mX = new TH1D("x", "; X (cm)", 100, 0, 3.5);
106	TH1D *mY = new TH1D("y", "; y / chamber width (%)", 100, -1, 1);
107	TH1D *mZ = new TH1D("z", "; Z (cm)", 200, -400, 400);
108
109	TH2D *mPlaneY = new TH2D("planeY", ";plane;Y (cm);#Clusters", 6, -0.5, 5.5, 100, -60, 60);
110	TH2D *mYX = new TH2D("yx", ";Y (cm);X (cm);#Clusters", 100, -60, 60, 120, 280, 400);
111	#endif
112
113	TH1D *mChrg = new TH1D("chrg", ";Charge(no pad-filter)", 200, 0, 200);
114	TH1D *mChrg2 = new TH1D("chrg2", ";Charge(pads<9)", 200, 0, 200);
115	TH1D *mChrg3 = new TH1D("chrg3", ";Charge(2<pads<9)", 200, 0, 200);
116
117	TH2D *mQvsPads = new TH2D("QvsPads", ";Number of pads;Q;#Clusters", 15, -0.5, 14.5, 100, 0, 1000);
118
119	// There are 540 detectors
120	TH2D *mQvsDetector = new TH2D("QvsDet", ";Number of detector;Q;#Clusters", 540, -0.5, 539.5, 100, 0, 1000);
121
122	// Supermoduleid=int(detectorNum / 30)
123	TH2D *mQvsSModId = new TH2D("QvsSModId", ";Supermodule id;Q;#Clusters", 18, -0.5, 17.5, 100, 0, 1000);
124	// Histograms for the already installed supermodules
125	TH1D *mSMod0 = new TH1D("SMod0", ";Charge", 200, 0, 200);
126	TH1D *mSMod8 = new TH1D("SMod8", ";Charge", 200, 0, 200);
127	TH1D *mSMod9 = new TH1D("SMod9", ";Charge", 200, 0, 200);
128	TH1D *mSMod17 = new TH1D("SMod17", ";Charge", 200, 0, 200);
129
130	// Index + 2 corresponds to the number of activated pads
131	TH1D *mPads[14];
132	char name[8];
133	for (Int_t ind = 0; ind < 14; ind++)
134	{
135	sprintf(name,"%dPads",ind + 2);
136	mPads[ind] = new TH1D(name, ";Charge", 200, 0, 200);
137	}
138
139	// Amplitude distribution will be displayed event by event for the first 100 events
140	TH1D *mChrgEBE = new TH1D("chrgEBE", ";Charge", 200, 0, 200);
141
142	// Keep track of the number of clusters per event
143	TH1D *mClusters = new TH1D("Clusters", ";Number of clusters per event", 300, -0.5, 299.5);
144
145	// Keep track of the number of clusters per chamber with Q < kQThreshold
146	char clusTitle[55];
147	sprintf(clusTitle,";Number of clusters with Q < %d vs detector", kQThreshold);
148	TH1D *mClusCham = new TH1D("ClusCham", clusTitle, 540, -0.5, 539.5);
149	// Same but this time noisy detectors excluded
150	sprintf(clusTitle,";Number of clusters with Q < %d vs detector", kQThreshold);
151	TH1D *mClusChamEx = new TH1D("ClusChamEx", clusTitle, 540, -0.5, 539.5);
152
153	TH2D *mColRow = new TH2D("ColRow", ";Column;Row;#Clusters", 144, -0.5, 143.5, 18, -0.5, 17.5);
154	TH2D *mColRowFiltered = new TH2D("ColRowFiltered", ";Column;Row;#Clusters", 144, -0.5, 143.5, 18, -0.5, 17.5);
155
156	TH2D *mColTime = new TH2D("ColTime", ";TimeBin;Column", 18, -0.5, 17.5, 144, -0.5, 143.5);
157	TH2D *mColTimeFiltered = new TH2D("ColTimeFiltered", ";TimeBin;Column;#Clusters", 18, -0.5, 17.5, 144, -0.5, 143.5);
158
159	// YX histograms for the already installed supermodules
160	TH2D mTrSMod0 = new TH2D("TrSMod0", ";Col;time+layer20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
161	TH2D mTrSMod8 = new TH2D("TrSMod8", ";Col;time+layer20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
162	TH2D mTrSMod9 = new TH2D("TrSMod9", ";Col;time+layer20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
163	TH2D mTrSMod17 = new TH2D("TrSMod17", ";Col;time+layer20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
164
165
166	#ifdef EVENT_BY_EVENT
167	TCanvas* eventc = new TCanvas("eventc","Amplitude distribution for event 0",700,500);
168	char title[40];
169	#endif
170
171	#ifdef EVENT_BY_EVENT2
172	TCanvas* sTrc = new TCanvas("sTrc","Track - Supermodules",800,800);
173	sTrc->SetLogz();
174	sTrc->Divide(2,2);
175	char title2[40];
176	#endif
177
178
179	int nEvents = runLoader->GetNumberOfEvents();
180
181	AliTRDcluster *cls = 0;
182
183	for(Int_t ev = 0; ev < nEvents - 1; ev++)
184	{
185	clusters = 0;
186
187	TTree *tree = runLoader->GetTreeR("TRD", 0);
188	tree->SetBranchAddress("TRDcluster", &module);
189
190	int N = tree->GetEntries();
191
192	// Check number of clusters
193	for(Int_t ind = 0; ind < N; ind++)
194	{
195	tree->GetEntry(ind);
196	Int_t m = module->GetEntries();
197
198	for (Int_t j = 0; j < m; j++)
199	{
200	if (cls != 0) delete cls;
201	cls = (AliTRDcluster*)module->At(j);
202	if (!isBadCluster(cls->GetDetector(), (Int_t)cls->GetQ(), cls->GetPadCol(), cls->GetPadRow())) clusters++;
203	else
204	{
205	// Uncomment the following lines AND comment the next if-clause to display the histograms for the
206	// EXCLUDED clusters
207
208	//mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
209	//mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
210	//mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());
211	//mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
212	//mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());
213	//mClusCham->Fill(cls->GetDetector());
214	//mYX->Fill(cls->GetY(), cls->GetX());
215	}
216	}
217	}
218
219	// Process event, if number of clusters is in the correct range
220	if (clusters >= kThresholdLower && clusters <= kThresholdUpper)
221	{
222	clusters = 0;
223	numPassed++;
224
225	for (Int_t i = 0; i < N; i++)
226	{
227	tree->GetEntry(i);
228	int m = module->GetEntries();
229
230	#ifdef PRINT_OTHER
231	mNCls->Fill(m);
232	#endif
233
234	for(Int_t j = 0; j < m; j++)
235	{
236	fromNoisyDet = 0;
237	clusters++;
238
239	if (cls != 0) delete cls;
240	cls = (AliTRDcluster*)module->At(j);
241
242	// Filter the clusters of the events, too
243	if (cls->GetQ() < 20) continue;
244
245	mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
246
247	if (cls->GetQ() <= kQThreshold)
248	{
249	mClusCham->Fill(cls->GetDetector());
250	}
251
252	// Exclude noisy detectors
253	for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
254	{
255	if (cls->GetDetector() == exDets[indDet])
256	{
257	fromNoisyDet = 1;
258	break;
259	}
260	}
261
262	// Exclude I
263	if (cls->GetPadCol() == 142 \|\| cls->GetPadCol() == 136 \|\| cls->GetPadCol() == 66 \|\| cls->GetPadCol() == 65 \|\| cls->GetPadCol() <= 2) continue;
264	//if (cls->GetPadCol() == 142 \|\| cls->GetPadCol() == 66 \|\| cls->GetPadCol() == 65 \|\| cls->GetPadCol() <= 2) continue;
265	//if (cls->GetPadCol() == 142 \|\| cls->GetPadCol() <= 2) continue;
266
267	if (fromNoisyDet)
268	{
269	mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());
270
271	// Exclude II
272	//if (cls->GetPadRow() == 15 \|\| (cls->GetPadCol() >= 134 && cls->GetPadRow() <= 11) \|\| (cls->GetPadCol() >= 119 && cls->GetPadCol() <= 123 && cls->GetPadRow() == 11)) continue;
273
274	if ((cls->GetPadCol() >= 134 && cls->GetPadRow() <= 5) \|\| (cls->GetPadCol() == 141 && cls->GetPadRow() >= 10)) continue;
275
276	mColRowFiltered->Fill(cls->GetPadCol(), cls->GetPadRow());
277	}
278
279	mColTimeFiltered->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
280
281	if (cls->GetQ() <= kQThreshold)
282	{
283	mClusChamEx->Fill(cls->GetDetector());
284	}
285
286	#ifdef PRINT_OTHER
287	mX->Fill(cls->GetX());
288	mZ->Fill(cls->GetZ());
289	mYX->Fill(cls->GetY(), cls->GetX());
290	mResY->Fill(TMath::Sqrt(cls->GetSigmaY2()));
291	mResZ->Fill(TMath::Sqrt(cls->GetSigmaZ2()));
292	mNPads->Fill(cls->GetNPads());
293	mCenter->Fill(cls->GetCenter());
294	mTime->Fill(cls->GetLocalTimeBin());
295	#endif
296
297	if (cls->GetNPads() > 2 && cls->GetNPads() < 9) mChrg3->Fill(cls->GetQ());
298	if (cls->GetNPads() < 9) mChrg2->Fill(cls->GetQ());
299	mChrg->Fill(cls->GetQ());
300
301	mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
302	mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
303	mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());
304
305	Int_t superModule = (Int_t)(cls->GetDetector() / 30.);
306
307	// Position in supermodule = cls->GetDetector() - superModule * 30
308	// There are 6 layers per supermodule:
309	Int_t layer = (cls->GetDetector() - superModule * 30) % 6;
310	switch (superModule)
311	{
312	case 0:
313	mSMod0->Fill(cls->GetQ());
314	mTrSMod0->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
315	break;
316	case 8:
317	mSMod8->Fill(cls->GetQ());
318	mTrSMod8->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
319	break;
320	case 9:
321	mSMod9->Fill(cls->GetQ());
322	mTrSMod9->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
323	break;
324	case 17:
325	mSMod17->Fill(cls->GetQ());
326	mTrSMod17->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
327	break;
328	default:
329	break;
330	}
331
332	// Only consider events with maximal 15 pads here
333	if (cls->GetNPads() - 2 < 14 && cls->GetNPads() - 2 >= 0)
334	{
335	mPads[cls->GetNPads() - 2]->Fill(cls->GetQ());
336	}
337
338	mChrgEBE->Fill(cls->GetQ());
339	}
340	}
341
342
343	#ifdef EVENT_BY_EVENT
344	// Event by event display for the first 100 events
345	if (numPassed < 100)
346	{
347	sprintf(title,"Charge distribution for event %d", ev);
348	eventc->SetTitle(title);
349	mChrgEBE->Draw();
350	eventc->Update();
351	gSystem->Sleep(400);
352	mChrgEBE->Reset();
353	}
354	#endif
355
356	#ifdef EVENT_BY_EVENT2
357	// Event by event display
358	sprintf(title2,"Supermodule tracks for event %d", ev);
359	sTrc->SetTitle(title2);
360	sTrc->cd(1);
361	mTrSMod0->Draw("colz");
362	sTrc->cd(2);
363	mTrSMod8->Draw("colz");
364	sTrc->cd(3);
365	mTrSMod9->Draw("colz");
366	sTrc->cd(4);
367	mTrSMod17->Draw("colz");
368	sTrc->Update();
369	gSystem->Sleep(400);
370	mTrSMod0->Reset();
371	mTrSMod8->Reset();
372	mTrSMod9->Reset();
373	mTrSMod17->Reset();
374	#endif
375
376
377	mClusters->Fill(clusters);
378	}
379
380	runLoader->GetNextEvent();
381	}
382
383	#ifdef PRINT_OTHER
384	new TCanvas();
385	gPad->SetLogy();
386	mNCls->Draw();
387
388	new TCanvas();
389	gPad->SetLogy();
390	mResY->Draw();
391
392	new TCanvas();
393	gPad->SetLogy();
394	mResZ->Draw();
395	#endif
396
397	new TCanvas("chrgc", "Charge(no pad-filter)", 700, 500);
398	mChrg->Draw();
399	// Try landau-fit
400	TF1* lFit = new TF1("lFit", "landau", 20, 200);
401	cout << "\n\nFit (no pad-filter):\n";
402	mChrg->Fit(lFit, "", "", 20, 200);
403
404	new TCanvas("chrg2c", "Charge(pads<9)", 700, 500);
405	mChrg2->Draw();
406	// Try landau-fit
407	TF1* lFit2 = new TF1("lFit2", "landau", 20, 200);
408	cout << "\n\nFit (pads<9):\n";
409	mChrg2->Fit(lFit2, "", "", 20, 200);
410
411	new TCanvas("chrg3c", "Charge(2<pads<9)", 700, 500);
412	mChrg3->Draw();
413	// Try landau-fit
414	TF1* lFit3 = new TF1("lFit3", "landau", 20, 200);
415	cout << "\n\nFit (2<pads<9):\n";
416	mChrg3->Fit(lFit3, "", "", 20, 200);
417
418	#ifdef PRINT_OTHER
419	new TCanvas();
420	mNPads->Draw();
421
422	new TCanvas();
423	mTime->Draw();
424
425	new TCanvas();
426	mX->Draw();
427
428	new TCanvas();
429	mY->Draw();
430
431	new TCanvas();
432	mZ->Draw();
433
434	new TCanvas();
435	mCenter->Draw();
436
437	new TCanvas();
438	mYX->Draw("colz");
439
440	new TCanvas();
441	mPlane->Draw();
442
443	new TCanvas();
444	mPlaneY->Draw("colz");
445	#endif
446
447	new TCanvas();
448	gPad->SetLogz();
449	mQvsPads->Draw("colz");
450
451	new TCanvas();
452	gPad->SetLogz();
453	mQvsDetector->Draw("colz");
454	TLine* tl = 0;
455	for (Int_t num = 30; num < 540; num += 30)
456	{
457	// Draw vertical lines every 30 bins to sketch the layer modules
458	tl = new TLine(num - 0.5, 0, num - 0.5, 1000);
459	tl->Draw();
460	}
461
462	new TCanvas();
463	gPad->SetLogz();
464	mQvsSModId->Draw("colz");
465
466	TCanvas* sc = new TCanvas("sc","Supermodules", 800, 800);
467	sc->SetLogz();
468	sc->Divide(2,2);
469	sc->cd(1);
470	mSMod0->Draw();
471	sc->cd(2);
472	mSMod8->Draw();
473	sc->cd(3);
474	mSMod9->Draw();
475	sc->cd(4);
476	mSMod17->Draw();
477
478
479	TCanvas* padc = new TCanvas("padc","Pads", 1000, 1000);
480	padc->Divide(4,4);
481	for (Int_t ind = 0; ind < 14; ind++)
482	{
483	padc->cd(ind + 1);
484	mPads[ind]->Draw();
485	}
486
487	new TCanvas();
488	mClusters->Draw();
489
490	TCanvas* clusChamc = new TCanvas("clusChamc", "Clusters against detector with threshold", 800, 800);
491	clusChamc->Divide(1,2);
492	clusChamc->cd(1);
493	mClusCham->Draw();
494	clusChamc->cd(2);
495	mClusChamEx->Draw();
496
497	TCanvas* noisyc = new TCanvas("noisyc", "Col/Row of noisy detectors");
498	noisyc->Divide(1,2);
499	noisyc->cd(1);
500	gPad->SetLogz();
501	mColRow->Draw("colz");
502	noisyc->cd(2);
503	gPad->SetLogz();
504	mColRowFiltered->Draw("colz");
505
506	TCanvas* ColTimec = new TCanvas("ColTimec", "Col/Time");
507	ColTimec->Divide(1,2);
508	ColTimec->cd(1);
509	gPad->SetLogz();
510	mColTime->Draw("colz");
511	ColTimec->cd(2);
512	gPad->SetLogz();
513	mColTimeFiltered->Draw("colz");
514
515	cout << numPassed << " out of " << nEvents << " events passed the criterion (number of 'good' clusters >= ";
516	cout << kThresholdLower << " && clusters <= " << kThresholdUpper << ").\n";
517	}