ITS/AliITSNeuralFit.C

   1 void AliITSNeuralFit (Int_t field_factor = 1, Bool_t exact_pid = kTRUE,
   2                       const char *fname_chains = "its_chains.root",
   3                       const char *fname_save   = "its_neural.root",
   4                       const char *fname_points = "its_recpoints_v1.root")
   5 {
   6         Int_t i, k;
   7
   8         TStopwatch timer;
   9
  10         // Look for magnetic field
  11         // if 'factor' option is 0, the field is read from gAlice,
  12         // else it is converted into double.
  13         // if 'exact_pid' is true, the "galice.root" file is opened
  14         // to get informations about particles
  15
  16         Double_t factor = (Double_t)field_factor;
  17         Int_t nparticles = 0;
  18         TFile *fileAlice = 0;
  19         AliITSPid *pid = new AliITSPid(1);
  20
  21         // Remove an already existing Run Loader
  22         if (gAlice) {
  23                 delete gAlice->GetRunLoader();
  24                 delete gAlice;
  25                 gAlice=0;
  26         }
  27
  28         // Instance the new Run Loader
  29         AliRunLoader* rl = AliRunLoader::Open("galice.root");
  30         if (rl == 0x0) {
  31                 cerr<<"AliITSFindTracks.C : Can not open session RL=NULL"<< endl;
  32                 return 3;
  33         }
  34
  35         // Instance the ITS Loader
  36         AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
  37         if (itsl == 0x0) {
  38                 cerr<<"AliITSFindTracksV2.C : Can not get ITS loader"<<endl;
  39                 return 4;
  40         }
  41
  42         // Load the gAlice
  43         if (rl->LoadgAlice()) {
  44                 cerr<<"AliITSFindTracksV2.C : LoadgAlice returned error"<<endl;
  45                 delete rl;
  46                 return 3;
  47         }
  48
  49         rl->LoadKinematics();
  50         AliStack* stack = rl->Stack();
  51
  52         if (field_factor == 0) {
  53                 factor = rl->GetAliRun()->Field()->SolenoidField() / 2.;
  54         }
  55         else {
  56                 factor = (Double_t)field_factor;
  57         }
  58
  59         // Open files and trees and set branches
  60         TFile *fileP = new TFile(fname_points);
  61         TFile *fileC = new TFile(fname_chains);
  62         TTree *treeP = (TTree*)fileP->Get("TreeP");
  63         TTree *treeC = (TTree*)fileC->Get("TreeC");
  64         AliITSglobalRecPoint *pnt = 0, *newpt = 0;
  65         treeP->SetBranchAddress("Points", &pnt);
  66         Int_t l[6];
  67         treeC->SetBranchAddress("l0", &l[0]);
  68         treeC->SetBranchAddress("l1", &l[1]);
  69         treeC->SetBranchAddress("l2", &l[2]);
  70         treeC->SetBranchAddress("l3", &l[3]);
  71         treeC->SetBranchAddress("l4", &l[4]);
  72         treeC->SetBranchAddress("l5", &l[5]);
  73         Int_t nP = (Int_t)treeP->GetEntries();
  74         Int_t nC = (Int_t)treeC->GetEntries();
  75
  76
  77         // Store points into a TClonesArray
  78         TClonesArray arrayP("AliITSglobalRecPoint", nP);
  79         for (i = 0; i < nP; i++) {
  80                 treeP->GetEntry(i);
  81                 new(arrayP[i]) AliITSglobalRecPoint(*pnt);
  82         }
  83
  84
  85         // Generate new file of fitted tracks
  86         Double_t mass, dedx, mom;
  87         Int_t label, pdg;
  88         AliITSNeuralTrack *track = 0;
  89         TFile *fileT = new TFile(fname_save, "recreate");
  90         TTree *treeT = new TTree("TreeT", "Neural found tracks (fitted)");
  91         treeT->Branch("Tracks", "AliITSNeuralTrack", &track);
  92
  93         timer.Start();
  94
  95
  96         // Scan chains tree and fit tracks
  97         cout << "Selected field factor = " << factor << endl;
  98         cout << "Fitting..." << endl;
  99         for (i = 0; i < nC; i++) {
 100
 101                 // Get track
 102                 treeC->GetEntry(i);
 103                 track = new AliITSNeuralTrack;
 104                 for (k = 0; k < 6; k++) {
 105                         newpt = (AliITSglobalRecPoint*)arrayP[l[k]];
 106                         track->Insert(newpt);
 107                 }
 108
 109                 // Cook label and preliminary fit
 110                 track->AssignLabel();
 111                 track->SetFieldFactor(factor);
 112                 track->RiemannFit();
 113
 114                 // PID
 115                 pdg = TMath::Abs(((TParticle*)stack->Particle(lab))->GetPdgCode());
 116                 if (pdg == 0) pdg = 211;
 117                 switch (pdg) {
 118                         case 11: mass = 0.000511; break;
 119                         case 211: mass = 0.13957018; break;
 120                         case 321: mass = 0.493677; break;
 121                         case 2212: mass = 0.938272; break;
 122                         default:
 123                                 cout << "PDG code " << pdg << " not recognized. Skipping track..." << endl;
 124                                 continue;
 125                 }
 126                 track->SetPDGcode(pdg);
 127                 track->SetMass(mass);
 128
 129                 // Kalman Fit
 130                 track->SeedCovariance();
 131                 track->KalmanFit();
 132                 if (!track->PropagateTo(3.0)) continue;
 133                 treeT->Fill();
 134         }
 135         cerr << "DONE" << endl;
 136
 137         timer.Stop();
 138         timer.Print();
 139
 140         treeT->Write();
 141         fileT->Close();
 142 }
 143