TSamuraiFDC2Physics.cxx 14.3 KB
Newer Older
Adrien Matta's avatar
Adrien Matta committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
/*****************************************************************************
 * Copyright (C) 2009-2016   this file is part of the NPTool Project         *
 *                                                                           *
 * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
 * For the list of contributors see $NPTOOL/Licence/Contributors             *
 *****************************************************************************/

/*****************************************************************************
 * Original Author: Adrien MATTA  contact address: matta@lpccaen.in2p3.fr    *
 *                                                                           *
 * Creation Date  : October 2020                                             *
 * Last update    :                                                          *
 *---------------------------------------------------------------------------*
 * Decription:                                                               *
 *  This class hold SamuraiFDC2 treated data                                 *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * Comment:                                                                  *
 *                                                                           *
 *****************************************************************************/
#include "TSamuraiFDC2Physics.h"

//   STL
#include <sstream>
#include <iostream>
#include <cmath>
#include <stdlib.h>
#include <limits>

//   NPL
#include "RootInput.h"
#include "RootOutput.h"
#include "TAsciiFile.h"
#include "NPOptionManager.h"
#include "NPDetectorFactory.h"
Adrien Matta's avatar
Adrien Matta committed
36
#include "NPSystemOfUnits.h"
Adrien Matta's avatar
Adrien Matta committed
37
//   ROOT
Adrien Matta's avatar
Adrien Matta committed
38
using namespace NPUNITS;
Adrien Matta's avatar
Adrien Matta committed
39 40 41 42 43 44 45 46 47
///////////////////////////////////////////////////////////////////////////

ClassImp(TSamuraiFDC2Physics)
  ///////////////////////////////////////////////////////////////////////////
  TSamuraiFDC2Physics::TSamuraiFDC2Physics(){
    m_EventData         = new TSamuraiFDC2Data ;
    m_PreTreatedData    = new TSamuraiFDC2Data ;
    m_EventPhysics      = this ;
    //m_Spectra           = NULL;
48
    ToTThreshold = 180;
Adrien Matta's avatar
Adrien Matta committed
49 50 51 52 53 54 55 56 57 58
  }

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::BuildSimplePhysicalEvent(){
  BuildPhysicalEvent();
}

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::BuildPhysicalEvent(){
  PreTreat();
Adrien Matta's avatar
Adrien Matta committed
59 60
  RemoveNoise();

Adrien Matta's avatar
Adrien Matta committed
61
  // Map[detector & plane angle, vector of spatial information]
Adrien Matta's avatar
Adrien Matta committed
62 63 64 65
  static map<std::pair<unsigned int,double>, vector<double> > X ; 
  static map<std::pair<unsigned int,double>, vector<double> > Z ; 
  static map<std::pair<unsigned int,double>, vector<double> > R ; 
  X.clear();Z.clear();R.clear();
Adrien Matta's avatar
Adrien Matta committed
66 67
  unsigned int size = Detector.size();
  for(unsigned int i = 0 ; i < size ; i++){
68
    if(DriftLength[i]>0.1){
Adrien Matta's avatar
Adrien Matta committed
69 70 71 72
      int det = Detector[i];
      int layer = Layer[i];
      int wire = Wire[i]; 
      SamuraiDCIndex idx(det,layer,wire);
73
      std::pair<unsigned int, double> p(det,Wire_Angle[idx]);
Adrien Matta's avatar
Adrien Matta committed
74 75 76 77 78 79 80
      X[p].push_back(Wire_X[idx]); 
      Z[p].push_back(Wire_Z[idx]); 
      R[p].push_back(DriftLength[i]); 
    }
  }

  // Reconstruct the vector for each of the plane of each of the detector
81
  static double X0,X100,a,b; // store the BuildTrack2D results
Adrien Matta's avatar
Adrien Matta committed
82 83
  static map<std::pair<unsigned int,double>, TVector3 > VX0 ;  
  static map<std::pair<unsigned int,double>, TVector3 > VX100 ;  
84
  static map<std::pair<unsigned int,double>, int > MultPlane ;  
Adrien Matta's avatar
Adrien Matta committed
85
  VX0.clear();VX100.clear();
Adrien Matta's avatar
Adrien Matta committed
86
  for(auto it = X.begin();it!=X.end();++it){
87
    m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b); 
88 89 90 91
    // very small a means track perpendicular to the chamber, what happen when there is pile up

    if(abs(a)>1000)
      PileUp++;
Adrien Matta's avatar
Adrien Matta committed
92

93 94 95 96 97 98 99 100 101 102 103 104
    MultPlane[it->first] = X[it->first].size() ;
    Mult+=X[it->first].size();
    // Position at z=0
    TVector3 P(X0,0,0);
    P.RotateZ(it->first.second);
    VX0[it->first]=P;
    // Direction of the vector in the plane
    TVector3 D = TVector3(X100,0,0);
    D.RotateZ(it->first.second);
    VX100[it->first]=D;

  }
