--- /dev/null
+Int_t AliTOFanalyzeHits()
+{
+ //
+ // Analyzes the hits and fills QA-histograms
+ //
+
+ Int_t rc = 0;
+
+ if (!gAlice) {
+ cout << "<AliTOFanalyzeHits> No AliRun object found" << endl;
+ rc = 1;
+ return rc;
+ }
+ gAlice->GetEvent(0);
+
+ // Get the pointer to the TOF detector
+ AliTOF *tof = (AliTOF *) gAlice->GetDetector("TOF");
+ if (!tof) {
+ cout << "<AliTOFanalyzeHits> No TOF detector found" << endl;
+ rc = 2;
+ return rc;
+ }
+
+ // 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.);
+
+ // Get the pointer hit tree
+ TTree *hitTree = gAlice->TreeH();
+ if (!hitTree) {
+ cout << "<AliTOFanalyzeHits> No hit tree found" << endl;
+ rc = 4;
+ return rc;
+ }
+
+ Int_t countHits = 0;
+ Int_t nBytes = 0;
+
+ // 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++) {
+
+ gAlice->ResetHits();
+ nBytes += hitTree->GetEvent(iTrack);
+
+
+ // Loop through the TOF hits
+ Int_t iHit = 0;
+ AliTOFhitT0 *hit = (AliTOFhitT0 *) tof->FirstHit(-1);
+ while (hit) {
+
+ countHits++;
+ iHit++;
+
+ Float_t x = hit->X();
+ Float_t y = hit->Y();
+ Float_t z = hit->Z();
+ Float_t circleLen=TMath::Sqrt(x*x+y*y)*TMath::ATan2(y,x);
+ h2hitMap->Fill(circleLen,z);
+
+ 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 = 3;
+ return rc;
+ }
+ // 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();
+
+ }
+
+ }
+
+ 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;
+
+
+ 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,400,400);
+ 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();
+ fout->Close();
+
+ return rc;
+
+}
+
--- /dev/null
+static Int_t eventsPerRun = 100;
+void Config()
+{
+ // 7-DEC-2000 09:00
+ // Switch on Transition Radiation simulation. 6/12/00 18:00
+ // iZDC=1 7/12/00 09:00
+ // ThetaRange is (0., 180.). It was (0.28,179.72) 7/12/00 09:00
+ // Theta range given through pseudorapidity limits 22/6/2001
+
+ // Set Random Number seed
+ // gRandom->SetSeed(12345);
+
+ new AliGeant3("C++ Interface to Geant3");
+
+ if (!gSystem->Getenv("CONFIG_FILE"))
+ {
+ TFile *rootfile = new TFile("TOF_test.root", "recreate");
+
+ rootfile->SetCompressionLevel(2);
+ }
+
+ TGeant3 *geant3 = (TGeant3 *) gMC;
+
+ //
+ // Set External decayer
+ AliDecayer *decayer = new AliDecayerPythia();
+
+ decayer->SetForceDecay(kAll);
+ decayer->Init();
+ gMC->SetExternalDecayer(decayer);
+ //
+ //
+ //=======================================================================
+ // ******* GEANT STEERING parameters FOR ALICE SIMULATION *******
+ geant3->SetTRIG(1); //Number of events to be processed
+ geant3->SetSWIT(4, 10);
+ geant3->SetDEBU(0, 0, 1);
+ //geant3->SetSWIT(2,2);
+ geant3->SetDCAY(1);
+ geant3->SetPAIR(1);
+ geant3->SetCOMP(1);
+ geant3->SetPHOT(1);
+ geant3->SetPFIS(0);
+ geant3->SetDRAY(0);
+ geant3->SetANNI(1);
+ geant3->SetBREM(1);
+ geant3->SetMUNU(1);
+ geant3->SetCKOV(1);
+ geant3->SetHADR(1); //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3)
+ geant3->SetLOSS(2);
+ geant3->SetMULS(1);
+ geant3->SetRAYL(1);
+ geant3->SetAUTO(1); //Select automatic STMIN etc... calc. (AUTO 1) or manual (AUTO 0)
+ geant3->SetABAN(0); //Restore 3.16 behaviour for abandoned tracks
+ geant3->SetOPTI(2); //Select optimisation level for GEANT geometry searches (0,1,2)
+ geant3->SetERAN(5.e-7);
+
+ Float_t cut = 1.e-3; // 1MeV cut by default
+ Float_t tofmax = 1.e10;
+
+ // GAM ELEC NHAD CHAD MUON EBREM MUHAB EDEL MUDEL MUPA TOFMAX
+ geant3->SetCUTS(cut, cut, cut, cut, cut, cut, cut, cut, cut, cut,
+ tofmax);
+
+ // #### AliGenHIJING generation######################################
+ //
+ AliGenHijing *gener = new AliGenHijing(1);
+ gener->SetEnergyCMS(5500);
+ gener->SetReferenceFrame("CMS");
+ gener->SetProjectile("A", 208,82);
+ gener->SetTarget("A", 208,82);
+ gener->SetImpactParameterRange(0.,3.);
+ gener->SetJetQuenching(1);
+ gener->SetShadowing(1);
+ gener->SetDecaysOff(0);
+ gener->SetSelectAll(0);
+ gener->SetOrigin(0.,0.,0.);
+ gener->SetSpectators(0);
+ gener->SetThetaRange(85.,95.);
+ gener->SetTrackingFlag(1);
+ gener->Init();
+
+ //
+ // Activate this line if you want the vertex smearing to happen
+ // track by track
+ //
+ //gener->SetVertexSmear(perTrack);
+
+ gAlice->SetField(-999, 2); //Specify maximum magnetic field in Tesla (neg. ==> default field)
+ //gAlice->SetField(-999,2,2.);//Scale factor=2 ==> field=0.4T
+
+ Int_t iABSO = 1;
+ Int_t iCASTOR = 1;
+ Int_t iDIPO = 1;
+ Int_t iFMD = 1;
+ Int_t iFRAME = 1;
+ Int_t iHALL = 1;
+ Int_t iITS = 1;
+ Int_t iMAG = 1;
+ Int_t iMUON = 1;
+ Int_t iPHOS = 1;
+ Int_t iPIPE = 1;
+ Int_t iPMD = 1;
+ Int_t iRICH = 1;
+ Int_t iSHIL = 1;
+ Int_t iSTART = 1;
+ Int_t iTOF = 1;
+ Int_t iTPC = 1;
+ Int_t iTRD = 1;
+ Int_t iZDC = 1;
+ Int_t iEMCAL = 0;
+
+ //=================== Alice BODY parameters =============================
+ AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
+
+
+ if (iMAG)
+ {
+ //=================== MAG parameters ============================
+ // --- Start with Magnet since detector layouts may be depending ---
+ // --- on the selected Magnet dimensions ---
+ AliMAG *MAG = new AliMAG("MAG", "Magnet");
+ }
+
+
+ if (iABSO)
+ {
+ //=================== ABSO parameters ============================
+ AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
+ }
+
+ if (iDIPO)
+ {
+ //=================== DIPO parameters ============================
+
+ AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
+ }
+
+ if (iHALL)
+ {
+ //=================== HALL parameters ============================
+
+ AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
+ }
+
+
+ if (iFRAME)
+ {
+ //=================== FRAME parameters ============================
+
+ AliFRAME *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
+
+ }
+
+ if (iSHIL)
+ {
+ //=================== SHIL parameters ============================
+
+ AliSHIL *SHIL = new AliSHILvF("SHIL", "Shielding");
+ }
+
+
+ if (iPIPE)
+ {
+ //=================== PIPE parameters ============================
+
+ AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
+ }
+
+ if(iITS) {
+
+//=================== ITS parameters ============================
+ //
+ // As the innermost detector in ALICE, the Inner Tracking System "impacts" on
+ // almost all other detectors. This involves the fact that the ITS geometry
+ // still has several options to be followed in parallel in order to determine
+ // the best set-up which minimizes the induced background. All the geometries
+ // available to date are described in the following. Read carefully the comments
+ // and use the default version (the only one uncommented) unless you are making
+ // comparisons and you know what you are doing. In this case just uncomment the
+ // ITS geometry you want to use and run Aliroot.
+ //
+ // Detailed geometries:
+ //
+ //
+ //AliITS *ITS = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
+ //
+ //AliITS *ITS = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
+ //
+ AliITSvPPRasymm *ITS = new AliITSvPPRasymm("ITS","New ITS PPR detailed version with asymmetric services");
+ ITS->SetMinorVersion(2); // don't touch this parameter if you're not an ITS developer
+ ITS->SetReadDet(kFALSE); // don't touch this parameter if you're not an ITS developer
+ // ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det"); // don't touch this parameter if you're not an ITS developer
+ ITS->SetThicknessDet1(200.); // detector thickness on layer 1 must be in the range [100,300]
+ ITS->SetThicknessDet2(200.); // detector thickness on layer 2 must be in the range [100,300]
+ ITS->SetThicknessChip1(200.); // chip thickness on layer 1 must be in the range [150,300]
+ ITS->SetThicknessChip2(200.); // chip thickness on layer 2 must be in the range [150,300]
+ ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ ITS->SetCoolingFluid(1); // 1 --> water ; 0 --> freon
+ //
+ //AliITSvPPRsymm *ITS = new AliITSvPPRsymm("ITS","New ITS PPR detailed version with symmetric services");
+ //ITS->SetMinorVersion(2); // don't touch this parameter if you're not an ITS developer
+ //ITS->SetReadDet(kFALSE); // don't touch this parameter if you're not an ITS developer
+ //ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRsymm2.det"); // don't touch this parameter if you're not an ITS developer
+ //ITS->SetThicknessDet1(200.); // detector thickness on layer 1 must be in the range [100,300]
+ //ITS->SetThicknessDet2(200.); // detector thickness on layer 2 must be in the range [100,300]
+ //ITS->SetThicknessChip1(200.); // chip thickness on layer 1 must be in the range [150,300]
+ //ITS->SetThicknessChip2(200.); // chip thickness on layer 2 must be in the range [150,300]
+ //ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ //ITS->SetCoolingFluid(1); // 1 --> water ; 0 --> freon
+ //
+ //
+ // Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful
+ // for reconstruction !):
+ //
+ //
+ //AliITSvPPRcoarseasymm *ITS = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services");
+ //ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ //ITS->SetSupportMaterial(0); // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
+ //
+ //AliITS *ITS = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services");
+ //ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ //ITS->SetSupportMaterial(0); // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
+ //
+ //
+ //
+ // Geant3 <-> EUCLID conversion
+ // ============================
+ //
+ // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
+ // media to two ASCII files (called by default ITSgeometry.euc and
+ // ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
+ // The default (=0) means that you dont want to use this facility.
+ //
+ ITS->SetEUCLID(0);
+ }
+
+
+ if (iTPC)
+ {
+ //============================ TPC parameters ================================
+ // --- This allows the user to specify sectors for the SLOW (TPC geometry 2)
+ // --- Simulator. SecAL (SecAU) <0 means that ALL lower (upper)
+ // --- sectors are specified, any value other than that requires at least one
+ // --- sector (lower or upper)to be specified!
+ // --- Reminder: sectors 1-24 are lower sectors (1-12 -> z>0, 13-24 -> z<0)
+ // --- sectors 25-72 are the upper ones (25-48 -> z>0, 49-72 -> z<0)
+ // --- SecLows - number of lower sectors specified (up to 6)
+ // --- SecUps - number of upper sectors specified (up to 12)
+ // --- Sens - sensitive strips for the Slow Simulator !!!
+ // --- This does NOT work if all S or L-sectors are specified, i.e.
+ // --- if SecAL or SecAU < 0
+ //
+ //
+ //-----------------------------------------------------------------------------
+
+ // gROOT->LoadMacro("SetTPCParam.C");
+ // AliTPCParam *param = SetTPCParam();
+ //AliTPC *TPC = new AliTPCv2("TPC", "Default");
+ AliTPC *TPC = new AliTPCv1("TPC","Default"); //fast simulation (Pierella)
+
+ // All sectors included
+ TPC->SetSecAL(-1);
+ TPC->SetSecAU(-1);
+
+ }
+
+ if (iTOF)
+ {
+ //=================== TOF parameters ============================
+ //AliTOF *TOF = new AliTOFv2("TOF", "normal TOF");
+ AliTOF *TOF = new AliTOFv4T0("TOF", "normal TOF"); // default for time zero determination
+ }
+
+ if (iRICH)
+ {
+ //=================== RICH parameters ===========================
+ AliRICH *RICH = new AliRICHv1("RICH", "normal RICH");
+
+ }
+
+
+ if (iZDC)
+ {
+ //=================== ZDC parameters ============================
+
+ AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC");
+ }
+
+ if (iCASTOR)
+ {
+ //=================== CASTOR parameters ============================
+
+ AliCASTOR *CASTOR = new AliCASTORv1("CASTOR", "normal CASTOR");
+ }
+
+ if (iTRD)
+ {
+ //=================== TRD parameters ============================
+
+ //AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
+ AliTRD *TRD = new AliTRDv0("TRD","TRD fast simulator");// Pierella
+
+ // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 90% Xe + 10% CO2)
+ TRD->SetGasMix(1);
+
+ // With hole in front of PHOS
+ TRD->SetPHOShole();
+ // With hole in front of RICH
+ TRD->SetRICHhole();
+ // Switch on TR
+ AliTRDsim *TRDsim = TRD->CreateTR();
+ }
+
+ if (iFMD)
+ {
+ //=================== FMD parameters ============================
+
+ AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
+ FMD->SetRingsSi1(256);
+ FMD->SetRingsSi2(64);
+ FMD->SetSectorsSi1(20);
+ FMD->SetSectorsSi2(24);
+ }
+
+ if (iMUON)
+ {
+ //=================== MUON parameters ===========================
+
+ AliMUON *MUON = new AliMUONv1("MUON", "default");
+ }
+ //=================== PHOS parameters ===========================
+
+ if (iPHOS)
+ {
+ AliPHOS *PHOS = new AliPHOSv1("PHOS", "GPS2");
+ }
+
+
+ if (iPMD)
+ {
+ //=================== PMD parameters ============================
+
+ AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
+
+ PMD->SetPAR(1., 1., 0.8, 0.02);
+ PMD->SetIN(6., 18., -580., 27., 27.);
+ PMD->SetGEO(0.0, 0.2, 4.);
+ PMD->SetPadSize(0.8, 1.0, 1.0, 1.5);
+
+ }
+
+ if (iEMCAL && !iRICH)
+ {
+ //=================== EMCAL parameters ============================
+ AliEMCAL *EMCAL = new AliEMCALv1("EMCAL", "EMCALArch1a");
+ }
+
+ if (iSTART)
+ {
+ //=================== START parameters ============================
+ AliSTART *START = new AliSTARTv1("START", "START Detector");
+ }
+
+
+}
+
+Float_t EtaToTheta(Float_t arg){
+ return (180./TMath::Pi())*2.*atan(exp(-arg));
+}
--- /dev/null
+Int_t AliTOFtest()
+{
+ //
+ // Test macro for the TOF code
+ // report bug to Fabrizio.Pierella@cern.ch
+ // Use case:
+ // start aliroot
+ // root [0] .L AliTOFtest.C
+ // root [1] AliTOFtest()
+
+ Int_t rc = 0;
+
+ // Initialize the test setup
+
+ //gAlice->Init("$(ALICE_ROOT)/TOF/AliTOFconfig.C");
+ gAlice->Init("AliTOFconfig.C");
+
+ // Run one central Hijing event and create the hits (time required:
+ // some minuts)
+
+ gAlice->SetDebug(2);
+ gAlice->Run(1);
+
+ if (gAlice) delete gAlice;
+ TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TOF_test.root");
+ gAlice = (AliRun *) file->Get("gAlice");
+
+ // Analyze the TOF hits
+ if (rc = AliTOFanalyzeHits()) return rc;
+
+ return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTOFanalyzeHits()
+{
+ //
+ // Analyzes the hits and fills QA-histograms
+ //
+
+ Int_t rc = 0;
+
+ if (!gAlice) {
+ cout << "<AliTOFanalyzeHits> No AliRun object found" << endl;
+ rc = 1;
+ return rc;
+ }
+ gAlice->GetEvent(0);
+
+ // Get the pointer to the TOF detector
+ AliTOF *tof = (AliTOF *) gAlice->GetDetector("TOF");
+ if (!tof) {
+ cout << "<AliTOFanalyzeHits> No TOF detector found" << endl;
+ rc = 2;
+ return rc;
+ }
+
+ // 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.);
+
+ // Get the pointer hit tree
+ TTree *hitTree = gAlice->TreeH();
+ if (!hitTree) {
+ cout << "<AliTOFanalyzeHits> No hit tree found" << endl;
+ rc = 4;
+ return rc;
+ }
+
+ Int_t countHits = 0;
+ Int_t nBytes = 0;
+
+ // 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++) {
+
+ gAlice->ResetHits();
+ nBytes += hitTree->GetEvent(iTrack);
+
+
+ // Loop through the TOF hits
+ Int_t iHit = 0;
+ AliTOFhitT0 *hit = (AliTOFhitT0 *) tof->FirstHit(-1);
+ while (hit) {
+
+ countHits++;
+ iHit++;
+
+ Float_t x = hit->X();
+ Float_t y = hit->Y();
+ Float_t z = hit->Z();
+ Float_t circleLen=TMath::Sqrt(x*x+y*y)*TMath::ATan2(y,x);
+ h2hitMap->Fill(circleLen,z);
+
+ 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 = 3;
+ return rc;
+ }
+ // 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();
+
+ }
+
+ }
+
+ 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;
+
+
+ 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();
+ fout->Close();
+
+ return rc;
+
+}
//
// Calculate the t0 distribution
- //
+ // selecting 15 tracks
// (primary particles reaching TOF)
- //
-
-
+ // NB: see the description of the analogous AliTOFT0 class
+ // NB II: apply this macro to a realistic Pb-Pb event (e.g. Hijing event)
+ // Use case
+ // - start root
+ // root [0] .L tzero15tracksopt.C
+ // // exec the macro with the default time resolution 120 ps
+ // root [1] tzero15tracksopt("hij25evCentralTOFv4T0Set10.04T.root")
+
+ // // exec the macro with 150 ps global time resolution
+ // // NB: time resolution has to be given in [s]
+ // root [1] tzero15tracksopt("hij25evCentralTOFv4T0Set10.04T.root", 1.5e-10))
+
+
// Dynamically link some shared libs
if (gClassTable->GetID("AliRun") < 0) {
gROOT->LoadMacro("loadlibs.C");
Bool_t isgoodpart=(abspdg==211 || abspdg==2212 || abspdg==321);
time*=1.e+9; // tof given in nanoseconds
- if (particle->GetFirstMother() < 0 && isgoodpart && mom<=1.7 && mom>=1.25){
+ if (particle->GetFirstMother() < 0 && isgoodpart && mom<=1.75 && mom>=1.25){
selected+=1;
istop=selected;
selected=0;
cout << "starting t0 calculation for current set" << endl;
for (Int_t i1=0; i1<3;i1++) {
+ beta[0]=momentum[0]/TMath::Sqrt(massarray[i1]*massarray[i1]+momentum[0]*momentum[0]);
for (Int_t i2=0; i2<3;i2++) {
+ beta[1]=momentum[1]/TMath::Sqrt(massarray[i2]*massarray[i2]+momentum[1]*momentum[1]);
for (Int_t i3=0; i3<3;i3++) {
+ beta[2]=momentum[2]/TMath::Sqrt(massarray[i3]*massarray[i3]+momentum[2]*momentum[2]);
for (Int_t i4=0; i4<3;i4++) {
+ beta[3]=momentum[3]/TMath::Sqrt(massarray[i4]*massarray[i4]+momentum[3]*momentum[3]);
for (Int_t i5=0; i5<3;i5++) {
+ beta[4]=momentum[4]/TMath::Sqrt(massarray[i5]*massarray[i5]+momentum[4]*momentum[4]);
for (Int_t i6=0; i6<3;i6++) {
+ beta[5]=momentum[5]/TMath::Sqrt(massarray[i6]*massarray[i6]+momentum[5]*momentum[5]);
for (Int_t i7=0; i7<3;i7++) {
+ beta[6]=momentum[6]/TMath::Sqrt(massarray[i7]*massarray[i7]+momentum[6]*momentum[6]);
for (Int_t i8=0; i8<3;i8++) {
+ beta[7]=momentum[7]/TMath::Sqrt(massarray[i8]*massarray[i8]+momentum[7]*momentum[7]);
for (Int_t i9=0; i9<3;i9++) {
+ beta[8]=momentum[8]/TMath::Sqrt(massarray[i9]*massarray[i9]+momentum[8]*momentum[8]);
for (Int_t i10=0; i10<3;i10++) {
+ beta[9]=momentum[9]/TMath::Sqrt(massarray[i10]*massarray[i10]+momentum[9]*momentum[9]);
for (Int_t i11=0; i11<3;i11++) {
+ beta[10]=momentum[10]/TMath::Sqrt(massarray[i11]*massarray[i11]+momentum[10]*momentum[10]);
for (Int_t i12=0; i12<3;i12++) {
+ beta[11]=momentum[11]/TMath::Sqrt(massarray[i12]*massarray[i12]+momentum[11]*momentum[11]);
for (Int_t i13=0; i13<3;i13++) {
+ beta[12]=momentum[12]/TMath::Sqrt(massarray[i13]*massarray[i13]+momentum[12]*momentum[12]);
for (Int_t i14=0; i14<3;i14++) {
- for (Int_t i15=0; i15<3;i15++) {
- beta[0]=momentum[0]/sqrt(massarray[i1]*massarray[i1]+momentum[0]*momentum[0]);
- beta[1]=momentum[1]/sqrt(massarray[i2]*massarray[i2]+momentum[1]*momentum[1]);
- beta[2]=momentum[2]/sqrt(massarray[i3]*massarray[i3]+momentum[2]*momentum[2]);
- beta[3]=momentum[3]/sqrt(massarray[i4]*massarray[i4]+momentum[3]*momentum[3]);
- beta[4]=momentum[4]/sqrt(massarray[i5]*massarray[i5]+momentum[4]*momentum[4]);
- beta[5]=momentum[5]/sqrt(massarray[i6]*massarray[i6]+momentum[5]*momentum[5]);
- beta[6]=momentum[6]/sqrt(massarray[i7]*massarray[i7]+momentum[6]*momentum[6]);
- beta[7]=momentum[7]/sqrt(massarray[i8]*massarray[i8]+momentum[7]*momentum[7]);
- beta[8]=momentum[8]/sqrt(massarray[i9]*massarray[i9]+momentum[8]*momentum[8]);
- beta[9]=momentum[9]/sqrt(massarray[i10]*massarray[i10]+momentum[9]*momentum[9]);
- beta[10]=momentum[10]/sqrt(massarray[i11]*massarray[i11]+momentum[10]*momentum[10]);
- beta[11]=momentum[11]/sqrt(massarray[i12]*massarray[i12]+momentum[11]*momentum[11]);
- beta[12]=momentum[12]/sqrt(massarray[i13]*massarray[i13]+momentum[12]*momentum[12]);
- beta[13]=momentum[13]/sqrt(massarray[i14]*massarray[i14]+momentum[13]*momentum[13]);
- beta[14]=momentum[14]/sqrt(massarray[i15]*massarray[i15]+momentum[14]*momentum[14]);
-
-
+ beta[13]=momentum[13]/TMath::Sqrt(massarray[i14]*massarray[i14]+momentum[13]*momentum[13]);
+ for (Int_t i15=0; i15<3;i15++) {
+ beta[14]=momentum[14]/TMath::Sqrt(massarray[i15]*massarray[i15]+momentum[14]*momentum[14]);
+
Float_t meantzero=0.;
Float_t sumAllweights=0.;
for (Int_t itz=0; itz<15;itz++) {
sqMomError[itz]=((1.-beta[itz]*beta[itz])*0.025)*((1.-beta[itz]*beta[itz])*0.025)*(tracktoflen[itz]/(0.299792*beta[itz]))*(tracktoflen[itz]/(0.299792*beta[itz])); // this gives the square of the momentum error in nanoseconds
sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); // total error for the current track
sumAllweights+=1./sqTrackError[itz];
- // redefining beta, it is useful in order to calculate t zero
- beta[itz]*=0.299792;
- timezero[itz]=(tracktoflen[itz]/beta[itz])-timeofflight[itz];
- weightedtimezero[itz]=((tracktoflen[itz]/beta[itz])-timeofflight[itz])/sqTrackError[itz];// weighted time zero for current track
+ timezero[itz]=(tracktoflen[itz]/(beta[itz]*0.299792))-timeofflight[itz];
+ weightedtimezero[itz]=((tracktoflen[itz]/(beta[itz]*0.299792))-timeofflight[itz])/sqTrackError[itz];// weighted time zero for current track
meantzero+=weightedtimezero[itz];
}
+
meantzero=meantzero/sumAllweights; // it is given in [ns]
Float_t dummychisquare=0.;
}
}
- cout << "true" << truparticle[0] << truparticle[1] << truparticle[2] << truparticle[3] << truparticle[4] << truparticle[5] << truparticle[6] << truparticle[7] << truparticle[8] << truparticle[9] <<truparticle[10] << truparticle[11] << truparticle[12] << truparticle[13] << truparticle[14] <<endl;
- cout << "best" << assparticle[0] << assparticle[1] << assparticle[2] << assparticle[3] << assparticle[4] << assparticle[5] << assparticle[6] << assparticle[7] << assparticle[8] << assparticle[9] <<assparticle[10] << assparticle[11] << assparticle[12] << assparticle[13] << assparticle[14] << endl;
- if(truparticle[0]==assparticle[0] && truparticle[1]==assparticle[1] && truparticle[2]==assparticle[2] && truparticle[3]==assparticle[3] && truparticle[4]==assparticle[4]&& truparticle[5]==assparticle[5] && truparticle[6]==assparticle[6] && truparticle[7]==assparticle[7] && truparticle[8]==assparticle[8] && truparticle[9]==assparticle[9] &&truparticle[10]==assparticle[10] && truparticle[11]==assparticle[11] && truparticle[12]==assparticle[12] && truparticle[13]==assparticle[13] && truparticle[14]==assparticle[14]) ngood+=1;
+ cout << "true" << truparticle[0] << truparticle[1] << truparticle[2] << truparticle[3] << truparticle[4] << truparticle[5] << truparticle[6] << truparticle[7] << truparticle[8] << truparticle[9] <<truparticle[10] << truparticle[11] << truparticle[1
+2] << truparticle[13] << truparticle[14] <<endl;
+ cout << "best" << assparticle[0] << assparticle[1] << assparticle[2] << assparticle[3] << assparticle[4] << assparticle[5] << assparticle[6] << assparticle[7] << assparticle[8] << assparticle[9] <<assparticle[10] << assparticle[11] << assparticle[1
+2] << assparticle[13] << assparticle[14] << endl;
+ if(truparticle[0]==assparticle[0] && truparticle[1]==assparticle[1] && truparticle[2]==assparticle[2] && truparticle[3]==assparticle[3] && truparticle[4]==assparticle[4]&& truparticle[5]==assparticle[5] && truparticle[6]==assparticle[6] && trupart
+icle[7]==assparticle[7] && truparticle[8]==assparticle[8] && truparticle[9]==assparticle[9] &&truparticle[10]==assparticle[10] && truparticle[11]==assparticle[11] && truparticle[12]==assparticle[12] && truparticle[13]==assparticle[13] && truparticle[14]
+==assparticle[14]) ngood+=1;
if(truparticle[0]!=assparticle[0]) nmisidentified0+=1;
if(truparticle[1]!=assparticle[1]) nmisidentified1+=1;
if(truparticle[2]!=assparticle[2]) nmisidentified2+=1;
if(truparticle[13]!=assparticle[13]) nmisidentified13+=1;
if(truparticle[14]!=assparticle[14]) nmisidentified14+=1;
cout << "chisquare for current set" << chisquare << endl;
- //Float_t weight=1./chisquare;
- //htzerobest->Fill(t0best,weight);
htzerobest->Fill(t0best);
hchibest->Fill(chisquare);
Double_t dblechisquare=(Double_t)chisquare;
chisquare=999.;
t0best=999.;
- } // end for the current set. close if(istop==5)
+ } // end for the current set. close if(istop==15)
} // end condition on ipartold
}//end loop on events
- nmisidentified=(nmisidentified0+nmisidentified1+nmisidentified2+nmisidentified3+nmisidentified4+nmisidentified5+nmisidentified6+nmisidentified7+nmisidentified8+nmisidentified9+nmisidentified10+nmisidentified11+nmisidentified12+nmisidentified13+nmisidentified14);
+ nmisidentified=(nmisidentified0+nmisidentified1+nmisidentified2+nmisidentified3+nmisidentified4+nmisidentified5+nmisidentified6+nmisidentified7+nmisidentified8+nmisidentified9+nmisidentified10+nmisidentified11+nmisidentified12+nmisidentified13+nmisident
+ified14);
cout << "total misidentified " << nmisidentified << endl;
TFile *houtfile = new TFile(outFileName,"recreate");
houtfile->cd();
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff