QA check for TOF hits updated
authorvicinanz <vicinanz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Apr 2002 10:05:17 +0000 (10:05 +0000)
committervicinanz <vicinanz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Apr 2002 10:05:17 +0000 (10:05 +0000)
TOF/AliTOFanalyzeHits.C [new file with mode: 0644]
TOF/AliTOFconfig.C [new file with mode: 0644]
TOF/AliTOFtest.C [new file with mode: 0644]
TOF/tzero15tracks.C

diff --git a/TOF/AliTOFanalyzeHits.C b/TOF/AliTOFanalyzeHits.C
new file mode 100644 (file)
index 0000000..afc5b99
--- /dev/null
@@ -0,0 +1,219 @@
+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;
+
+}
+
diff --git a/TOF/AliTOFconfig.C b/TOF/AliTOFconfig.C
new file mode 100644 (file)
index 0000000..0ca0219
--- /dev/null
@@ -0,0 +1,370 @@
+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));
+}
diff --git a/TOF/AliTOFtest.C b/TOF/AliTOFtest.C
new file mode 100644 (file)
index 0000000..9ccaf48
--- /dev/null
@@ -0,0 +1,253 @@
+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;
+
+}
index a11668b..577838b 100644 (file)
@@ -3,11 +3,21 @@ void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10)
   
   //
   // 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");
@@ -136,7 +146,7 @@ void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10)
            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;
@@ -167,49 +177,47 @@ void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10)
              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.;
@@ -255,9 +263,13 @@ void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10)
                }
              } 
              
-             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;
@@ -274,8 +286,6 @@ void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10)
              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;
@@ -285,7 +295,7 @@ void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10)
              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
@@ -300,7 +310,8 @@ void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10)
     
   }//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();