Adrien Matta's avatar
Adrien Matta committed
105
  // Reconstruct the central position (z=0) for each detector
Adrien Matta's avatar
Adrien Matta committed
106 107
  static map<unsigned int,vector<TVector3> > C ;  
  C.clear();
Adrien Matta's avatar
Adrien Matta committed
108
  TVector3 P;
109

Adrien Matta's avatar
Adrien Matta committed
110 111
  for(auto it1 = VX0.begin();it1!=VX0.end();++it1){
    for(auto it2 = it1;it2!=VX0.end();++it2){
Adrien Matta's avatar
Adrien Matta committed
112
      if(it1!=it2 && it1->first.first==it2->first.first){// different plane, same detector
113
        m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
114 115
        if(P.X()!=-10000)
          C[it1->first.first].push_back(P);
Adrien Matta's avatar
Adrien Matta committed
116
      }
117
    }
Adrien Matta's avatar
Adrien Matta committed
118
  }
119
  // Reconstruct the position at z=100 for each detector
Adrien Matta's avatar
Adrien Matta committed
120 121 122 123 124
  static map<unsigned int,vector<TVector3> > C100 ;  
  C100.clear();
  for(auto it1 = VX100.begin();it1!=VX100.end();++it1){
    for(auto it2 = it1;it2!=VX100.end();++it2){
      if(it1!=it2 && it1->first.first==it2->first.first){// different plane, same detector
125
        m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
126 127 128

        if(P.X()!=-10000)
          C100[it1->first.first].push_back(P);
Adrien Matta's avatar
Adrien Matta committed
129
      }
130
    }
Adrien Matta's avatar
Adrien Matta committed
131
  }
Adrien Matta's avatar
Adrien Matta committed
132 133
  // Build the Reference position by averaging all possible pair 
  size = C[2].size();
Adrien Matta's avatar
Adrien Matta committed
134
  double PosX100,PosY100;
Adrien Matta's avatar
Adrien Matta committed
135 136
  if(size){
    PosX=0;
Adrien Matta's avatar
Adrien Matta committed
137
    PosY=0;
Adrien Matta's avatar
Adrien Matta committed
138 139
    PosX100=0;
    PosY100=0;
Adrien Matta's avatar
Adrien Matta committed
140 141
    for(unsigned int i = 0 ; i < size ; i++){
      PosX+= C[2][i].X(); 
Adrien Matta's avatar
Adrien Matta committed
142
      PosY+= C[2][i].Y(); 
Adrien Matta's avatar
Adrien Matta committed
143 144
      PosX100+= C100[2][i].X(); 
      PosY100+= C100[2][i].Y(); 
145
    //       cout << C[2][i].X() << " (" << C[2][i].Y() << ") ";
Adrien Matta's avatar
Adrien Matta committed
146
    } 
147 148
   // cout << endl;
    MultMean=size;
Adrien Matta's avatar
Adrien Matta committed
149
    // Mean position at Z=0
150 151
    PosX=PosX/size; 
    PosY=PosY/size; 
Adrien Matta's avatar
Adrien Matta committed
152
    // Mean position at Z=100
153 154
    PosX100=PosX100/size; 
    PosY100=PosY100/size; 
Adrien Matta's avatar
Adrien Matta committed
155 156 157 158 159
    // Compute ThetaX, angle between the Direction vector projection in XZ with
    // the Z axis
    ThetaX=atan((PosX100-PosX)/100.);
    // Compute PhiY, angle between the Direction vector projection in YZ with
    // the Z axis
160
    PhiY=atan((PosY100-PosY)/100.);
Adrien Matta's avatar
Adrien Matta committed
161
    Dir=TVector3(PosX100-PosX,PosY100-PosY,100).Unit();
Adrien Matta's avatar
Adrien Matta committed
162
  }
Adrien Matta's avatar
Adrien Matta committed
163

