ITS/oldmacros/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 "AliRunLoader.h"
   9   #include "AliTPCLoader.h"
  10   #include "AliITSLoader.h"
  11   #include "AliITS.h"
  12   #include "AliITSgeom.h"
  13   #include "AliITStrackerANN.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 AliRunLoader::Instance();
  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         AliESD* event = new AliESD;
 157         TTree* esdTree = (TTree*) tpcf->Get("esdTree");
 158         if (!esdTree) {
 159           cerr<<"no ESD tree found !\n";
 160           return 1;
 161         }
 162         esdTree->SetBranchAddress("ESD", &event);
 163
 164         // Loop on events
 165         TStopwatch timer;
 166         for (Int_t i = 0; i < nev; i++) {
 167                 cerr << "Processing event number: " << i << endl;
 168                 esdTree->GetEvent(i);
 169
 170                 rl->GetEvent(i);
 171
 172                 // Get clusters from file
 173                 TTree *cTree=itsl->TreeR();
 174                 if (!cTree) {
 175                         cerr<<"AliITSFindTracksANN.C : NO clusters tree!" << endl;
 176                         return 4;
 177                 }
 178
 179                 // Load clusters in tracker
 180                 tracker.LoadClusters(cTree);
 181
 182                 // Create array structure and arrange points in it
 183                 tracker.CreateArrayStructure(nsectors);
 184                 tracker.ArrangePoints("its.recpoints.txt");
 185                 tracker.StoreOverallMatches();
 186                 //tracker.PrintMatches(kFALSE);
 187
 188                 // ***************************
 189                 //  NEURAL TRACKING OPERATION
 190                 // ***************************
 191
 192                 Bool_t isStable = kFALSE;
 193                 Int_t i, nUnits = 0, nLinks = 0, nTracks = 0;
 194                 Int_t sectTracks = 0, totTracks = 0;
 195
 196                 // tracking through sectors
 197                 cout << endl;
 198                 Int_t sector, curv;
 199 //              nUnits = tracker.CreateNetwork(0, ncuts - 1);
 200 //              cout << endl << nUnits << " NEURONS CREATED" << endl;
 201 //              return;
 202
 203                 for (sector = 0; sector < nsectors; sector++) {
 204                         sectTracks = 0;
 205                         for(curv = 0; curv < ncuts; curv++) {
 206                                 // units creation
 207                                 nUnits = tracker.CreateNetwork(sector, curv);
 208                                 if (!nUnits) {
 209                                         cout << "no neurons --> skipped" << endl;
 210                                         continue;
 211                                 }
 212                                 // units connection
 213                                 nLinks = tracker.LinkNeurons();
 214                                 if (!nLinks) {
 215                                         cout << "no connections --> skipped" << endl;
 216                                         continue;
 217                                 }
 218                                 // neural updating
 219                                 for (i = 0;; i++) {
 220                                         isStable = tracker.Update();
 221                                         if (isStable) break;
 222                                 }
 223                                 // tracks saving
 224                                 tracker.CleanNetwork();
 225                                 nTracks = tracker.StoreTracks();
 226                                 cout << nUnits << " units, ";
 227                                 cout << nLinks << " connections, ";
 228                                 cout << i << " cycles --> ";
 229                                 cout << nTracks << " tracks" << endl;
 230                                 sectTracks += nTracks;
 231                         }
 232                         cout << "\n >>> Total tracks found in sector: " << sectTracks << endl;
 233                         totTracks += sectTracks;
 234                 }
 235                 tracker.ExportTracks("ITS.Neural.root");
 236                 cout << "\n\n--- Totally found " << totTracks << " tracks ---\n\n" << endl;
 237
 238                 // ***************************
 239
 240                 // End of operations: unload clusters
 241                 tracker.UnloadClusters();
 242
 243         }
 244         timer.Stop();
 245         timer.Print();
 246
 247         // Close files & delete objects
 248         delete event;
 249         tpcf->Close();
 250         itsf->Close();
 251         delete rl;
 252         return rc;
 253 }
 254
 255