ITS/AliITSStoreFindableTracksCompiled.C

   1 #include <fstream.h>
   2
   3 #include <TClassTable.h>
   4 #include <TFile.h>
   5 #include <TTree.h>
   6 #include <TString.h>
   7 #include <TClonesArray.h>
   8 #include <TParticle.h>
   9
  10 #include "AliRun.h"
  11 #include "AliITS.h"
  12 #include "AliITSgeom.h"
  13 #include "AliITSRecPoint.h"
  14
  15 Int_t AliITSStoreFindableTracksCompiled
  16 (Int_t nMinClusters = 5, const Text_t *evname = "galice", Int_t evnum = 0)
  17 {
  18         // Make sure that ALICE objects are loaded
  19         if (gAlice) {
  20                 delete gAlice;
  21                 gAlice = 0;
  22         }
  23
  24         // Define the names of all involved files
  25         TString strEventFile(evname);
  26         TString strOutputFile(evname);
  27         strEventFile.Append(".root");
  28         strOutputFile.Append("_tracks_");
  29         strOutputFile += evnum;
  30         strOutputFile.Append(".root");
  31
  32         // Connect the Root Galice file containing Geometry, Kine and Hits
  33         TFile *fileEvent = (TFile*)gROOT->GetListOfFiles()->FindObject(strEventFile);
  34         if (!fileEvent) fileEvent = new TFile(strEventFile,"UPDATE");
  35
  36         // Get AliRun object from file
  37         gAlice = (AliRun*)fileEvent->Get("gAlice");
  38         if (gAlice) cout << "OK, found an AliRun object in file" << endl;
  39
  40         // Get ITS related objects and data
  41         AliITS* ITS =(AliITS *)gAlice->GetDetector("ITS");
  42         if (!ITS) {
  43                 cerr << "ITS object not found!" << endl;
  44                 return 1;
  45         }
  46         AliITSgeom *geometry = ITS->GetITSgeom();
  47         if (!geometry) {
  48                 cerr << "ITS geometry object not found!" << endl;
  49                 return 2;
  50         }
  51
  52         // Count the number of modules per layer
  53         Int_t nLadders, nDetectors, mod_min[6], mod_max[6];
  54         for(Int_t i = 0; i < 6; i++) {
  55                 nLadders = geometry->GetNladders(i + 1);
  56                 nDetectors = geometry->GetNdetectors(i + 1);
  57                 mod_min[i] = geometry->GetModuleIndex(i + 1, 1, 1);
  58                 mod_max[i] = geometry->GetModuleIndex(i + 1, nLadders, nDetectors);
  59         }
  60
  61         // Load event and ITS recpoints
  62         Int_t nParticles = gAlice->GetEvent(evnum);
  63         cout << "Event number: " << evnum << endl;
  64         cout << "# particles : " << nParticles <<endl;
  65         if (nParticles <= 0) {
  66                 cerr << "Can't have <= 0 particles!" << endl;
  67                 return 3;
  68         }
  69         AliITSRecPoint *recp = 0;
  70         TClonesArray  *recPoints = ITS->RecPoints();
  71         TObjArray     *particles = gAlice->Particles();
  72         Int_t nTracks = gAlice->GetNtrack(); //FCA correction
  73         Bool_t *hitITSLayer[6];
  74         for (Int_t i = 0; i < 6; i++) {
  75                 hitITSLayer[i] = new Bool_t[nTracks];
  76                 for (Int_t j = 0; j < nTracks; j++) hitITSLayer[i][j] = kFALSE;
  77         }
  78
  79         // Load recpoints in event
  80         TTree *TR = gAlice->TreeR();
  81         if (!TR) {
  82                 cerr << "TreeR object not found!" << endl;
  83                 return 4;
  84         }
  85
  86         // Scan recpoints and define findable tracks
  87         Int_t nModules = (Int_t)TR->GetEntries(), nPoints = 0, nEmpty = 0;
  88         cout << "Found " << nModules;
  89         cout << " entries in the TreeR (must be one per module!)" << endl;
  90         for (Int_t layer = 1; layer <= 6; layer++) {
  91                 for (Int_t mod = mod_min[layer - 1]; mod <= mod_max[layer - 1]; mod++) {
  92                         ITS->ResetRecPoints();
  93                         TR->GetEntry(mod);
  94                         nPoints = recPoints->GetEntries();
  95                         if(!nPoints) {
  96                                 nEmpty++;
  97                                 continue;
  98                         }
  99                         for (Int_t point = 0; point < nPoints; point++) {
 100                                 recp = (AliITSRecPoint*)recPoints->UncheckedAt(point);
 101                                 for (Int_t it = 0; it < 3; it++) {
 102                                         Int_t track = recp->GetLabel(it);
 103                                         if(track < 0) continue;
 104                                         if(track > nTracks) {
 105                                                 cout << "Found track index " << track;
 106                                                 cout << " whilw gAlice->GetNtrack() = " << nTracks << endl;
 107                                                 continue;
 108                                         }
 109                                         hitITSLayer[layer - 1][track] = kTRUE;
 110                                 } // loop over recpoint labels
 111                         } //loop over points
 112                 } //loop over modules
 113         } //loop over layers
 114         cout << "Found " << nEmpty << " empty modules" << endl;
 115
 116         // Scan the file of tracks in TPC to retrieve the findable TPC tracks
 117         TString strLabelsTPC;
 118         Int_t label, pdg_code, nFindablesTPC = 0;
 119         Double_t dummy;
 120         ifstream tpc("good_tracks_tpc");
 121         while (tpc >> label >> pdg_code) {
 122                 for (Int_t i = 0; i < 6; i++) tpc >> dummy;
 123                 nFindablesTPC++;
 124                 strLabelsTPC.Append(Form("[%d]", label));
 125         }
 126
 127         // Define the TTree with tracks data by means of a set of variables
 128         Int_t nFindablesITS = 0, nFindablesITSTPC = 0;
 129         Int_t nhits, tpc_ok, mother, entry = 0;
 130         Double_t vx, vy, vz;
 131         Double_t px, py, pz, pt;
 132
 133         TTree *tree = new TTree("Tracks", "Findable tracks in ITS");
 134
 135         tree->Branch("vx", &vx, "vx/D");
 136         tree->Branch("vy", &vy, "vy/D");
 137         tree->Branch("vz", &vz, "vz/D");
 138         tree->Branch("px", &px, "px/D");
 139         tree->Branch("py", &py, "py/D");
 140         tree->Branch("pz", &pz, "pz/D");
 141         tree->Branch("pt", &pt, "pt/D");
 142         tree->Branch("label", &label, "label/I");
 143         tree->Branch("entry", &entry, "entry/I");
 144         tree->Branch("mother", &mother, "mother/I");
 145         tree->Branch("pdg_code", &pdg_code, "pdg_code/I");
 146         tree->Branch("nhits", &nhits, "nhits/I");
 147         tree->Branch("tpc_ok", &tpc_ok, "tpc_ok/I");
 148
 149         // Fill the tree
 150         cout << endl;
 151         TParticle *p = 0;
 152         for (Int_t i = 0; i < nTracks; i++) {
 153                 nhits = 0;
 154                 for (Int_t j = 0; j < 6; j++) if (hitITSLayer[j][i]) nhits++;
 155                 if (nhits < nMinClusters) continue;
 156                 p = gAlice->Particle(i);
 157                 px = p->Px();
 158                 py = p->Py();
 159                 pz = p->Pz();
 160                 pt = p->Pt();
 161                 vx = p->Vx();
 162                 vy = p->Vy();
 163                 vz = p->Vz();
 164                 mother = p->GetFirstMother();
 165                 cout << "Track " << i << " stored\r" << flush;
 166                 tpc_ok = (strLabelsTPC.Contains(Form("[%d]", i)));
 167                 pdg_code = p->GetPdgCode();
 168                 label = i;
 169                 nFindablesITS++;
 170                 if (tpc_ok) nFindablesITSTPC++;
 171                 tree->Fill();
 172                 entry++;
 173         }
 174
 175         // Save into a file
 176         TFile *fileOutput = new TFile(strOutputFile, "recreate");
 177         tree->Write();
 178         fileOutput->Close();
 179
 180         cout << "# findable tracks in TPC     : " << nFindablesTPC << endl;
 181         cout << "# findable tracks in ITS     : " << nFindablesITS << endl;
 182         cout << "# findable tracks in ITS+TPC : " << nFindablesITSTPC << endl;
 183
 184         return 0;
 185 }