Adrien Matta's avatar
Adrien Matta committed
164 165 166 167 168
  return;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::PreTreat(){
  ClearPreTreatedData();
Adrien Matta's avatar
Adrien Matta committed
169 170
  static CalibrationManager* Cal = CalibrationManager::getInstance();
  static string channel;
171
  // one map per detector
Adrien Matta's avatar
Adrien Matta committed
172 173 174 175
  map<unsigned int, vector<double> > X ; 
  map<unsigned int, vector<double> > Z ; 
  map<unsigned int, vector<double> > R ; 

Adrien Matta's avatar
Adrien Matta committed
176 177
  unsigned int size = m_EventData->Mult();
  for(unsigned int i = 0 ; i < size ; i++){
Adrien Matta's avatar
Adrien Matta committed
178 179 180
    // EDGE=1 is the leading edge, IE, the real time.
    // EDGE=0 is the trailing edge, so it helps build Tot
    if(m_EventData->GetEdge(i)==1){
Adrien Matta's avatar
Adrien Matta committed
181 182 183 184 185 186 187
      int det   = m_EventData->GetDetectorNbr(i); 
      int layer = m_EventData->GetLayerNbr(i); 
      int wire  = m_EventData->GetWireNbr(i); 
      double time = m_EventData->GetTime(i);
      double etime = 0;
      // look for matching trailing edge   
      for(unsigned int j = 0 ; j < size ; j++){
Adrien Matta's avatar
Adrien Matta committed
188
        if(m_EventData->GetEdge(j)==0){
Adrien Matta's avatar
Adrien Matta committed
189 190 191 192 193 194 195 196
          int edet   = m_EventData->GetDetectorNbr(j); 
          int elayer = m_EventData->GetLayerNbr(j); 
          int ewire  = m_EventData->GetWireNbr(j); 
          // same wire
          if(wire==ewire && layer==elayer && det==edet){
            etime = m_EventData->GetTime(j); 
          }    
        }
Adrien Matta's avatar
Adrien Matta committed
197 198 199 200
        if(etime && etime>time)
          break;
        else
          etime=0;
Adrien Matta's avatar
Adrien Matta committed
201 202
      }
      // a valid wire must have an edge
Adrien Matta's avatar
Adrien Matta committed
203
      if(etime && time && etime-time>ToTThreshold){
Adrien Matta's avatar
Adrien Matta committed
204 205 206 207
        Detector.push_back(det);
        Layer.push_back(layer);       
        Wire.push_back(wire);
        Time.push_back(time);
Adrien Matta's avatar
Adrien Matta committed
208
        ToT.push_back(etime-time);
Adrien Matta's avatar
Adrien Matta committed
209 210
        channel="SamuraiFDC2/L" + NPL::itoa(layer);
        DriftLength.push_back(Cal->ApplySigmoid(channel,etime));
Adrien Matta's avatar
Adrien Matta committed
211
      }
Adrien Matta's avatar
Adrien Matta committed
212 213 214 215 216
    }

  }
  return;
}
Adrien Matta's avatar
Adrien Matta committed
217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::RemoveNoise(){
  // Remove the noise by looking if a matching wire exist in the adjacent layer
  // this done by looking at the closest plane with the same orientation
  unsigned int size = Detector.size(); 
  for(unsigned int i = 0 ; i < size ; i++){
    bool match=false;
    int det = Detector[i];
    int layer = Layer[i];
    int wire = Wire[i];
    // look for matching adjacent wire   

    for(unsigned int j = 0 ; j < size ; j++){
      int adet = Detector[j];
      int alayer = Layer[j];
      int awire = Wire[j];
      bool blayer = false;
      if(layer%2==0){
        if(layer+1==alayer)
          blayer=true;
      }

      else{
        if(layer-1==alayer)
          blayer=true;
      }
Adrien Matta's avatar
Adrien Matta committed
243

Adrien Matta's avatar
Adrien Matta committed
244 245 246 247 248
      if(det==adet && blayer && abs(wire-awire)<=1){
        match=true;
        break;
      }
    }
Adrien Matta's avatar
Adrien Matta committed
249 250 251 252 253 254

    if(match)
      Matched.push_back(true);
    else
      Matched.push_back(false);
  }
Adrien Matta's avatar
Adrien Matta committed
255 256
  return;
}
Adrien Matta's avatar
Adrien Matta committed
257 258
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::Clear(){
259 260
  MultMean=0;
  PileUp=0;
Adrien Matta's avatar
Adrien Matta committed
261 262
  Mult=0;
  PosX=PosY=-10000;
Adrien Matta's avatar
Adrien Matta committed
263
  DriftLength.clear();
Adrien Matta's avatar
Adrien Matta committed
264 265 266 267 268 269 270
  Detector.clear();
  Layer.clear();
  Wire.clear();
  Time.clear();
  ToT.clear();
  ParticleDirection.clear();
  MiddlePosition.clear();
Adrien Matta's avatar
Adrien Matta committed
271
  Matched.clear();
Adrien Matta's avatar
Adrien Matta committed
272 273 274 275 276 277 278
}
///////////////////////////////////////////////////////////////////////////

