[u/mrichter/AliRoot.git] / TOF / AliTOFanalyzeHits.C

Int_t AliTOFanalyzeHits(Int_t nevents = 1, Bool_t drawing = kFALSE)
{
  //
  // Analyzes the hits and fills QA-histograms 
  // Use case:
  // start aliroot
  // root [0] .L AliTOFanalyzeHits.C
  // root [1] AliTOFanalyzeHits()
  //
  // Updated to the new I/O by: A. De Caro, C. Zampolli

  Int_t rc = 0;
  
  // Define the histograms
  // x,y,z, rho, tof, padx, padz, sector, plate, strip, (x vs y)

  // hit-map in a plane
  TH2F *h2hitMap = new TH2F("h2hitMap","Hit Map (projection on the plane)",2500,-12500.,12500.,800,-400.,400.);
  // time of flight distribution for primaries and secondaries
  TH1F *htofp    = new TH1F("htofp","Time of Flight (primaries)",800,0.,80.);
  TH1F *htofs    = new TH1F("htofs","Time of Flight (secondaries)",800,0.,80.);
  // momentum when striking the TOF for primaries and secondaries
  TH1F *htofmomp = new TH1F("htofmomp","Momentum at TOF (primaries)",100,0.,10.);
  TH1F *htofmoms = new TH1F("htofmoms","Momentum at TOF (secondaries)",100,0.,10.);
  // TOF hit volumes
  TH1F *hsector  = new TH1F("hsector","Sector",20,0.,20.);
  TH1F *hplate   = new TH1F("hplate","Plate ", 6,0., 6.);
  TH1F *hstrip   = new TH1F("hstrip","Strip ",25,0.,25.);
  TH1F *hpadz    = new TH1F("hpadz","Pad along z ",3,0.,3.);
  TH1F *hpadx    = new TH1F("hpadx","Pad along x",50,0.,50.);
  // track length when striking the TOF (used by AliTOFT0)
  TH1F *htrackLenp= new TH1F("htrackLenp","Track Length on TOF for Primaries",800,0.,800.);

  // Histograms added to control the right TOF element numbering:
  // it should be increasing with the azimuthal and polar angles

  TH2F *hmoduleVStheta = new TH2F("hmoduleVStheta", "hmoduleVStheta", 180,0.,180.,6,0,6);
  TH2F *hsectorVSphi   = new TH2F("hsectorVSphi", "hsectorVSphi", 360,0.,360.,19,0,19);
  TH2F *hstripVStheta   = new TH2F("hstripVStheta", "hstripVStheta", 180,0.,180.,25,0,25);
  TH2F *hpadzVStheta   = new TH2F("hpadzVStheta", "hpadzVStheta", 180,0.,180.,3,0,3);
  TH2F *hpadxVSphi     = new TH2F("hpadxVSphi", "hpadxVSphi", 360,0.,360.,49,0,49);

  // Dynamically link some shared libs
  if (gClassTable->GetID("AliRun") < 0) {
    gROOT->LoadMacro("loadlibs.C");
    loadlibs();
  }

  if (gAlice)
    {
      delete gAlice->GetRunLoader();
      delete gAlice;
      gAlice = 0x0;
    }
  
  AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::fgkDefaultEventFolderName,"read");
  if (!rl)
    {
      cerr<<"Can't load RunLoader from file"<<"!\n";
      rc = 1;
      return rc;
    }
  
  rl->LoadgAlice();
  gAlice = rl->GetAliRun();

  if (!gAlice)
    {
      cerr << "<AliTOFanalyzeHits> AliRun object not found on file\n ";
      rc = 2;
      return rc;
    }

  // Get the pointer to the TOF detector
  AliTOF *tof = (AliTOF *) gAlice->GetDetector("TOF");
  AliLoader *tofl = rl->GetLoader("TOFLoader");
  if (tof == 0x0 || tofl == 0x0) {
    cerr << "<AliTOFanalyzeHits> Can not find TOF or TOFLoader\n";
    rc = 3;
    return rc;
  }

  Int_t countHits = 0;

  rl->LoadHeader();

  for (Int_t ievent=0; ievent<nevents; ievent++) {
    printf ("Processing event %d \n", ievent);
    rl->GetEvent(ievent);

    // Get the pointer Hit tree
    tofl->LoadHits();
    TTree *hitTree = tofl->TreeH();
    tof->SetTreeAddress();
    if (!hitTree) {
      cout << "<AliTOFanalyzeHits> No TreeH found" << endl;
      rc = 4;
      return rc;
    }

    rl->LoadKinematics();
    //AliStack* stack = rl->Stack(); // it is not necessary to use the stack!

    // Get the number of entries in the hit tree
    // (Number of primary particles creating a hit somewhere)
    Int_t nTrack = (Int_t) hitTree->GetEntries();
    cout << "<AliTOFanalyzeHits> Found " << nTrack 
         << " primary particles with hits" << endl;

    Int_t nPrimaryOnTof = 0;
    Int_t nSecondaryOnTof = 0;
    Int_t nelectron  = 0;
    Int_t npion      = 0;
    Int_t nkaon      = 0;
    Int_t nproton    = 0;    
    Int_t nmuon      = 0;
  
    // Loop through all entries in the tree
    for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {

      tof->ResetHits();
      hitTree->GetEvent(iTrack);

      // Loop through the TOF hits  
      AliTOFhitT0 *hit = (AliTOFhitT0 *) tof->FirstHit(-1);
      while (hit) {

	countHits++;
	
        Float_t x     = hit->X();
        Float_t y     = hit->Y();
        Float_t z     = hit->Z();
        Float_t phiAngle=TMath::Pi() + TMath::ATan2(-y,-x);
        Float_t rhoRadius=TMath::Sqrt(x*x+y*y);
        Float_t thetaAngle=TMath::Pi() + TMath::ATan2(-rhoRadius,-z);
        Float_t dummy=rhoRadius*phiAngle;
        h2hitMap->Fill(dummy,z);

        phiAngle*=180./TMath::Pi();
        thetaAngle*=180./TMath::Pi();

        Float_t flightTime = hit->GetTof(); // [s]
        flightTime *= 1.e+09; // convert in [ns]
        Float_t angle = hit->GetIncA();
        Float_t tofmom = hit->GetMom(); // [GeV/c]
        Float_t trackLen = hit->GetLen(); // [cm]

        // TOF hit volumes
        Int_t sector = hit->GetSector(); // range [1-18]
        Int_t plate  = hit->GetPlate();  // range [1- 5]
        Int_t strip  = hit->GetStrip();  // range [1-20]
        Int_t padz   = hit->GetPadz();   // range [1- 2]
        Int_t padx   = hit->GetPadx();   // range [1-48]
        // it is QA, then I perform QA!
        Bool_t isHitBad = (sector<1 || sector>18 || plate<1 || plate >5 || padz<1 || padz>2 || padx<1 || padx>48);

        if (isHitBad) {
          cout << "<AliTOFanalyzeHits>  strange hit found" << endl;
	  rc = 5;
	  return rc;
        }

        hmoduleVStheta->Fill(thetaAngle,plate);
	hstripVStheta->Fill(thetaAngle,strip);
        hsectorVSphi->Fill(phiAngle,sector);
        hpadxVSphi->Fill(phiAngle,padx);
        hpadzVStheta->Fill(thetaAngle,padz);

        // filling hit volume histos
        hsector->Fill(sector);
        hplate->Fill(plate);
        hstrip->Fill(strip);
        hpadx->Fill(padx);
        hpadz->Fill(padz);

        Int_t track = hit->Track();
        TParticle *part = gAlice->Particle(track);

        // getting MC info for the current track
        if (part->GetFirstMother()<0){
          Int_t icod = TMath::Abs(part->GetPdgCode());
 	  switch (icod) {
	  case 211:
	    npion++;
	    break ;
	  case 321:
	    nkaon++;
	    break ;
	  case 2212:
	    nproton++;
	    break ;
	  case 11:
	    nelectron++;
	    break ;
	  case 13:
	    nmuon++;
	    break ;
	  }
	  htofp->Fill(flightTime);
	  htofmomp->Fill(tofmom);
          htrackLenp->Fill(trackLen);
        } else {
	  htofs->Fill(flightTime);
	  htofmoms->Fill(tofmom);
        }

        // go to next hit
        hit = (AliTOFhitT0 *) tof->NextHit();

      }

    }

    tofl->UnloadHits();
    rl->UnloadKinematics();

    cout << "<AliTOFanalyzeHits> Found " << countHits << " hits in total" << endl;
    cout << npion     << " primary pions reached the TOF detector"     << endl;
    cout << nkaon     << " primary kaons reached the TOF detector"     << endl;
    cout << nproton   << " primary protons reached the TOF detector"   << endl;
    cout << nelectron << " primary electrons reached the TOF detector" << endl;
    cout << nmuon     << " primary muons reached the TOF detector"     << endl;

  }

  rl->UnloadHeader();
  rl->UnloadgAlice();

  if (drawing) {  
    TCanvas *cHits = new TCanvas("cHits","AliTOFanalyzeHits hit volumes",50,50,900,900);
    cHits->Divide(3,2);
    cHits->cd(1);
    hsector->Draw();
    cHits->cd(2);
    hplate->Draw();
    cHits->cd(3);
    hstrip->Draw();
    cHits->cd(4);
    hpadz->Draw();
    cHits->cd(5);
    hpadx->Draw();
    
    TCanvas *chitmap = new TCanvas("chitmap","AliTOFanalyzeHits Hit Map",50,50,600,600);
    chitmap->cd();
    h2hitMap->Draw();
    
    TCanvas *ctrackLen = new TCanvas("ctrackLen","AliTOFanalyzeHits Track Length for primaries on TOF",50,50,400,400);
    ctrackLen->cd();
    htrackLenp->Draw();
    
    TCanvas *ctofmom = new TCanvas("ctofmom","AliTOFanalyzeHits flight times",50,50,700,700);
    ctofmom->Divide(2,2);
    ctofmom->cd(1);
    gPad->SetLogy();
    htofp->Draw();
    ctofmom->cd(2);
    gPad->SetLogy();
    htofs->Draw();
    ctofmom->cd(3);
    gPad->SetLogy();
    htofmomp->Draw();
    ctofmom->cd(4);
    gPad->SetLogy();
    htofmoms->Draw();
  }
  
  // save histos into file TOF_hitsQA.root
  TFile *fout = new TFile("TOF_hitsQA.root","RECREATE");
  h2hitMap->Write();
  htofp->Write();
  htofs->Write();
  htofmomp->Write();
  htofmoms->Write();
  hsector->Write();
  hplate->Write();
  hstrip->Write();
  hpadz->Write();
  hpadx->Write();
  htrackLenp->Write();

  hmoduleVStheta->Write();
  hsectorVSphi->Write();
  hstripVStheta->Write();
  hpadzVStheta->Write();
  hpadxVSphi->Write();

  fout->Close(); 

  cout << " Finished AliTOFanalizeHits \n";

  if (gAlice)
    {
      delete gAlice->GetRunLoader();
      delete gAlice;
      gAlice = 0x0;
    }
  
  return rc;

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