ITS/AliITSplotSPDClusters.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2
   3 // Standard includes
   4 #include <iostream.h>
   5 // Root includes
   6 #include <TFile.h>
   7 #include <TObjArray.h>
   8 #include <TClonesArray.h>
   9 #include <TTree.h>
  10 #include <TBranch.h>
  11 #include <TH1D.h>
  12 // AliRoot includes
  13 #include "STEER/AliRun.h"
  14 #include "ITS/AliITS.h"
  15 #include "ITS/AliITSRawCluster.h"
  16
  17 #endif
  18
  19 void AliITSplotSPDClusters(const char* filename="galice_80.root"){
  20     //------------------------------------------------------------------
  21     // This macro will read the TreeC produced during reconstruction and
  22     // plot the number and type of clusters found in the SPD.
  23     // root [0] .L AliITSplotSPDClusters.C  //this loads the macro in memory
  24     // root [1] AliITSplotSPDClusters();    //by default process first event
  25     // or
  26     // root [1] AliITSplotSPDClusters("galice2.root"); // process all events
  27     //                                                    from the file
  28     //                                                    galice2.root.
  29     // or
  30     // root [0] .x AliITSplotSPDClusters.C("galice2.root") // will all of the
  31     //                                                        events from the
  32     //                                                        galice2.root
  33     //------------------------------------------------------------------
  34
  35     // delete the current gAlice object, the one from input file will be used
  36     if(gAlice){
  37         delete gAlice;
  38         gAlice = 0;
  39     }else{ // Dynamically link some shared libs
  40         if(gClassTable->GetID("AliRun") < 0) {
  41             gROOT->LoadMacro("loadlibs.C");
  42             loadlibs();
  43         } // end if
  44     } // end if gAlice
  45
  46     // Connect the Root Galice file containing Geometry, Kine and Hits
  47     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
  48     if(!file) file = new TFile(filename);
  49     // Get AliRun object from file or create it if not on file
  50     if(!gAlice) {
  51         gAlice = (AliRun*)file->Get("gAlice");
  52         if(gAlice) cout << "AliRun object found on file" << endl;
  53         if(!gAlice) gAlice = new AliRun("gAlice","SPD Clulster analysis");
  54     } // end if !gAlice
  55     Int_t nevents = gAlice->GetEventsPerRun();
  56
  57     // Get pointers to the ITS Alice detectors.
  58     AliITS *ITS = (AliITS*) gAlice->GetDetector("ITS");
  59     if(ITS==0){
  60         cout << "ITS not found. Exiting." << endl;
  61         return;
  62     } // end if ITS==0
  63     AliITSresponseSPDdubna *resp0 = ITS->DetType(0)->GetResponseModel();
  64     Float_t diffCoeff= resp0->DistanceOverVoltage();//Get Value of Diffusion Coefficient parameter d/v.
  65     diffCoeff *= resp0->Temperature();
  66     diffCoeff = TMath::Sqrt(diffCoeff*8.6173376e-5);// K_Boltzman/e_Coulumb
  67     char    aDiffCoeff[100];            //Character array for sprintf
  68
  69     sprintf(aDiffCoeff,"Number of SPD Clusters in x, layer 1, DiffCoeff=%f #sqrt{cm}",
  70             diffCoeff);
  71     TH1D *hclx1 = new TH1D("hclx1",aDiffCoeff,15,0.5,15.5);
  72     sprintf(aDiffCoeff,"Number of SPD Clusters in z, layer 1, DiffCoeff=%f #sqrt{cm}",
  73             diffCoeff);
  74     TH1D *hclz1 = new TH1D("hclz1",aDiffCoeff,5,0.5,5.5);
  75     sprintf(aDiffCoeff,"Number of SPD Clusters in x, layer 2, DiffCoeff=%f #sqrt{cm}",
  76             diffCoeff);
  77     TH1D *hclx2 = new TH1D("hclx2",aDiffCoeff,15,0.5,15.5);
  78     sprintf(aDiffCoeff,"Number of SPD Clusters in z, layer 2, DiffCoeff=%f #sqrt{cm}",
  79             diffCoeff);
  80     TH1D *hclz2 = new TH1D("hclz2",aDiffCoeff,5,0.5,5.5);
  81     // Create Arrays with clusters from:  Data, Ba/Sa Model, old version of
  82     // Dubna
  83     Float_t dataX1[9] = {0.493, 0.493, 0.0273, 0.00617, 0.00112, 0.000257,
  84                          0.000257,  0.000257, 0.000257};
  85                         //clusters from the data in the x-dimension with
  86                         //standard paramaters
  87     Float_t baSaX1[11] = {0.942, 0.0180, 0.0112, 0.00389, 0.00203, 0.000257,
  88                           0.001, 0.000257, 0.000257, 0.000257, 0.001};
  89                         //same for Ba/Sa model
  90     Float_t dubnaX1[7] =  {0.889, 0.0789, 0.0151, 0.00389, 0.001, 0.001,
  91                            0.001};
  92                         //same for old version of Dubna model
  93     Float_t dataX2[9] = {0.482, 0.482, 0.0325, 0.00791, 0.00157, 0.000275,
  94                          0.000275, 0.000275, 0.000275};
  95                         //clusters from the data in the x-dimension with
  96                         //optimized paramaters
  97     Float_t baSaX2[11] = {0.946, 0.0175, 0.0245, 0.00482, 0.001, 0.001,
  98                           0.000275,  0.001, 0.000275, 0.000275, 0.001};
  99                         //same for Ba/Sa model
 100     Float_t dubnaX2[8] = {0.638, 0.325, 0.0275, 0.01, 0.00431, 0.000275,
 101                           0.001, 0.001};
 102                         //same for old version of Dubna model
 103     Float_t dataZ1[4] = {0.889, 0.0624, 0.000257, 0.000257};
 104                         //clusters from the data in the z-dimension with
 105                         //standard paramaters
 106     Float_t baSaZ1[3] = {1., 0.0160, 0.000257}; //same for Ba/Sa model
 107     Float_t dubnaZ1[3] = {0.889, 0.0180, 0.000257};
 108                         //same for old version of Dubna model
 109     Float_t dataZ2[4] = {0.889, 0.0624, 0.000257, 0.000257};
 110                         //clusters from the data in the z-dimension with
 111                         //optimized paramaters
 112     Float_t baSaZ2[4] = {1., 0.0160, 0.00215, 0.000257}; //same for Ba/Sa model
 113     Float_t dubnaZ2[3] = {0.889, 0.0412, 0.000257};
 114                         //same for old version of Dubna model
 115
 116     TH1F *hdataX1  = new TH1F("hdataX1","Number of SPD Clusters in x, layer 1",
 117                               9, 0.5, 9.5);
 118     hdataX1->SetMarkerColor(2);
 119     hdataX1->SetMarkerStyle(20);
 120     hdataX1->SetMarkerSize(0.7);
 121     TH1F *hbaSaX1  = new TH1F("hbaSaX1", "Ba/Sa", 11, 0.5, 11.5);
 122     TH1F *hdubnaX1 = new TH1F("hdubnaX1", "Old Dubna", 7, 0.5, 7.5);
 123     TH1F *hdataX2  = new TH1F("hdataX2","Number of SPD Clusters in x, layer 2",
 124                               9, 0.5, 9.5);
 125     hdataX2->SetMarkerColor(2);
 126     hdataX2->SetMarkerStyle(20);
 127     hdataX2->SetMarkerSize(0.7);
 128     TH1F *hbaSaX2  = new TH1F("hbaSaX2", "Ba/Sa", 11, 0.5, 11.5);
 129     TH1F *hdubnaX2 = new TH1F("hdubnaX2", "Old Dubna", 8, 0.5, 8.5);
 130     TH1F *hdataZ1  = new TH1F("hdataZ1","Number of SPD Clusters in z, layer 1",
 131                               4, 0.5, 4.5);
 132     hdataZ1->SetMarkerColor(2);
 133     hdataZ1->SetMarkerStyle(20);
 134     hdataZ1->SetMarkerSize(0.7);
 135     TH1F *hbaSaZ1  = new TH1F("hbaSaZ1", "Ba/Sa", 3, 0.5, 3.5);
 136     TH1F *hdubnaZ1 = new TH1F("hdubnaZ1", "Old Dubna", 3, 0.5, 3.5);
 137     TH1F *hdataZ2  = new TH1F("hdataZ2","Number of SPD Clusters in z, layer 2",
 138                               4, 0.5, 4.5);
 139     hdataZ2->SetMarkerColor(2);
 140     hdataZ2->SetMarkerStyle(20);
 141     hdataZ2->SetMarkerSize(0.7);
 142     TH1F *hbaSaZ2  = new TH1F("hbaSaZ2", "Ba/Sa", 4, 0.5, 4.5);
 143     TH1F *hdubnaZ2 = new TH1F("hdubnaZ2", "Old Dubna", 3, 0.5, 3.5);
 144
 145     Int_t j = 0; //j will be a counter for the loops to fill the histograms
 146                  //with the values for clusters for the Data, the Ba/Sa model
 147                  //and the old Dubna model
 148     for(j=0; j<9; j++){
 149         hdataX1->SetBinContent((j+1), (Double_t)dataX1[j]);
 150         hdataX1->SetBinError((j+1), 0.00001);
 151     }
 152     for(j=0; j<11; j++) hbaSaX1->Fill((Float_t)j +0.5, baSaX1[j]);
 153     for(j=0; j<7; j++)  hdubnaX1->Fill((Float_t)j +0.5, dubnaX1[j]);
 154     for(j=0; j<9; j++){
 155         hdataX2->SetBinContent((j+1), (Double_t)dataX2[j]);
 156         hdataX2->SetBinError((j+1), 0.00001);
 157     }
 158     for(j=0; j<11; j++) hbaSaX2->Fill((Float_t)j +0.5, baSaX2[j]);
 159     for(j=0; j<8; j++)  hdubnaX2->Fill((Float_t)j +0.5, dubnaX2[j]);
 160     for(j=0; j<4; j++){
 161         hdataZ1->SetBinContent((j+1), (Double_t)dataZ1[j]);
 162         hdataZ1->SetBinError((j+1), 0.00001);
 163     }
 164     for(j=0; j<3; j++)  hbaSaZ1->Fill((Float_t)j +0.5, baSaZ1[j]);
 165     for(j=0; j<3; j++)  hdubnaZ1->Fill((Float_t)j +0.5, dubnaZ1[j]);
 166     for(j=0; j<4; j++){
 167         hdataZ2->SetBinContent((j+1), (Double_t)dataZ2[j]);
 168         hdataZ2->SetBinError((j+1), 0.00001);
 169     }
 170     for(j=0; j<4; j++)  hbaSaZ2->Fill((Float_t)j +0.5, baSaZ2[j]);
 171     for(j=0; j<3; j++)  hdubnaZ2->Fill((Float_t)j +0.5, dubnaZ2[j]);
 172
 173     TTree *tc = 0;
 174     TBranch *spdBranch = 0;
 175     TClonesArray *spdClustcA = new TClonesArray("AliITSRawClusterSPD",1000);
 176     AliITSRawClusterSPD *spdClust = 0;
 177     char tn[20];
 178     Int_t evnt,i,n,nc;
 179     Float_t nclx = 0,nclz = 0;
 180     for(evnt=0;evnt<nevents;evnt++){
 181         sprintf(tn,"TreeC%d",evnt);
 182         tc = (TTree*) file->Get(tn);
 183         spdBranch = tc->GetBranch("ITSClustersSPD");
 184         spdBranch->SetAddress(&spdClustcA);
 185         n = (Int_t) tc->GetEntries();
 186         for(i=0;i<n;i++){
 187             tc->GetEvent(i);
 188             nc = spdClustcA->GetEntries();
 189             if(nc>240) nc = 240;
 190             for(j=0;j<nc;j++){
 191                 spdClust = (AliITSRawClusterSPD*) spdClustcA->At(j);
 192                 nclx = spdClust->NclX();
 193                 nclz = spdClust->NclZ();
 194                 if(spdClust->Module()<80){
 195                     hclx1->Fill(nclx-0.5);
 196                     hclz1->Fill(nclz-0.5);
 197                 }else if(spdClust->Module()<240){
 198                     hclx2->Fill(nclx-0.5);
 199                     hclz2->Fill(nclz-0.5);
 200                 } //endif
 201             } // end for j
 202         } // end for i
 203         spdClustcA->Clear();
 204         delete spdBranch; spdBranch = 0;
 205         delete tc; tc = 0;
 206     } // end for evnt
 207
 208     //Begin Process of normalizing histograms
 209     Double_t integral = (Double_t) hclx1->Integral();
 210     if(integral>0.0) hclx1->Scale(1./integral);
 211     integral = hclz1->Integral();
 212     if(integral>0.0) hclz1->Scale(1./integral);
 213     integral = hclx2->Integral();
 214     if(integral>0.0) hclx2->Scale(1./integral);
 215     integral = hclz2->Integral();
 216     if(integral>0.0) hclz2->Scale(1./integral);
 217
 218     hclx1->SetMinimum(0.000257);
 219     hclx1->SetMaximum(1.01);
 220     hclz1->SetMinimum(0.000274);
 221     hclz1->SetMaximum(1.01);
 222     hclx2->SetMinimum(0.000257);
 223     hclx2->SetMaximum(1.01);
 224     hclz2->SetMinimum(0.000257);
 225     hclz2->SetMaximum(1.01);
 226
 227     sprintf(aDiffCoeff,"SPD Clusters with Diffusion Coefficent=%f",diffCoeff);
 228     TCanvas *cSPDclusters = new TCanvas("cSPDclusters",aDiffCoeff,
 229                                        400,10,600,776);
 230     cSPDclusters->Divide(2, 2);
 231
 232     cSPDclusters->cd(1);
 233     cSPDclusters->SetLogy(1);
 234     hclx1->Draw("hist");
 235     hdataX1->Draw("same e1p");
 236     hbaSaX1->SetLineColor(4);
 237     hbaSaX1->SetLineStyle(2);
 238     hbaSaX1->Draw("same");
 239     hdubnaX1->SetLineColor(3);
 240     hdubnaX1->SetLineStyle(4);
 241     hdubnaX1->Draw("same");
 242     TLegend *l1 = new TLegend(0.55,0.65,0.76,0.82);
 243     l1->AddEntry(hclx1,"New simulation","l");
 244     l1->AddEntry(hdataX1,"Test Beam Results","p");
 245     l1->AddEntry(hbaSaX1,"Bari/Selero simulation","l");
 246     l1->AddEntry(hdubnaX1,"Dubna simulation","l");
 247     l1->Draw();
 248
 249     cSPDclusters->cd(2);
 250     cSPDclusters->SetLogy(1);
 251     hclz1->Draw("hist");
 252     hdataZ1->Draw("same e1p");
 253     hbaSaZ2->SetLineColor(4);
 254     hbaSaZ1->SetLineStyle(2);
 255     hbaSaZ1->Draw("same");
 256     hdubnaZ1->SetLineColor(3);
 257     hdubnaZ1->SetLineStyle(4);
 258     hdubnaZ1->Draw("same");
 259     TLegend *l2 = new TLegend(0.55,0.65,0.76,0.82);
 260     l2->AddEntry(hclz1,"New simulation","l");
 261     l2->AddEntry(hdataZ1,"Test Beam Results","p");
 262     l2->AddEntry(hbaSaZ1,"Bari/Selero simulation","l");
 263     l2->AddEntry(hdubnaZ1,"Dubna simulation","l");
 264     l2->Draw();
 265
 266     cSPDclusters->cd(3);
 267     cSPDclusters->SetLogy(1);
 268     hclx2->Draw("hist");
 269     TLegend *l3 = new TLegend(0.55,0.65,0.76,0.82);
 270     l3->AddEntry(hclx2,"New simulation","l");
 271     l3->Draw();
 272 /*
 273     hdataX2->Draw("e1p");
 274     hbaSaX2->SetLineColor(4);
 275     hbaSaX2->SetLineStyle(2);
 276     hbaSaX2->Draw("same");
 277     hdubnaX2->SetLineColor(3);
 278     hdubnaX2->SetLineStyle(4);
 279     hdubnaX2->Draw("same");
 280 */
 281     cSPDclusters->cd(4);
 282     cSPDclusters->SetLogy(1);
 283     hclz2->Draw("hist");
 284     TLegend *l4 = new TLegend(0.55,0.65,0.76,0.82);
 285     l4->AddEntry(hclz2,"New simulation","l");
 286     l4->Draw();
 287 /*
 288     hdataZ2->Draw("e1p");
 289     hbaSaZ2->SetLineColor(4);
 290     hbaSaZ2->SetLineStyle(2);
 291     hbaSaZ2->Draw("same");
 292     hdubnaZ2->SetLineColor(3);
 293     hdubnaZ2->SetLineStyle(4);
 294     hdubnaZ2->Draw("same");
 295 */
 296     cSPDclusters->Update();
 297
 298     if(gROOT->IsBatch()){
 299         cSPDclusters->Print("SPDClusters.eps");
 300     } // end if gROOT->IsBatch()
 301
 302 }