ITS/ITSstoreFindableTracks.C

   1 void ITSstoreFindableTracks(const char *nfile = "galice", Int_t evnum = 0)
   2 {
   3         TFile *froot = new TFile(Form("%s.root", nfile), "READ");
   4
   5         gAlice = (AliRun*)froot->Get("gAlice");
   6         if (!gAlice) {
   7                 cout << "gAlice not found in file!!!" << endl;
   8                 return;
   9         }
  10
  11         AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
  12
  13         Int_t nparts = gAlice->GetEvent(evnum);
  14         cout << "Tracks saved in event " << evnum << ": " << nparts << endl;
  15
  16         TClonesArray *recPoints = ITS->RecPoints();
  17         TTree        *treeRec   = gAlice->TreeR();
  18    Int_t         ntracks   = gAlice->GetNtrack(); //FCA correction
  19         Int_t         nmodules  = treeRec->GetEntries();
  20         Int_t         modmin[6];
  21         Int_t         modmax[6];
  22
  23         Int_t layer, nlads, ndets;
  24         AliITSgeom *gm = ITS->GetITSgeom();
  25         for (layer = 0; layer < 6; layer++) {
  26                 nlads  = gm->GetNladders(layer+1);
  27                 ndets  = gm->GetNdetectors(layer+1);
  28                 modmin[layer] = gm->GetModuleIndex(layer+1, 1, 1);
  29                 modmax[layer] = gm->GetModuleIndex(layer+1, nlads, ndets);
  30         }
  31
  32         Int_t track, *goodITS[6];
  33         for (layer = 0; layer < 6; layer++) {
  34                 goodITS[layer] = new Int_t[ntracks];
  35                 for(track = 0; track < ntracks; track++) goodITS[layer][track] = 0;
  36         }
  37
  38         Int_t irec, index, npoints = 0;
  39         for (Int_t layer = 1; layer <= 6; layer++) {
  40                 for (Int_t mod = modmin[layer-1]; mod <= modmax[layer-1]; mod++) {
  41                         ITS->ResetRecPoints();
  42                         treeRec->GetEvent(mod);
  43                         cout << "\rLayer " << layer << ", module: " << mod << flush;
  44                         TObjArrayIter *iter = (TObjArrayIter*)recPoints->MakeIterator();
  45                         while ((recp = (AliITSRecPoint*)iter.Next())) {
  46                                 for (index = 0; index < 3; index++) {
  47                                         track = recp->GetLabel(index);
  48                                         if(track < 0) continue;
  49                                         if(track > ntracks) {
  50                                                 cerr << " Error on track number \n";
  51                                                 continue;
  52                                         }
  53                                         goodITS[layer-1][track] = 1;
  54                                 } //loop over tracks
  55                         } //loop over points
  56                 } //loop over modules
  57                 cout << endl;
  58         } //loop over layers
  59
  60         // histos filled with neural results
  61         TFile     *fnew = new TFile(Form("%s_fnd.root", nfile), "RECREATE");
  62         TParticle *part = 0;
  63         TTree     *tree = new TTree("tree", "Tree of tracks");
  64
  65         Int_t count = 0, prim = 0;
  66         Double_t pt = 0.0;
  67         tree->Branch("count", &count, "count/I");
  68         tree->Branch("prim", &prim, "prim/I");
  69         tree->Branch("pt", &pt, "pt/D");
  70
  71         for(track = 0; track < ntracks; track++) {
  72                 prim = 0;
  73                 count = 0;
  74                 part = gAlice->Particle(track);
  75                 if (!(track % 10000)) cout << "\rStoring track " << track << flush;
  76                 if (part->GetFirstMother() == -1) prim = 1;
  77                 for (layer = 0; layer < 6; layer++) count += goodITS[layer][track];
  78                 pt = part->Pt();
  79                 tree->Fill();
  80         }
  81
  82         fnew->cd();
  83         tree->Write("tree", TObject::kOverwrite);
  84         fnew->Close();
  85 }