TRD/Macros/AliTRDanalyzeTestBeam.C

   1
   2 void AliTRDanalyzeTestBeam(Int_t run, Int_t begin, Int_t end) {
   3
   4   gROOT->SetStyle("Plain");
   5   //gStyle->SetPadTopMargin(0.02);
   6   //gStyle->SetPadRightMargin(0.02);
   7   //gStyle->SetPadLeftMargin(0.1);
   8
   9   gStyle->SetPalette(1);
  10   gStyle->SetOptStat(0);
  11   gStyle->SetOptDate();
  12
  13   TGaxis::SetMaxDigits(3);
  14
  15   TStopwatch st;
  16   st.Start();
  17
  18   const Int_t N = 640;
  19
  20   // declare histograms
  21   const int nbins = 100;
  22   const int nstart = 0;
  23
  24   TH1D *mSi1L = new TH1D("mSi1L", ";number of Si1 fired pads", 50, nstart-0.5, nstart+nbins-0.5);
  25   TH1D *mSi2L = new TH1D("mS21L", ";number of Si2 fired pads", 50, nstart-0.5, nstart+nbins-0.5);
  26
  27   TProfile *mSi1ChP = new TProfile("mSi1ChP",";Si1 pad number;signal", 1280, -0.5, 1279.5, 0, 200, "s");
  28   TProfile *mSi2ChP = new TProfile("mS21ChP",";Si2 pad number;signal", 1280, -0.5, 1279.5, 0, 200, "s");
  29
  30   TH1D *mSi1N = new TH1D("mSi1N", "Noise Dist Si1;ADC", 100, 0, 50);
  31   TH1D *mSi2N = new TH1D("mSi2N", "Noise Dist Si2;ADC", 100, 0, 50);
  32
  33   TH1D *mSiCh[4];
  34   mSiCh[0] = new TH1D("mSi1ChX", ";Si1X pad amplitude (ADC)", 250, -0.5, 249.5);
  35   mSiCh[1] = new TH1D("mSi1ChY", ";Si1Y pad amplitude (ADC)", 250, -0.5, 249.5);
  36   mSiCh[2] = new TH1D("mSi2ChX", ";Si2X pad amplitude (ADC)", 250, -0.5, 249.5);
  37   mSiCh[3] = new TH1D("mSi2ChY", ";Si2Y pad amplitude (ADC)", 250, -0.5, 249.5);
  38
  39   TH1D *mSiFullCh[4];
  40   mSiFullCh[0] = new TH1D("mSi1fChX", "Si1X;max amplitude (ADC)", 300, -0.5, 299.5);
  41   mSiFullCh[1] = new TH1D("mSi1fChY", "Si1Y;max amplitude (ADC)", 300, -0.5, 299.5);
  42   mSiFullCh[2] = new TH1D("mSi2fChX", "Si2X;max amplitude (ADC)", 300, -0.5, 299.5);
  43   mSiFullCh[3] = new TH1D("mSi2fChY", "Si2Y;max amplitude (ADC)", 300, -0.5, 299.5);
  44
  45   TH2D *mPos[2];
  46   mPos[0] = new TH2D("posSi1", ";x Si1 (mm);y Si1 (mm)", 128, 0, 32, 128, 0, 32);
  47   mPos[1] = new TH2D("posSi2", ";x Si2 (mm);y Si2 (mm)", 128, 0, 32, 128, 0, 32);
  48
  49   TH2D *mCor[2];
  50   mCor[0] = new TH2D("corX", ";Si1 X (mm);Si2 X (mm)", 128, 0, 32, 128, 0, 32);
  51   mCor[1] = new TH2D("corY", ";Si1 Y (mm);Si2 Y (mm)", 128, 0, 32, 128, 0, 32);
  52
  53   TH2D *mChCor[2];
  54   mChCor[0] = new TH2D("ChCorSi1", ";Si1 amp X;Si1 amp Y", 100, 0, 200, 100, 0, 200);
  55   mChCor[1] = new TH2D("ChCorSi2", ";Si2 amp X;Si2 amp Y", 100, 0, 200, 100, 0, 200);
  56
  57   gStyle->SetOptStat(11);
  58   TH1D *mPb   = new TH1D("mPb", ";Amp Pb (ADC)", 150, -0.5, 4499.5);
  59   TH1D *mCher = new TH1D("mCher", ";Amp Cherenkov (ADC)", 150, -0.5, 4499.5);
  60   TH2D *mPbCher = new TH2D("mPbCher", ";amp Cherenkov;amp Pb", 150, -0.5, 4499.5, 150, -0.5, 4599.5);
  61   // gStyle->SetOptStat(0);
  62
  63   // needed by the AliTRDRawStreamTB
  64   AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  65   AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
  66   calib->SetRun(0);
  67
  68   // TRD data monitoring
  69   TH1D *mDet = new TH1D("mDet", ";chamber", 20, -0.5, 19.5);
  70   TH1D *mROB = new TH1D("mROB", ";rob", 20, -0.5, 19.5);
  71   TH1D *mTRDsig = new TH1D("mTRDsig", ";trdSignal", 100, -0.5, 299.5);
  72
  73   //AliLog::SetClassDebugLevel("AliTRDRawStreamTB", 10);
  74
  75   int counter = 0;
  76   // for(Int_t run = 365; run < 386; run++) {
  77   //for(Int_t run = 369; run < 382; run++) {
  78   //for(int run = 387; run < 389; run++) {
  79
  80   // if (run == 389) continue;
  81     cout << run << endl;
  82
  83     for(Int_t fn=begin; fn<end; fn++) {
  84
  85       // connect to data
  86       //const char *base="/Users/radomski/data/1GeV/";
  87       const char *base = "./";
  88       const char *filename = Form("%s/run%d_gdc_daq09.%03d.raw",base,run,fn);
  89       if (gSystem->AccessPathName(filename)) continue;
  90       cout << filename << endl;
  91
  92       AliTRDtestBeam *data = new AliTRDtestBeam(filename);
  93
  94       // process data
  95       while (data->NextEvent()) {
  96
  97         if (!(counter%1000)) cout << "Event = " << counter << endl;
  98         counter++;
  99
 100         /*
 101         AliTRDRawStreamTB *tb = data->GetTRDrawStream();
 102         while(tb->Next()) {
 103           mROB->Fill(tb->GetROB());
 104           mDet->Fill(tb->GetDet());
 105           int *sig = tb->GetSignals();
 106           mTRDsig->Fill(sig[0]);
 107           mTRDsig->Fill(sig[1]);
 108           mTRDsig->Fill(sig[2]);
 109         }
 110         delete tb;
 111         */
 112
 113         mCher->Fill(data->GetCher());
 114         mPb->Fill(data->GetPb());
 115         mPbCher->Fill(data->GetCher(), data->GetPb());
 116
 117         mSi1L->Fill(data->GetNSi1());
 118         mSi2L->Fill(data->GetNSi2());
 119
 120         for(int i=0; i<data->GetNSi1(); i++) {
 121           Int_t q = data->GetSi1Charge(i);
 122           Int_t a = data->GetSi1Address(i);
 123           if (a == 0) continue; // noisy channels
 124           mSi1ChP->Fill(a, q);
 125           if (a < N) mSiCh[0]->Fill(q);
 126           else mSiCh[1]->Fill(q);
 127         }
 128
 129         for(int i=0; i<data->GetNSi2(); i++) {
 130           Int_t q = data->GetSi2Charge(i);
 131           Int_t a = data->GetSi2Address(i);
 132           if (a == 0 || a == 1279) continue; // noisy channels
 133           mSi2ChP->Fill(a, q);
 134           if (a < N) mSiCh[2]->Fill(q);
 135           else mSiCh[3]->Fill(q);
 136         }
 137
 138         mSiFullCh[0]->Fill(data->GetQx(0));
 139         mSiFullCh[1]->Fill(data->GetQy(0));
 140         mSiFullCh[2]->Fill(data->GetQx(1));
 141         mSiFullCh[3]->Fill(data->GetQy(1));
 142
 143         for(int k=0; k<2; k++)
 144           mChCor[k]->Fill(data->GetQx(k), data->GetQy(k));
 145
 146         /*
 147           if (data->GetQx(0) < 20) continue;
 148           if (data->GetQx(1) < 20) continue;
 149           if (data->GetQy(0) < 20) continue;
 150           if (data->GetQy(1) < 20) continue;
 151         */
 152
 153         for(int k=0; k<2; k++)
 154           mPos[k]->Fill(data->GetX(k), data->GetY(k));
 155
 156         mCor[0]->Fill(data->GetX(0), data->GetX(1));
 157         mCor[1]->Fill(data->GetY(0), data->GetY(1));
 158       }
 159
 160       delete data;
 161     }
 162     //}
 163
 164   // process histograms
 165   for(int i=1; i<1281; i++) mSi1N->Fill(mSi1ChP->GetBinError(i));
 166   for(int i=1; i<1281; i++) mSi2N->Fill(mSi2ChP->GetBinError(i));
 167
 168   // display
 169   cout << "Number of Events = " << counter << endl;
 170
 171   /**/
 172   TCanvas *c = new TCanvas("siliconSignal", "silicon signal");
 173   c->Divide(2,2, 0.01, 0.01);
 174   c->cd(1);
 175   mSi1L->Draw();
 176   c->cd(2);
 177   mSi2L->Draw();
 178   c->cd(3);
 179   mSi1ChP->Draw();
 180   c->cd(4);
 181   mSi2ChP->Draw();
 182   /* */
 183
 184   /**/
 185   c = new TCanvas();
 186   c->Divide(2,2,0.01,0.01);
 187   c->cd(1);
 188   mSi1N->Draw();
 189   c->cd(2);
 190   mSi2N->Draw();
 191   /**/
 192
 193   /**/
 194   // pads
 195   c = new TCanvas("siPads", "Silicon Pads");
 196   c->Divide(2,2, 0.01, 0.01);
 197   for(int i=0; i<4; i++) {
 198     c->cd(i+1);
 199     gPad->SetLogy();
 200     mSiCh[i]->Draw();
 201   }
 202
 203   // clusters
 204   c = new TCanvas("siCluster", "silicon clusters");
 205   c->Divide(2,2, 0.01, 0.01);
 206   for(int i=0; i<4; i++) {
 207     c->cd(i+1);
 208     gPad->SetLogy();
 209     mSiFullCh[i]->Draw();
 210   }
 211
 212   // position and correlation
 213   c = new TCanvas("siPosition", "reconstructed position");
 214   c->Divide(2,2, 0.01, 0.01);
 215   for(int i=0; i<2; i++) {
 216     c->cd(1+i);
 217     mPos[i]->Draw("col");
 218     c->cd(3+i);
 219     mCor[i]->Draw("col");
 220   }
 221
 222   c = new TCanvas("siCharge", "si charge correlation");
 223   c->Divide(2,2, 0.01, 0.01);
 224   c->cd(1);
 225   gPad->SetLogz();
 226   mChCor[0]->Draw("col");
 227   c->cd(2);
 228   gPad->SetLogz();
 229   mChCor[1]->Draw("col");
 230   /**/
 231
 232   new TCanvas();
 233   gPad->SetLogy();
 234   mCher->Draw();
 235
 236   // electron sample
 237   int bin = mCher->FindBin(500.);
 238   double ef = (mCher->Integral(bin, 151)/ mCher->GetSum());
 239   TLatex *l = new TLatex(2e3, 0.02*mCher->GetSum(), Form("Electron fraction = %.2f ", ef));
 240   l->Draw();
 241
 242   new TCanvas();
 243   gPad->SetLogy();
 244   mPb->Draw();
 245
 246   new TCanvas();
 247   gPad->SetLogz();
 248   mPbCher->Draw("colz");
 249
 250   /*
 251   c = new TCanvas();
 252   c->Divide(2,2,0.01, 0.01);
 253   c->cd(1);
 254   gPad->SetLogy();
 255   mTRDsig->Draw();
 256
 257   c->cd(2);
 258   mROB->Draw();
 259
 260   c->cd(3);
 261   mDet->Draw();
 262   */
 263
 264   st.Stop();
 265   st.Print();
 266 }
 267
 268 //Int_t SelectEvent() {
 269 //}