ITS/AliITSFindTracksANN.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2   #include "Riostream.h"
   3   #include "TFile.h"
   4   #include "TTree.h"
   5   #include "TStopwatch.h"
   6
   7   #include "AliRun.h"
   8   #include "AliMagF.h"
   9   #include "AliRunLoader.h"
  10   #include "AliTPCLoader.h"
  11   #include "AliITSLoader.h"
  12   #include "AliITS.h"
  13   #include "AliITSgeom.h"
  14   #include "AliITStrackerANN.h"
  15   #include "AliESD.h"
  16 #endif
  17
  18 extern AliRun *gAlice;
  19
  20 Int_t AliITSFindTracksANN
  21 (Int_t nsectors = 10,      // number of azimutal sectors
  22  Int_t nev = 5)            // number of events to process
  23 {
  24
  25         // ==================================
  26         // ==== EVENT READING ===============
  27         // ==================================
  28
  29         // Remove an already existing Run Loader
  30         if (gAlice) {
  31                 delete gAlice->GetRunLoader();
  32                 delete gAlice;
  33                 gAlice=0;
  34         }
  35
  36         // Instance the new Run Loader
  37         AliRunLoader* rl = AliRunLoader::Open("galice.root");
  38         if (rl == 0x0) {
  39                 cerr<<"AliITSFindTracks.C : Can not open session RL=NULL"<< endl;
  40                 return 3;
  41         }
  42
  43         // Instance the ITS Loader
  44         AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
  45         if (itsl == 0x0) {
  46                 cerr<<"AliITSFindTracksV2.C : Can not get ITS loader"<<endl;
  47                 return 4;
  48         }
  49
  50         // Load the gAlice
  51         if (rl->LoadgAlice()) {
  52                 cerr<<"AliITSFindTracksV2.C : LoadgAlice returned error"<<endl;
  53                 delete rl;
  54                 return 3;
  55         }
  56
  57         // Set NECESSARY conversion constant for magnetic field
  58         AliKalmanTrack::SetConvConst
  59         (
  60                 1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()
  61         );
  62
  63         // Get the ITS data & geometry & recpoints
  64         AliITS *dITS = (AliITS*)rl->GetAliRun()->GetDetector("ITS");
  65         if (!dITS) {
  66                 cerr<<"AliITSFindClusters.C : Can not find the ITS detector !"<<endl;
  67                 return 6;
  68         }
  69         AliITSgeom *geom = dITS->GetITSgeom();
  70         itsl->LoadRecPoints("read");
  71
  72         // ==================================
  73         // ==== CURVATURE CUT DEFINITION ====
  74         // ==================================
  75
  76         // These values define the curvature cuts for all steps
  77         // within a sector.
  78         // For a greater clarity, the cuts are indicated in units
  79         // of transverse momentum (GeV/c) but these value have no
  80         // exact physical meaning, but are useful to understand
  81         // well what means a choice in the value of a certain
  82         // curvature constraint
  83         // NOTE: be careful to make sure that the 'ncuts' variable
  84         //       have the same value of the dimension of the allocated arrays
  85
  86         Int_t ncuts;
  87         Double_t *p, *cut;
  88
  89         ncuts = 6;
  90         p = new Double_t[6];
  91         cut = new Double_t[6];
  92         cut[ 0] = 0.0006;
  93         cut[ 1] = 0.0010;
  94         cut[ 2] = 0.0013;
  95         cut[ 3] = 0.0018;
  96         cut[ 4] = 0.0020;
  97         cut[ 5] = 0.0030;
  98
  99
 100         // ==========================
 101         // ==== OTHER PARAMETERS ====
 102         // ==========================
 103
 104
 105         Double_t helix[5]   = { 0.03, 0.02, 0.01, 0.03, 0.2 };
 106         Double_t theta2D[5] = { 1.5, 1.0, 1.0, 3.0, 10.0 };
 107         Double_t theta3D[5] = { 200., 200., 200., 200.0, 200.0 };
 108
 109
 110         /*
 111         Double_t helix[5]   = { 0.3, 0.3, 0.3, 0.5, 0.5 };
 112         Double_t theta2D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
 113         Double_t theta3D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
 114         */
 115
 116         Double_t temp   = 1.0;     // temperature parameter
 117         Double_t var    = 0.00001; // stabilization threshold
 118
 119         Double_t exp    = 20.0;    // straight-line excitator
 120         Double_t gtoc   = 6.0;     // gain/cost contribution ratio
 121
 122         Double_t min    = 0.4;     // minimum in random activation initialization
 123         Double_t max    = 0.6;     // maximum in random activation initialization
 124         Double_t actmin = 0.55;    // activation threshold for binary map conversion
 125
 126         // Instance the Tracker
 127         AliITStrackerANN tracker(geom, 2);
 128
 129         // Set cuts
 130         tracker.SetVertex(0.0, 0.0, 0.0);
 131
 132         tracker.SetCuts(ncuts, cut, theta2D, theta3D, helix);
 133         tracker.SetTemperature(temp);
 134         tracker.SetVariationLimit(var);
 135         tracker.SetGainToCostRatio(gtoc);
 136         tracker.SetWeightExponent(exp);
 137         tracker.SetInitInterval(min, max);
 138         tracker.SetActThreshold(actmin);
 139
 140         tracker.SetPolarInterval(45.0);
 141
 142         // Set number of events
 143         if (nev > rl->GetNumberOfEvents()) nev = rl->GetNumberOfEvents();
 144         Int_t rc = 0;
 145
 146         // Get ESD files
 147         TFile *itsf=TFile::Open("AliESDits.root","RECREATE");
 148         if ((!itsf)||(!itsf->IsOpen())) {
 149                 cerr << "Can't AliESDits.root !\n";
 150                 return 1;
 151         }
 152         TFile *tpcf=TFile::Open("AliESDtpc.root");
 153         if ((!tpcf)||(!tpcf->IsOpen())) {
 154                 cerr<<"Can't AliESDtpc.root !\n";
 155                 return 1;
 156         }
 157         AliESD* event = new AliESD;
 158         TTree* esdTree = (TTree*) tpcf->Get("esdTree");
 159         if (!esdTree) {
 160           cerr<<"no ESD tree found !\n";
 161           return 1;
 162         }
 163         esdTree->SetBranchAddress("ESD", &event);
 164
 165         // Loop on events
 166         TStopwatch timer;
 167         for (Int_t i = 0; i < nev; i++) {
 168                 cerr << "Processing event number: " << i << endl;
 169                 esdTree->GetEvent(i);
 170
 171                 rl->GetEvent(i);
 172
 173                 // Get clusters from file
 174                 TTree *cTree=itsl->TreeR();
 175                 if (!cTree) {
 176                         cerr<<"AliITSFindTracksANN.C : NO clusters tree!" << endl;
 177                         return 4;
 178                 }
 179
 180                 // Load clusters in tracker
 181                 tracker.LoadClusters(cTree);
 182
 183                 // Create array structure and arrange points in it
 184                 tracker.CreateArrayStructure(nsectors);
 185                 tracker.ArrangePoints("its.recpoints.txt");
 186                 tracker.StoreOverallMatches();
 187                 //tracker.PrintMatches(kFALSE);
 188
 189                 // ***************************
 190                 //  NEURAL TRACKING OPERATION
 191                 // ***************************
 192
 193                 Bool_t isStable = kFALSE;
 194                 Int_t i, nUnits = 0, nLinks = 0, nTracks = 0;
 195                 Int_t sectTracks = 0, totTracks = 0;
 196
 197                 // tracking through sectors
 198                 cout << endl;
 199                 Int_t sector, curv;
 200 //              nUnits = tracker.CreateNetwork(0, ncuts - 1);
 201 //              cout << endl << nUnits << " NEURONS CREATED" << endl;
 202 //              return;
 203
 204                 for (sector = 0; sector < nsectors; sector++) {
 205                         sectTracks = 0;
 206                         for(curv = 0; curv < ncuts; curv++) {
 207                                 // units creation
 208                                 nUnits = tracker.CreateNetwork(sector, curv);
 209                                 if (!nUnits) {
 210                                         cout << "no neurons --> skipped" << endl;
 211                                         continue;
 212                                 }
 213                                 // units connection
 214                                 nLinks = tracker.LinkNeurons();
 215                                 if (!nLinks) {
 216                                         cout << "no connections --> skipped" << endl;
 217                                         continue;
 218                                 }
 219                                 // neural updating
 220                                 for (i = 0;; i++) {
 221                                         isStable = tracker.Update();
 222                                         if (isStable) break;
 223                                 }
 224                                 // tracks saving
 225                                 tracker.CleanNetwork();
 226                                 nTracks = tracker.StoreTracks();
 227                                 cout << nUnits << " units, ";
 228                                 cout << nLinks << " connections, ";
 229                                 cout << i << " cycles --> ";
 230                                 cout << nTracks << " tracks" << endl;
 231                                 sectTracks += nTracks;
 232                         }
 233                         cout << "\n >>> Total tracks found in sector: " << sectTracks << endl;
 234                         totTracks += sectTracks;
 235                 }
 236                 tracker.ExportTracks("ITS.Neural.root");
 237                 cout << "\n\n--- Totally found " << totTracks << " tracks ---\n\n" << endl;
 238
 239                 // ***************************
 240
 241                 // End of operations: unload clusters
 242                 tracker.UnloadClusters();
 243
 244         }
 245         timer.Stop();
 246         timer.Print();
 247
 248         // Close files & delete objects
 249         delete event;
 250         tpcf->Close();
 251         itsf->Close();
 252         delete rl;
 253         return rc;
 254 }
 255
 256