ITS/AliITSFindTracksANN.C

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