////   Innherited from VDetector Class   ////

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::ReadConfiguration(NPL::InputParser parser){
Adrien Matta's avatar
Adrien Matta committed
279 280 281 282 283 284 285 286 287 288 289 290 291
  vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SAMURAIFDC2");
  if(NPOptionManager::getInstance()->GetVerboseLevel())
    cout << "//// " << blocks.size() << " detector(s) found " << endl; 

  vector<string> token= {"XML"};

  for(unsigned int i = 0 ; i < blocks.size() ; i++){
    cout << endl << "////  Samurai FDC2 (" << i+1 << ")" << endl;
    string xmlpath = blocks[i]->GetString("XML");
    NPL::XmlParser xml;
    xml.LoadFile(xmlpath);
    AddDC("SAMURAIFDC2",xml);
  }
Adrien Matta's avatar
Adrien Matta committed
292
}
Adrien Matta's avatar
Adrien Matta committed
293 294 295

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::AddDC(string name, NPL::XmlParser& xml){
296
  std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName(name);  
Adrien Matta's avatar
Adrien Matta committed
297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
  // FDC2 case
  if(name=="SAMURAIFDC2"){
    unsigned int det=2;
    unsigned int size = b.size();
    for(unsigned int i = 0 ; i < size ; i++){
      unsigned int layer = b[i]->AsInt("layer"); 
      unsigned int wire  = b[i]->AsInt("wireid"); 
      double X = b[i]->AsDouble("wirepos");  
      double Z = b[i]->AsDouble("wirez");  
      string sDir = b[i]->AsString("anodedir");
      double T=0;
      if(sDir=="X")
        T=0*deg;
      else if(sDir=="Y")
        T=90*deg;
      else if(sDir=="U")
        T=30*deg;
      else if(sDir=="V")
        T=-30*deg;
      else{
        cout << "ERROR: Unknown layer orientation for Samurai FDC2"<< endl;
        exit(1);
      }
Adrien Matta's avatar
Adrien Matta committed
320 321 322
      SamuraiDCIndex idx(det,layer,wire);
      Wire_X[idx]=X;
      Wire_Z[idx]=Z;
323
      Wire_Angle[idx]=T;
Adrien Matta's avatar
Adrien Matta committed
324 325 326 327
    }
  }
}

Adrien Matta's avatar
Adrien Matta committed
328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::InitSpectra(){  
  //m_Spectra = new TSamuraiFDC2Spectra(m_NumberOfDetector);
}

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::FillSpectra(){  
  //  m_Spectra -> FillRawSpectra(m_EventData);
  //  m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
  //  m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::CheckSpectra(){  
  //  m_Spectra->CheckSpectra();  
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::ClearSpectra(){  
  // To be done
}
///////////////////////////////////////////////////////////////////////////
map< string , TH1*> TSamuraiFDC2Physics::GetSpectra() {
  /*  if(m_Spectra)
      return m_Spectra->GetMapHisto();
      else{
      map< string , TH1*> empty;
      return empty;
      }*/
  map< string , TH1*> empty;
  return empty;

} 

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::WriteSpectra(){
  // m_Spectra->WriteSpectra();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::AddParameterToCalibrationManager(){
Adrien Matta's avatar
Adrien Matta committed
366 367 368 369 370 371 372
  CalibrationManager* Cal = CalibrationManager::getInstance();

  // each layer
  for( int l = 0 ; l < 14 ; ++l){
    Cal->AddParameter("SamuraiFDC2", "L"+ NPL::itoa(l),"FDC2_L"+ NPL::itoa(l));
  }

Adrien Matta's avatar
Adrien Matta committed
373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419
}

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::InitializeRootInputRaw(){
  TChain* inputChain = RootInput::getInstance()->GetChain()   ;
  inputChain->SetBranchStatus( "SamuraiFDC2" , true );
  // The following line is necessary only for system were the tree is splitted
  // (older root version). The found argument silenced the Branches not found
  // warning for non splitted tree.
  if(inputChain->FindBranch("fDC_*"))
    inputChain->SetBranchStatus( "fDC_*",true);
  inputChain->SetBranchAddress( "SamuraiFDC2" , &m_EventData );

}

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::InitializeRootInputPhysics(){
}

///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::InitializeRootOutput(){
  TTree* outputTree = RootOutput::getInstance()->GetTree();
  outputTree->Branch( "SamuraiFDC2" , "TSamuraiFDC2Physics" , &m_EventPhysics );
}

////////////////////////////////////////////////////////////////////////////////
//            Construct Method to be pass to the DetectorFactory              //
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TSamuraiFDC2Physics::Construct(){
  return (NPL::VDetector*) new TSamuraiFDC2Physics();
}

////////////////////////////////////////////////////////////////////////////////
//            Registering the construct method to the factory                 //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
  class proxy_samuraiFDC2{
    public:
      proxy_samuraiFDC2(){
        NPL::DetectorFactory::getInstance()->AddToken("SAMURAIFDC2","Samurai");
        NPL::DetectorFactory::getInstance()->AddDetector("SAMURAIFDC2",TSamuraiFDC2Physics::Construct);
      }
  };

  proxy_samuraiFDC2 p_samuraiFDC2;
}