TOF/AliTOFanalyzeHits.C

   1 Int_t AliTOFanalyzeHits(Int_t numberOfEvents=0, Bool_t drawing=kFALSE)
   2 {
   3
   4   /////////////////////////////////////////////////////////////////////////
   5   //
   6   // Analyzes TOF hits and fills QA-histograms
   7   // and writes the histograms in the TOF_hitsQA.root file
   8   //
   9   // Use case:
  10   // start aliroot
  11   // root [0] .L AliTOFanalyzeHits.C
  12   // root [1] AliTOFanalyzeHits()
  13   //
  14   // By default, it analyzes hits for all the events in the header file
  15   // and does not draw the histograms filled
  16   //
  17   // If you want analyze hits only the 1th event
  18   // you can use the following line:
  19   //
  20   // root[0] .L AliTOFanalyzeHits.C
  21   // root[1] AliTOFanalyzeHits(1)
  22   //
  23   // Updated to the new I/O: C. Zampolli
  24   //
  25   // Report problems to decaro@sa.infn.it
  26   //
  27   /////////////////////////////////////////////////////////////////////////
  28
  29   Int_t rc = 0;
  30
  31   // Define the histograms
  32   // x,y,z, rho, tof, padx, padz, sector, plate, strip, (x vs y)
  33
  34   // hit-map in a plane
  35   TH2F *h2hitMap = new TH2F("h2hitMap","Hit Map (projection on the plane)",2500,-12500.,12500.,800,-400.,400.);
  36   // time of flight distribution for primaries and secondaries
  37   TH1F *htofp    = new TH1F("htofp","Time of Flight (primaries)",800,0.,80.);
  38   TH1F *htofs    = new TH1F("htofs","Time of Flight (secondaries)",800,0.,80.);
  39   // momentum when striking the TOF for primaries and secondaries
  40   TH1F *htofmomp = new TH1F("htofmomp","Momentum at TOF (primaries)",100,0.,10.);
  41   TH1F *htofmoms = new TH1F("htofmoms","Momentum at TOF (secondaries)",100,0.,10.);
  42   // TOF hit volumes
  43   TH1F *hsector  = new TH1F("hsector","Sector",18,0.,18.);
  44   TH1F *hplate   = new TH1F("hplate","Plate ", 5,0., 5.);
  45   TH1F *hstrip   = new TH1F("hstrip","Strip ",20,0.,20.);
  46   TH1F *hpadz    = new TH1F("hpadz","Pad along z ",2,0.,2.);
  47   TH1F *hpadx    = new TH1F("hpadx","Pad along x",48,0.,48.);
  48   // track length when striking the TOF (used by AliTOFT0)
  49   TH1F *htrackLenp= new TH1F("htrackLenp","Track Length on TOF for Primaries",800,0.,800.);
  50
  51   // Histograms added to control the right TOF element numbering:
  52   // it should be increasing with the azimuthal and polar angles
  53
  54   TH2F *hmoduleVStheta     = new TH2F("hmoduleVStheta", "hmoduleVStheta", 180,0.,180.,5,0,5);
  55   TH2F *hsectorVSphi       = new TH2F("hsectorVSphi", "hsectorVSphi", 360,0.,360.,18,0,18);
  56   TH2F *hstripVStheta      = new TH2F("hstripVStheta", "hstripVStheta", 180,0.,180.,20,0,20);
  57   //TH2F *hpadzVStheta       = new TH2F("hpadzVStheta", "hpadzVStheta", 180,0.,180.,2,0,2);
  58   TH2F *hpadxVSphi         = new TH2F("hpadxVSphi", "hpadxVSphi", 360,0.,360.,48,0,48);
  59   TH2F *hpadz2stripVStheta = new TH2F("hpadz2stripVStheta", "hpadz2stripVStheta", 180,0.,180.,40,0,40);
  60   //TH2F *hdzVSpadz2strip    = new TH2F("hdzVSpadz2strip", "hdzVSpadz2strip",40,0,40,70,-3.5,3.5);
  61
  62   // Dynamically link some shared libs
  63   if (gClassTable->GetID("AliRun") < 0) {
  64     gROOT->LoadMacro("loadlibs.C");
  65     loadlibs();
  66   }
  67
  68   if (gAlice)
  69     {
  70       delete AliRunLoader::Instance();
  71       delete gAlice;
  72       gAlice = 0x0;
  73     }
  74
  75   AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
  76   if (!rl)
  77     {
  78       cerr<<"Can't load RunLoader from file! \n";
  79       rc = 1;
  80       return rc;
  81     }
  82
  83   rl->LoadgAlice();
  84   gAlice = rl->GetAliRun();
  85
  86   if (!gAlice)
  87     {
  88       cerr << "<AliTOFanalyzeHits> AliRun object not found on file \n";
  89       rc = 2;
  90       return rc;
  91     }
  92
  93   rl->LoadHeader();
  94
  95   // Get the pointer to the TOF detector
  96   AliLoader *tofl = rl->GetLoader("TOFLoader");
  97   AliTOF *tof = (AliTOF *) gAlice->GetDetector("TOF");
  98   if (tof == 0x0 || tofl == 0x0) {
  99     cerr << "<AliTOFanalyzeHits> Can not find TOF or TOFLoader \n";
 100     rc = 3;
 101     return rc;
 102   }
 103
 104   Int_t countHits = 0;
 105
 106   if (numberOfEvents==0) numberOfEvents=(Int_t)(rl->GetNumberOfEvents());
 107
 108   for (Int_t ievent=0; ievent<numberOfEvents; ievent++) {
 109     printf ("Processing event %d \n", ievent);
 110     rl->GetEvent(ievent);
 111
 112     // Get the pointer Hit tree
 113     tofl->LoadHits();
 114     TTree *hitTree = tofl->TreeH();
 115     tof->SetTreeAddress();
 116     if (!hitTree) {
 117       cout << "<AliTOFanalyzeHits> No TreeH found" << endl;
 118       rc = 4;
 119       return rc;
 120     }
 121
 122     rl->LoadKinematics();
 123     //AliStack* stack = rl->Stack(); // it is not necessary to use the stack!
 124
 125     // Get the number of entries in the hit tree
 126     // (Number of primary particles creating a hit somewhere)
 127     Int_t nTrack = (Int_t) hitTree->GetEntries();
 128     cout << "<AliTOFanalyzeHits> Found " << nTrack
 129          << " primary particles with hits \n";
 130
 131     Int_t nPrimaryOnTof = 0;
 132     Int_t nSecondaryOnTof = 0;
 133     Int_t nelectron  = 0;
 134     Int_t npion      = 0;
 135     Int_t nkaon      = 0;
 136     Int_t nproton    = 0;
 137     Int_t nmuon      = 0;
 138
 139     // Loop through all entries in the tree
 140     for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
 141
 142       tof->ResetHits();
 143       hitTree->GetEvent(iTrack);
 144
 145       // Loop through the TOF hits
 146       AliTOFhitT0 *hit = (AliTOFhitT0 *) tof->FirstHit(-1);
 147       while (hit) {
 148
 149         countHits++;
 150
 151         Float_t x     = hit->X();
 152         Float_t y     = hit->Y();
 153         Float_t z     = hit->Z();
 154         Float_t phiAngle=TMath::Pi() + TMath::ATan2(-y,-x);
 155         Float_t rhoRadius=TMath::Sqrt(x*x+y*y);
 156         Float_t thetaAngle=TMath::Pi() + TMath::ATan2(-rhoRadius,-z);
 157         Float_t dummy=rhoRadius*phiAngle;
 158         h2hitMap->Fill(dummy,z);
 159
 160         phiAngle*=180./TMath::Pi();
 161         thetaAngle*=180./TMath::Pi();
 162
 163         Float_t flightTime = hit->GetTof(); // [s]
 164         flightTime *= 1.e+09; // convert in [ns]
 165         Float_t angle = hit->GetIncA();
 166         Float_t tofmom = hit->GetMom(); // [GeV/c]
 167         Float_t trackLen = hit->GetLen(); // [cm]
 168
 169         // TOF hit volumes
 170         Int_t sector = hit->GetSector(); // range [1-18]
 171         Int_t plate  = hit->GetPlate();  // range [1- 5]
 172         Int_t strip  = hit->GetStrip();  // range [1-20]
 173         Int_t padz   = hit->GetPadz();   // range [1- 2]
 174         Int_t padx   = hit->GetPadx();   // range [1-48]
 175         // it is QA, then I perform QA!
 176
 177         Bool_t isHitBad = (sector<0 || sector>17 ||
 178                             plate<0 || plate>4   ||
 179                              padz<0 || padz>1    ||
 180                              padx<0 || padx>47   ||
 181                            ((strip<0 || strip>14) && plate == 2) ||
 182                            ((strip<0 || strip>18) && (plate == 1 || plate == 3)) ||
 183                            ((strip<0 || strip>19) && (plate == 0 || plate == 4)));
 184
 185         if (isHitBad) {
 186           cout << "<AliTOFanalyzeHits>  strange hit found \n";
 187           cout << "sector = " << sector <<
 188                   " plate = " << plate <<
 189                   " strip = " << strip <<
 190                   " padx = " << padx <<
 191                   " padz = " << padz << endl;
 192           rc = 5;
 193           return rc;
 194         }
 195
 196         hmoduleVStheta->Fill(thetaAngle,plate);
 197         hstripVStheta->Fill(thetaAngle,strip);
 198         hsectorVSphi->Fill(phiAngle,sector);
 199         //hpadzVStheta->Fill(thetaAngle,padx);
 200         hpadxVSphi->Fill(phiAngle,padx);
 201
 202         Float_t dummy2 = 2*strip + padz;
 203         hpadz2stripVStheta->Fill(thetaAngle,dummy2);
 204
 205         /*
 206           Float_t dummy3;
 207           if (hit->GetDz()<=0) dummy3 = hit->GetDz()-1.75;
 208           else dummy3 = hit->GetDz()+1.75;
 209           hdzVSpadz2strip->Fill(dummy2,dummy3);
 210         */
 211
 212         // filling hit volume histos
 213         hsector->Fill(sector);
 214         hplate->Fill(plate);
 215         hstrip->Fill(strip);
 216         hpadx->Fill(padx);
 217         hpadz->Fill(padz);
 218 /*
 219         Int_t track = hit->Track();
 220         TParticle *part = gAlice->GetMCApp()->Particle(track);
 221
 222         // getting MC info for the current track
 223         if (part->GetFirstMother()<0){
 224           Int_t icod = TMath::Abs(part->GetPdgCode());
 225           switch (icod) {
 226           case 211:
 227             npion++;
 228             break ;
 229           case 321:
 230             nkaon++;
 231             break ;
 232           case 2212:
 233             nproton++;
 234             break ;
 235           case 11:
 236             nelectron++;
 237             break ;
 238           case 13:
 239             nmuon++;
 240             break ;
 241           }
 242           htofp->Fill(flightTime);
 243           htofmomp->Fill(tofmom);
 244           htrackLenp->Fill(trackLen);
 245         } else {
 246           htofs->Fill(flightTime);
 247           htofmoms->Fill(tofmom);
 248         }
 249 */
 250         // go to next hit
 251         hit = (AliTOFhitT0 *) tof->NextHit();
 252
 253       }
 254
 255     }
 256
 257     tofl->UnloadHits();
 258     rl->UnloadKinematics();
 259
 260     cout << "<AliTOFanalyzeHits> Found " << countHits << " hits in total \n";
 261     cout << npion     << " primary pions reached the TOF detector \n";
 262     cout << nkaon     << " primary kaons reached the TOF detector \n";
 263     cout << nproton   << " primary protons reached the TOF detector \n";
 264     cout << nelectron << " primary electrons reached the TOF detector \n";
 265     cout << nmuon     << " primary muons reached the TOF detector \n";
 266
 267   }
 268
 269   cout << "hpadx->GetEntries() = " << hpadx->GetEntries() << endl;
 270
 271   rl->UnloadHeader();
 272   rl->UnloadgAlice();
 273
 274   if (drawing) {
 275     TCanvas *cHits = new TCanvas("cHits","AliTOFanalyzeHits hit volumes",50,50,900,900);
 276     cHits->Divide(3,2);
 277     cHits->cd(1);
 278     hsector->Draw();
 279     cHits->cd(2);
 280     hplate->Draw();
 281     cHits->cd(3);
 282     hstrip->Draw();
 283     cHits->cd(4);
 284     hpadz->Draw();
 285     cHits->cd(5);
 286     hpadx->Draw();
 287
 288     TCanvas *chitmap = new TCanvas("chitmap","AliTOFanalyzeHits Hit Map",50,50,600,600);
 289     chitmap->cd();
 290     h2hitMap->Draw();
 291
 292     TCanvas *ctrackLen = new TCanvas("ctrackLen","AliTOFanalyzeHits Track Length for primaries on TOF",50,50,400,400);
 293     ctrackLen->cd();
 294     htrackLenp->Draw();
 295
 296     TCanvas *ctofmom = new TCanvas("ctofmom","AliTOFanalyzeHits flight times",50,50,700,700);
 297     ctofmom->Divide(2,2);
 298     ctofmom->cd(1);
 299     gPad->SetLogy();
 300     htofp->Draw();
 301     ctofmom->cd(2);
 302     gPad->SetLogy();
 303     htofs->Draw();
 304     ctofmom->cd(3);
 305     gPad->SetLogy();
 306     htofmomp->Draw();
 307     ctofmom->cd(4);
 308     gPad->SetLogy();
 309     htofmoms->Draw();
 310   }
 311
 312   // save histos into file TOF_hitsQA.root
 313   TFile *fout = new TFile("TOF_hitsQA.root","RECREATE");
 314   h2hitMap->Write();
 315   htofp->Write();
 316   htofs->Write();
 317   htofmomp->Write();
 318   htofmoms->Write();
 319   hsector->Write();
 320   hplate->Write();
 321   hstrip->Write();
 322   hpadz->Write();
 323   hpadx->Write();
 324   htrackLenp->Write();
 325
 326   hmoduleVStheta->Write();
 327   hsectorVSphi->Write();
 328   hstripVStheta->Write();
 329   //hpadzVStheta->Write();
 330   hpadxVSphi->Write();
 331   hpadz2stripVStheta->Write();
 332   //hdzVSpadz2strip->Write();
 333
 334   fout->Close();
 335
 336   cout << " Finished AliTOFanalizeHits \n";
 337
 338   if (gAlice)
 339     {
 340       delete AliRunLoader::Instance();
 341       delete gAlice;
 342       gAlice = 0x0;
 343     }
 344
 345   return rc;
 346
 347 }