AMSDST
myTrTrackPlus.cxx
Go to the documentation of this file.
1 // Authors: M.Duranti - INFN di Perugia
2 #include "myTrTrackPlus.h"
3 #include "TF1.h"
4 #include "TClass.h"
5 
6 using namespace std;
7 
8 static int my_int_pow(int base, int exp);
9 
10 int my_int_pow(int base, int exp){
11  if (exp<0) return -1;
12  if (base==0) return 0;
13  int ss=1;
14  if (base<0 && exp%2==1)ss=-1;
15  int bb=abs(base);
16  int out=1;
17  for(int ii=0;ii<exp;ii++)
18  out*=bb;
19  return out*ss;
20 }
21 
22 //------------- funzione usata per fittare (1/R_g - 1/R_mc)*R_mc vs R_mc --------------------------------------------
23 //TF1* fkC = new TF1("user", "[0]/([1]+exp([2]*x))", 0.1, 100);
24 //TF1* fkC = new TF1("user", "-0.00181176/(-0.994949+exp(-0.0666901*x))", 0.1, 100);
25 TF1* fkC = new TF1("user", "1.0/((0.00394707/(-1.0417+exp(0.112602*x)))/x + 1.0/x)", 0.1, 100);//molto vecchia!!!!
26 //-------------------------------------------------------------------------------------------------------------------
27 
28 //--------------------------------------------------------------------
29 
30 ClassImp(myTrTrackPlus);
31 #ifdef _WITHGBATCH_
32 ClassImp(myTrTrackPlusFiller);
33 #endif //#ifdef _WITHGBATCH_
34 
35 //--------------------------------------------------------------------
36 
37 myTrTrackPlus::myTrTrackPlus(){
38 #ifdef PDEBUG
39  printf("In myTrTrackPlus::myTrTrackPlus\n");
40 #endif
41  PRINTDEBUG;
42  init();
43  PRINTDEBUG;
44 }
45 
46 myTrTrackPlus::~myTrTrackPlus(){
47 #ifdef PDEBUG
48  printf("In myTrTrackPlus::~myTrTrackPlus\n");
49 #endif
50  PRINTDEBUG;
51  Clear();
52 }
53 
54 void myTrTrackPlus::Clear(Option_t* option){
55 #ifdef PDEBUG
56  printf("In myTrTrackPlus::Clear\n");
57 #endif
58  PRINTDEBUG;
59  _fitresultmap.clear();
60  PRINTDEBUG;
61  BitPatternJ = 0;
62  kDef = -999999;
63  k89 = -999999;
64  k8 = -999999;
65  k9 = -999999;
66  kinn = -999999;
67  kDefNoMS = -999999;
68  k89NoMS = -999999;
69  k8NoMS = -999999;
70  k9NoMS = -999999;
71  kinnNoMS = -999999;
72  kDefSameWeight = -999999;
73  k89SameWeight = -999999;
74  k8SameWeight = -999999;
75  k9SameWeight = -999999;
76  kinnSameWeight = -999999;
77  kUp = -999999;
78  kLow = -999999;
79  kExt = -999999;
80  kAl = -999999;
81  kCIE = -999999;
82  kPGCIE = -999999;
83  Refit=0;
84  PRINTDEBUG;
85  for (int ii=0; ii<9; ii++) {
86  GlobalCoordinateAL[ii].setp(-999999, -999999, -999999);
87  GlobalCoordinateNL[ii].setp(-999999, -999999, -999999);
88  }
89  for (int ii=0; ii<5; ii++) NoiseMinDist[ii].setp(-999999, -999999, -999999);
90 #if defined _vdev_
91  fill_n(QStatusL, 9, 0);
92 #endif
93  PRINTDEBUG;
94  LH=-999999;
95  PRINTDEBUG;
96  TruncMeanX = -999999;
97  TruncMeanY = -999999;
98  TruncMeanInnX = -999999;
99  TruncMeanInnY = -999999;
100  PRINTDEBUG;
101  Q_NPoints = -999999;
102  Q_RMS = -999999;
103  InnerQ_RMS = -999999;
104  factor = -999999;
105  btempcor_status= -999999;
106  PRINTDEBUG;
107  distHMx_inn=-999999;
108  distHMy_inn=-999999;
109  distHMx_L1=-999999;
110  distHMy_L1=-999999;
111  distHMx_L9=-999999;
112  distHMy_L9=-999999;
113  distHMx_L8=-999999;
114  distHMy_L8=-999999;
115  PRINTDEBUG;
116  distNHx_inn=-999999;
117  distNHy_inn=-999999;
118  distNHx_L1=-999999;
119  distNHy_L1=-999999;
120  distNHx_L9=-999999;
121  distNHy_L9=-999999;
122  distNHx_L8=-999999;
123  distNHy_L8=-999999;
124  PRINTDEBUG;
125  PointTrackCogEcal.setp(-999999, -999999, -999999);
126  DirTrackCogEcal.setd(-999999, -999999);
127  PointTrackEntryEcal.setp(-999999, -999999, -999999);
128  DirTrackEntryEcal.setd(-999999, -999999);
129  PointTrackExitEcal.setp(-999999, -999999, -999999);
130  DirTrackExitEcal.setd(-999999, -999999);
131  PRINTDEBUG;
132  PointTrackTopTrd.setp(-999999, -999999, -999999);
133  DirTrackTopTrd.setd(-999999, -999999);
134  PointTrackCenterTrd.setp(-999999, -999999, -999999);
135  DirTrackCenterTrd.setd(-999999, -999999);
136  PointTrackBottomTrd.setp(-999999, -999999, -999999);
137  DirTrackBottomTrd.setd(-999999, -999999);
138  PRINTDEBUG;
139  // B=-999999;
140  // SenThetaX=-999999;
141  // SenThetaY=-999999;
142  // L=-999999;
143  // pMDR=-999999;
144  // pMDR_known=-999999;
145  PRINTDEBUG;
146  return;
147 }
148 
149 void myTrTrackPlus::init(){
150 #ifdef PDEBUG
151  printf("In myTrTrackPlus::init\n");
152 #endif
153  PRINTDEBUG;
154  Clear();
155  PRINTDEBUG;
156  return;
157 }
158 
159 Double32_t myTrTrackPlus::GetErrRinv(int ittp){
160 
161  static std::map<int, fitresult>::iterator it;
162  if (ittp>=0) {
163  it = _fitresultmap.find(ittp);
164  if (it!=_fitresultmap.end()) return it->second.ErrRinv;
165  }
166 
167  return 0.0;
168 }
169 
170 Double32_t myTrTrackPlus::GetChisq(int ittp){
171 
172  static std::map<int, fitresult>::iterator it;
173  if (ittp>=0) {
174  it = _fitresultmap.find(ittp);
175  if (it!=_fitresultmap.end()) return it->second.Chisq;
176  }
177 
178  return 0.0;
179 }
180 
181 Double32_t myTrTrackPlus::GetNormChisqX(int ittp){
182 
183  static std::map<int, fitresult>::iterator it;
184  if (ittp>=0) {
185  it = _fitresultmap.find(ittp);
186  if (it!=_fitresultmap.end()) return it->second.NormChisqX;
187  }
188 
189  return 0.0;
190 }
191 
192 Double32_t myTrTrackPlus::GetNormChisqY(int ittp){
193 
194  static std::map<int, fitresult>::iterator it;
195  if (ittp>=0) {
196  it = _fitresultmap.find(ittp);
197  if (it!=_fitresultmap.end()) return it->second.NormChisqY;
198  }
199 
200  return 0.0;
201 }
202 
203 Double32_t myTrTrackPlus::GetRigidity(int ittp){
204 
205  static std::map<int, fitresult>::iterator it;
206  if (ittp>=0) {
207  it = _fitresultmap.find(ittp);
208  if (it!=_fitresultmap.end()) return it->second.Rigidity;
209  }
210 
211  return 0.0;
212 }
213 
214 Double32_t myTrTrackPlus::GetRigidityCorrectedByMagnField(int ittp){
215 
216  static std::map<int, fitresult>::iterator it;
217  if (ittp>=0) {
218  it = _fitresultmap.find(ittp);
219  if (it!=_fitresultmap.end()){
220  if(btempcor_status==0)
221  return (it->second.Rigidity)*factor;
222  else
223  return 0.0;
224  }
225  }
226 
227  return 0.0;
228 }
229 
230 int myTrTrackPlus::iTrTrackPar(int algo, int pattern, int refit){
231 
232  int index = ComputeIndex(algo, pattern, refit);
233 
234  if (_fitresultmap.find(index)!=_fitresultmap.end()) return index;
235  else return -99;
236 }
237 
238 int myTrTrackPlus::ComputeIndex(int algo, int pattern, int refit) {
239 
240  int _refit=0;
241  if ((refit%2)>=2) //0, 1 and 2 should be considered equivalent. 3 means that a coordinare rebuild has been called and so it a different rigidity (most probably this refit will be not included in ntuples...)
242  _refit=refit;
243 
244  int _pattern=0;
245  if (pattern>9){ //it is a base10 hit pattern
246  for (int kk=0; kk<9; kk++){
247  if (((pattern/my_int_pow(10,kk))%10)==1){ //OLD scheme
248  _pattern|=(1<<kk);
249  }
250  else if (((pattern/my_int_pow(10,kk))%10)==9){ // J-Scheme
251  _pattern|=(1<<(kk+9));
252  }
253  }
254  }
255  else {
256  _pattern=(pattern<<18);
257  }
258  // printf("pattern=%d -> _pattern=%x\n", pattern, _pattern);
259  // printf("index=%x\n", (algo+(_refit<<5)+(_pattern<<10)));
260 
261  // it's enought to shift
262  // refit by 5 bits (algo is up to 24 and 5 bits allows up to 31)
263  // pattern by 10 bits (algo is up to 24 and refit is up to 23 and 10 bits allow up to 31 for each)
264  return (algo+(_refit<<5)+(_pattern<<10));
265 }
266 
267 Short_t myTrTrackPlus::span() {
268 
269  Short_t span=-999;
270 
271  PRINTDEBUG;
272 
273  bool HasL1 = ((BitPatternJ)&(1<<0))>>0;
274  bool HasL2 = ((BitPatternJ)&(1<<1))>>1;
275  bool HasL9 = ((BitPatternJ)&(1<<8))>>8;
276  // printf("%d %d ", HasL1, HasL2);//only for debug
277  // for (int ii=2; ii<8; ii++) printf("%d ", ((BitPatternJ)&(1<<ii))>>ii);//only for debug
278  // printf("%d\n", HasL9);//only for debug
279 
280  if (HasL1&HasL9) {
281  // printf("HasL1&HasL9\n");//only for debug
282  if (kDef>0) {
283  span=0;
284  }
285  }
286  else if (HasL1) {
287  // printf("HasL1\n");//only for debug
288  if (kDef>0) {
289  span=1;
290  }
291  }
292  else if (HasL9) {
293  // printf("HasL9\n");//only for debug
294  if (kDef>0) {
295  if (HasL2) {
296  // printf("HasL2\n");//only for debug
297  span=2;
298  }
299  else span=3;
300  }
301  }
302  else {
303  // printf("Inner\n");//only for debug
304  if (kDef>0) {
305  if (HasL2) {
306  // printf("HasL2\n");//only for debug
307  span=4;
308  }
309  else span=5;
310  }
311  }
312 
313  PRINTDEBUG;
314 
315  return span;
316 }
317 
318 //---------------------------------------------------------------------------------------
319 
320 #ifdef _WITHGBATCH_
321 
322 myTrTrackPlusFiller::myTrTrackPlusFiller(){
323 #ifdef PDEBUG
324  printf("In myTrTrackPlusFiller::myTrTrackPlusFiller\n");
325 #endif
326  PRINTDEBUG;
327  init();
328  PRINTDEBUG;
329 }
330 
331 myTrTrackPlusFiller::~myTrTrackPlusFiller(){
332 #ifdef PDEBUG
333  printf("In myTrTrackPlusFiller::~myTrTrackPlusFiller\n");
334 #endif
335  PRINTDEBUG;
336  this->Clear();
337 }
338 
339 void myTrTrackPlusFiller::Clear(Option_t* option){
340 #ifdef PDEBUG
341  printf("In myTrTrackPlusFiller::Clear\n");
342 #endif
343  PRINTDEBUG;
344 #ifdef _WITHGBATCH_
345  track=0;
346 #endif
347  PRINTDEBUG;
348  return;
349 }
350 
351 void myTrTrackPlusFiller::init(){
352 #ifdef PDEBUG
353  printf("In myTrTrackPlusFiller::init\n");
354 #endif
355  PRINTDEBUG;
356  Clear();
357  PRINTDEBUG;
358  return;
359 }
360 
361 void myTrTrackPlusFiller::Fill(TrTrackR* _track, short int ipart, EcalShowerR* shower, Double_t beta, bool kMC, bool heavy){
362 #ifdef PDEBUG
363  printf("In myTrTrackPlusFiller::Fill\n");
364 #endif
365  PRINTDEBUG;
366 
367  track=_track;
368 
369  int n_valid_hit=1;
370  /*
371  if(track){
372  for (int i = 0; i < track->GetNhits(); i++) {
373  TrRecHitR *hit = track->GetHit(i);
374  if(!hit) continue;
375  n_valid_hit++;
376  }
377  }
378  */
379  //--------------------TRACKER-----------------------------------
380  if (track && n_valid_hit>0) {
381 
382  BitPatternJ=track->GetBitPatternJ();
383 
384  PRINTDEBUG;
385 
386  bool HasL1 = ((BitPatternJ)&(1<<0))>>0;
387  bool HasL2 = ((BitPatternJ)&(1<<1))>>1;
388  bool HasL9 = ((BitPatternJ)&(1<<8))>>8;
389  // printf("%d %d ", HasL1, HasL2);//only for debug
390  // for (int ii=2; ii<8; ii++) printf("%d ", ((BitPatternJ)&(1<<ii))>>ii);//only for debug
391  // printf("%d\n", HasL9);//only for debug
392 
393  PRINTDEBUG;
394 
395  for (int ii=0; ii<track->GetNhits(); ii++) {
396  TrRecHitR* rh = track->pTrRecHit(ii);
397  if(!rh) continue; // this is needed
398  int lay=rh->GetLayerJ()-1;
399  int mult=rh->GetResolvedMultiplicity();
400  if (mult<0) mult=0;
401  // mp.tkrechits.push_back(rh);
402  AMSPoint AtrrechitposA = rh->GetGlobalCoordinateA(mult);
403  GlobalCoordinateAL[lay].CopyAMSPoint(AtrrechitposA);
404  AMSPoint AtrrechitposN = rh->GetGlobalCoordinateN(mult);
405  GlobalCoordinateNL[lay].CopyAMSPoint(AtrrechitposN);
406 
407 #if !defined _B524_ && !defined _B530_ && !defined _B538_ && !defined _B550_ && !defined _B572_ && !defined _B580_ && !defined _B584_ && !defined _B594_ && !defined _B598_ && !defined _B610_ && !defined _B620_ && !defined _B630_
408  QStatusL[lay]=rh->GetQStatus();
409 #endif
410  }
411 
412  GetMinDistanceNoiseHits(track, NoiseMinDist);
413 
414  PRINTDEBUG;
415 
416  Fits(beta, kMC, heavy);
417 
418  //-----------------------------------------------------------------------
419 
420  PRINTDEBUG;
421 
422  AMSPoint APointTrackEntryEcal;
423  AMSDir ADirTrackEntryEcal;
424  track->Interpolate(-142.792, APointTrackEntryEcal, ADirTrackEntryEcal, kDef);
425  PointTrackEntryEcal.CopyAMSPoint(APointTrackEntryEcal);
426  DirTrackEntryEcal.CopyAMSDir(ADirTrackEntryEcal);
427 
428  AMSPoint APointTrackExitEcal;
429  AMSDir ADirTrackExitEcal;
430  track->Interpolate(-158.457, APointTrackExitEcal, ADirTrackExitEcal, kDef);
431  PointTrackExitEcal.CopyAMSPoint(APointTrackExitEcal);
432  DirTrackExitEcal.CopyAMSDir(ADirTrackExitEcal);
433 
434  AMSPoint APointTrackCogEcal;
435  AMSDir ADirTrackCogEcal;
436  if (shower) track->Interpolate(shower->CofG[2], APointTrackCogEcal, ADirTrackCogEcal, kDef);
437  PointTrackCogEcal.CopyAMSPoint(APointTrackCogEcal);
438  DirTrackCogEcal.CopyAMSDir(ADirTrackCogEcal);
439 
440  AMSPoint APointTrackTopTrd;
441  AMSDir ADirTrackTopTrd;
442  track->Interpolate(141.825, APointTrackTopTrd, ADirTrackTopTrd, kDef);
443  PointTrackTopTrd.CopyAMSPoint(APointTrackTopTrd);
444  DirTrackTopTrd.CopyAMSDir(ADirTrackTopTrd);
445 
446  AMSPoint APointTrackCenterTrd;
447  AMSDir ADirTrackCenterTrd;
448  track->Interpolate(0.5*(141.825+86.725), APointTrackCenterTrd, ADirTrackCenterTrd, kDef);
449  PointTrackCenterTrd.CopyAMSPoint(APointTrackCenterTrd);
450  DirTrackCenterTrd.CopyAMSDir(ADirTrackCenterTrd);
451 
452  AMSPoint APointTrackBottomTrd;
453  AMSDir ADirTrackBottomTrd;
454  track->Interpolate(86.725, APointTrackBottomTrd, ADirTrackBottomTrd, kDef);
455  PointTrackBottomTrd.CopyAMSPoint(APointTrackBottomTrd);
456  DirTrackBottomTrd.CopyAMSDir(ADirTrackBottomTrd);
457 
458  // if (shower) LH = TrkLH::gethead()->GetLikelihoodRatioElectronWG(ipart); //no more present
459  // printf("lktrk = %f\n", LH);
460  // TrkLHVar varv = TrkLH::gethead()->GetTrkLHVar(ipart);
461  // printf("Energy: %f, TrkLHVar: %f\n", me.EnergyE, varv.energy);
462 
463  PRINTDEBUG;
464 
465 #ifdef PATHONB
466  MagField* magfield = MagField::GetPtr();
467 
468  PathOnB(track, kDef, magfield, B, SenThetaX, SenThetaY, L, pMDR, pMDR_known, 0.1, false);
469  // printf("B=%f\n", "B);
470 #endif
471 
472  }
473 
474  PRINTDEBUG;
475 
476  return;
477 }
478 
479 void myTrTrackPlusFiller::Fits(Double_t beta, bool kMC, bool heavy){
480  Double_t __beta=999;
481  Refit=1;
482 
483  if (heavy){
484 
485  if (!kMC) { //if not MC
486 #ifndef _B524_ //whatever but B524
487 #ifdef _B550_ //tag B550
488  TrExtAlignDB::ResetExtAlign();
489  TrExtAlignDB::ForceFromTDV = 1;
490 #endif
491  // Refit=4;//PG DEPRECATED
492  // Refit=5;//CIEMAT DEPRECATED
493  // Refit=3; //PG, rebuild also coo
494  Refit=2; //PG
495 #if defined _B610_ || defined _B620_ || defined _vdev_
496  // Refit=1;
497  Refit=3;
498 #endif
499  //tolti per via della correzione del B-field
500  // if (Rigidity>50.0) Refit=1;//if they're "high Rigidity" there's no need to Refit with a different mass
501 #else //B524
502  Refit=2;
503  if (Rigidity>50.0) Refit=1;//if they're "high Rigidity" there's no need to Refit with a different mass
504 #endif
505 
506  }
507 
508  __beta=beta;
509  if (fabs(__beta)>1) __beta=1;
510  }
511 
512  PRINTDEBUG;
513 
514  bool proceed=true;
515 
516  if (track) {
517 
518  PRINTDEBUG;
519 
520  bool HasL1 = ((BitPatternJ)&(1<<0))>>0;
521  bool HasL2 = ((BitPatternJ)&(1<<1))>>1;
522  bool HasL9 = ((BitPatternJ)&(1<<8))>>8;
523  // printf("%d %d ", HasL1, HasL2);//only for debug
524  // for (int ii=2; ii<8; ii++) printf("%d ", ((BitPatternJ)&(1<<ii))>>ii);//only for debug
525  // printf("%d\n", HasL9);//only for debug
526 
527  PRINTDEBUG;
528 
529  if (HasL1&HasL9) {
530  k89 = AddResults( 1, 7, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
531  k89NoMS = AddResults(11, 7, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
532  k89SameWeight = AddResults(21, 7, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
533  }
534 
535  PRINTDEBUG;
536 
537 #ifndef DOHOWEVERALLFITS
538  proceed=false;
539 #endif
540 
541  PRINTDEBUG;
542 
543  if (HasL9 && proceed) {
544  k9 = AddResults( 1, 6, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
545  k9NoMS = AddResults(11, 6, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
546  k9SameWeight = AddResults(21, 6, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
547  }
548 
549  PRINTDEBUG;
550 
551 #ifndef DOHOWEVERALLFITS
552  proceed=false;
553 #endif
554 
555  PRINTDEBUG;
556 
557  if (HasL1 && proceed) {
558  k8 = AddResults( 1, 5, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
559  k8NoMS = AddResults(11, 5, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
560  k8SameWeight = AddResults(21, 5, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
561  }
562 
563  PRINTDEBUG;
564 
565 #ifndef DOHOWEVERALLFITS
566  proceed=false;
567 #endif
568 
569  PRINTDEBUG;
570 
571  if (proceed) {
572  kinn = AddResults( 1, 3, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
573  kinnNoMS = AddResults(11, 3, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
574  kinnSameWeight = AddResults(21, 3, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
575  }
576 
577  PRINTDEBUG;
578 
579  //------------------assign default one-----------------------------------
580 
581  if (HasL1&HasL9 && k89>=0) {
582  kDef = k89;
583  kDefNoMS = k89NoMS;
584  kDefSameWeight = k89SameWeight;
585  }
586  else if (HasL1 && k8>=0) {
587  kDef = k8;
588  kDefNoMS = k8NoMS;
589  kDefSameWeight = k8SameWeight;
590  }
591  else if (HasL9 && k9>=0) {
592  kDef = k9;
593  kDefNoMS = k9NoMS;
594  kDefSameWeight = k9SameWeight;
595  }
596  else {
597  kDef = kinn;
598  kDefNoMS = kinnNoMS;
599  kDefSameWeight = kinnSameWeight;
600  }
601 
602  if (kDef>=0) {
603  fitresult fr;
604  fr.Rigidity = getRigidity(track, kDef, false);
605  fr.ErrRinv = track->GetErrRinv(kDef);
606  fr.Chisq = getChisq(track, kDef);
607  fr.NormChisqX = getNormChisqX(track, kDef);
608  fr.NormChisqY = getNormChisqY(track, kDef);
609  AddResults(fr, 0, 0, 0);
610  }
611  if (kDefNoMS>=0) {
612  fitresult fr;
613  fr.Rigidity = getRigidity(track, kDefNoMS, false);
614  fr.ErrRinv = track->GetErrRinv(kDefNoMS);
615  fr.Chisq = getChisq(track, kDefNoMS);
616  fr.NormChisqX = getNormChisqX(track, kDefNoMS);
617  fr.NormChisqY = getNormChisqY(track, kDefNoMS);
618  AddResults(fr, 10, 0, 0);
619  }
620  if (kDefSameWeight>=0) {
621  fitresult fr;
622  fr.Rigidity = getRigidity(track, kDefSameWeight, false);
623  fr.ErrRinv = track->GetErrRinv(kDefSameWeight);
624  fr.Chisq = getChisq(track, kDefSameWeight);
625  fr.NormChisqX = getNormChisqX(track, kDefSameWeight);
626  fr.NormChisqY = getNormChisqY(track, kDefSameWeight);
627  AddResults(fr, 20, 0, 0);
628  }
629  //------------------------------------------------------------------------------------------------
630 
631  kUp = AddResults(1, 1, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
632  kLow = AddResults(1, 2, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
633 
634  if (HasL1&HasL9) {
635  kExt = AddResults(1, 4, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
636  }
637 
638  kAl = AddResults(2, 0, Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
639  kCIE = AddResults(0, 0, 10+Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
640  kPGCIE = AddResults(0, 0, 20+Refit, TrTrackR::DefaultMass, TrTrackR::DefaultCharge, __beta);
641 
642  PRINTDEBUG;
644  float _factor=0.;
645  btempcor_status= (Int_t) MagnetVarp::btempcor( _factor);
646  factor = (Double32_t) _factor;
647  PRINTDEBUG;
648 
649  }
650 
651  return;
652 }
653 
654 void myTrTrackPlusFiller::FillCharge(TrTrackR* _track, Double_t beta){
655 #ifdef PDEBUG
656  printf("In myEventFiller::Once_trackercharge\n");
657 #endif
658  PRINTDEBUG;
659 
660  track=_track;
661 
662  //-------------------------------------------------------------------
663  if (track) {
664 
665 #if defined _B524_ || defined _B550_
666  const static int flag=TrCharge::kTruncMean|TrCharge::kAll;
667  Double_t sqrtmeanx = sqrt(TrCharge::GetMean(flag, track, TrCharge::kX, beta).Mean);
668  TruncMeanX = -0.129643+0.202016*sqrtmeanx-0.00208005*sqrtmeanx*sqrtmeanx+2.60621e-05*sqrtmeanx*sqrtmeanx*sqrtmeanx;
669  Double_t sqrtmeany = sqrt(TrCharge::GetMean(flag, track, TrCharge::kY, beta).Mean);
670  TruncMeanY = -0.129643+0.202016*sqrtmeany-0.00208005*sqrtmeany*sqrtmeany+2.60621e-05*sqrtmeany*sqrtmeany*sqrtmeany;//the function is for X so in Y could be inaccurate...
671  const static int flaginn=TrCharge::kTruncMean|TrCharge::kInner;
672  sqrtmeanx = sqrt(TrCharge::GetMean(flaginn, track, TrCharge::kX, beta).Mean);
673  TruncMeanInnX = -0.129643+0.202016*sqrtmeanx-0.00208005*sqrtmeanx*sqrtmeanx+2.60621e-05*sqrtmeanx*sqrtmeanx*sqrtmeanx;
674  sqrtmeany = sqrt(TrCharge::GetMean(flaginn, track, TrCharge::kY, beta).Mean);
675  TruncMeanInnY = -0.129643+0.202016*sqrtmeany-0.00208005*sqrtmeany*sqrtmeany+2.60621e-05*sqrtmeany*sqrtmeany*sqrtmeany;//the function is for X so in Y could be inaccurate
676 #else
677  const static int flag=TrCharge::kTruncMean|TrCharge::kAll|TrCharge::kSqrt;
678  TruncMeanX = TrCharge::GetMean(flag, track, TrCharge::kX, beta, -1, TrClusterR::DefaultChargeCorrOpt).Mean;
679  TruncMeanY = TrCharge::GetMean(flag, track, TrCharge::kY, beta, -1, TrClusterR::DefaultChargeCorrOpt).Mean;
680  const static int flaginn=TrCharge::kTruncMean|TrCharge::kInner|TrCharge::kSqrt;
681  TruncMeanInnX = TrCharge::GetMean(flaginn, track, TrCharge::kX, beta, -1, TrClusterR::DefaultChargeCorrOpt).Mean;
682  TruncMeanInnY = TrCharge::GetMean(flaginn, track, TrCharge::kY, beta, -1, TrClusterR::DefaultChargeCorrOpt).Mean;
683 #endif
684  // printf("TruncMeanX,Y = %f,%f\n", TruncMeanX, TruncMeanY);//only for debug
685 
686  Q_NPoints = track->GetQ_NPoints(beta);
687  Q_RMS = track->GetQ_RMS(beta);
688  InnerQ_RMS = track->GetInnerQ_RMS(beta);
689 
690  PRINTDEBUG;
691 
692  }
693 
694  return;
695 }
696 
697 void myTrTrackPlusFiller::FillMoreHeavy(TrTrackR* _track){
698 #ifdef PDEBUG
699  printf("In myTrTrack::FillMoreHeavy\n");
700 #endif
701  PRINTDEBUG;
702 
703  track=_track;
704 
705  AMSEventR* pev = AMSEventR::Head();
706 
707  //---------------------TRACKER----------------------------------
708  if (track) {
709  // //--------------- ONLY FOR DEBUG -----------------------
710  // int nclux=0;
711  // int ncluy=0;
712  // int totclux=0;
713  // int totcluy=0;
714  // int used=0;
715  // int ptrzero=0;
716 
717  // if (pev->nTrTrack()==1) {
718  // // printf("------------------------------------------------\n");
719  // // printf("Hits = %d\n", track->GetNhits());
720  // for (int layer=0; layer<9; layer++) {
721  // TrRecHitR* rechit;
722  // rechit = track->GetHitL(layer);
723  // if (rechit) {
724  // // TrClusterR* xclu = rechit->GetXCluster();
725  // TrClusterR* xclu = pev->pTrCluster(rechit->iTrCluster('x'));
726  // if (xclu) {
727  // // printf("X) Layer = %d\n", layer);
728  // nclux++;
729  // }
730  // // TrClusterR* yclu = rechit->GetYCluster();
731  // TrClusterR* yclu = pev->pTrCluster(rechit->iTrCluster('y'));
732  // if (yclu) {
733  // // printf("Y) Layer = %d\n", layer);
734  // ncluy++;
735  // }
736  // }
737  // }
738  // // printf("Hit clusters: X = %d, Y = %d -> %d\n", nclux, ncluy, nclux+ncluy);
739  // // printf("Total clusters = %d -> residual clusters (subtracted the ones used in track) = %d\n", pev->nTrCluster(), pev->nTrCluster()-nclux-ncluy);
740  // for (int ii=0; ii<pev->nTrCluster(); ii++) {
741  // TrClusterR* clu = pev->pTrCluster(ii);
742  // if (clu) {
743  // if (!(clu->Used())) {
744  // int layer;
745  // layer = clu->GetLayer();
746  // if (clu->GetSide()==0) {
747  // // printf("X) Layer = %d\n", layer);
748  // totclux++;
749  // }
750  // else {
751  // // printf("Y) Layer = %d\n", layer);
752  // totcluy++;
753  // }
754  // }
755  // else {
756  // used++;
757  // }
758  // }
759  // else {
760  // ptrzero++;
761  // }
762  // }
763  // // printf("Residual clusters: X = %d, Y = %d -> %d\n", totclux, totcluy, totclux+totcluy);
764  // static int diff=9999;
765  // diff=fabs(pev->nTrCluster()-nclux-ncluy-totclux-totcluy);
766  // if (diff!=0 && diff!=1) {
767  // printf("Not OK!!!\n");
768  // printf("Hit clusters: X = %d, Y = %d -> %d\n", nclux, ncluy, nclux+ncluy);
769  // printf("Total clusters = %d -> residual clusters (subtracted the ones used in track) = %d\n", pev->nTrCluster(), pev->nTrCluster()-nclux-ncluy);
770  // printf("Residual clusters: X = %d, Y = %d -> %d\n", totclux, totcluy, totclux+totcluy);
771  // printf("Used clusters = %d, Null pointers = %d -> %d -->> %d\n", used, ptrzero, used+ptrzero, totclux+totcluy+used+ptrzero);
772  // }
773  // }
774  // //-----------------------------------------------------------
775 
776  PRINTDEBUG;
777 
778  distHMx_inn=9999999;
779  distHMy_inn=9999999;
780  distHMx_L1=9999999;
781  distHMy_L1=9999999;
782  distHMx_L9=9999999;
783  distHMy_L9=9999999;
784  distHMx_L8=9999999;
785  distHMy_L8=9999999;
786 
787  distNHx_inn=9999999;
788  distNHy_inn=9999999;
789  distNHx_L1=9999999;
790  distNHy_L1=9999999;
791  distNHx_L9=9999999;
792  distNHy_L9=9999999;
793  distNHx_L8=9999999;
794  distNHy_L8=9999999;
795 
796  static Double_t minore;
797 
798  minore=999999999;
799  for (int ii=2; ii<8; ii++) if (fabs(HitCloseTrack(ii, 0, "HM", pev, track, kDef))<minore) minore=HitCloseTrack(ii, 0, "HM", pev, track, kDef);
800  distHMx_inn=minore;
801  minore=999999999;
802  for (int ii=2; ii<8; ii++) if (fabs(HitCloseTrack(ii, 1, "HM", pev, track, kDef))<minore) minore=HitCloseTrack(ii, 1, "HM", pev, track, kDef);
803  distHMy_inn=minore;
804  distHMx_L1=HitCloseTrack(1, 0, "HM", pev, track, kDef);
805  distHMy_L1=HitCloseTrack(1, 1, "HM", pev, track, kDef);
806  distHMx_L9=HitCloseTrack(9, 0, "HM", pev, track, kDef);
807  distHMy_L9=HitCloseTrack(9, 1, "HM", pev, track, kDef);
808  distHMx_L8=HitCloseTrack(8, 0, "HM", pev, track, kDef);
809  distHMy_L8=HitCloseTrack(8, 1, "HM", pev, track, kDef);
810 
811  minore=999999999;
812  for (int ii=2; ii<8; ii++) if (fabs(HitCloseTrack(ii, 0, "NH", pev, track, kDef))<minore) minore=HitCloseTrack(ii, 0, "NH", pev, track, kDef);
813  distNHx_inn=minore;
814  minore=999999999;
815  for (int ii=2; ii<8; ii++) if (fabs(HitCloseTrack(ii, 1, "NH", pev, track, kDef))<minore) minore=HitCloseTrack(ii, 1, "NH", pev, track, kDef);
816  distNHy_inn=minore;
817  distNHx_L1=HitCloseTrack(1, 0, "NH", pev, track, kDef);
818  distNHy_L1=HitCloseTrack(1, 1, "NH", pev, track, kDef);
819  distNHx_L9=HitCloseTrack(9, 0, "NH", pev, track, kDef);
820  distNHy_L9=HitCloseTrack(9, 1, "NH", pev, track, kDef);
821  distNHx_L8=HitCloseTrack(8, 0, "NH", pev, track, kDef);
822  distNHy_L8=HitCloseTrack(8, 1, "NH", pev, track, kDef);
823 
824  }
825 
826  PRINTDEBUG;
827 
828  return;
829 }
830 
831 int myTrTrackPlusFiller::AddResults(int algo, int pattern, int refit, float mass, float chrg, float beta, float fixrig){
832 
833  int kRet=0;
834 #ifdef PDEBUG
835  printf("Event %d Run %d \n",AMSEventR::Head()->Event(),AMSEventR::Head()->Run());
836 #endif
837  kRet = track->iTrTrackPar(algo, pattern, refit, mass, chrg, beta, fixrig);
838  if (kRet<0) {
839 #ifdef PDEBUG
840  printf("kRet(%d,%d,%d,%f,%f,%f,%f) = %d\n", algo, pattern, refit, mass, chrg, beta, fixrig, kRet);
841 #endif
842  }
843  else {
844 #ifdef PDEBUG
845  printf("kRet(%d,%d,%d,%f,%f,%f,%f) = %d\n", algo, pattern, refit, mass, chrg, beta, fixrig, kRet);
846 #endif
847  fitresult fr;
848  fr.Rigidity = getRigidity(track, kRet, false);
849  fr.ErrRinv = track->GetErrRinv(kRet);
850  fr.Chisq = getChisq(track, kRet);
851  fr.NormChisqX = getNormChisqX(track, kRet);
852  fr.NormChisqY = getNormChisqY(track, kRet);
853  AddResults(fr, algo, pattern, refit);
854  }
855 
856  return kRet;
857 }
858 
859 void myTrTrackPlusFiller::AddResults(fitresult& fr, int algo, int pattern, int refit){
860 
861  _fitresultmap.insert(pair<int, fitresult>(ComputeIndex(algo, pattern, refit), fr));
862 
863  return;
864 }
865 
866 Double_t getRigidity(TrTrackR* _track, int par, bool rigcorr){
867 
868  static Double_t _rigidity;
869  _rigidity=99999;
870 
871  if (!_track) return _rigidity;
872 
873  _rigidity=_track->GetRigidity(par);
874 
875  if (!rigcorr) return _rigidity;
876 
877  static Double_t sign;
878  sign=_rigidity/fabs(_rigidity);
879  // printf("R: %f %f ", _rigidity, sign);
880 
881  if (fabs(_rigidity)<0.1) _rigidity=0.0;
882  else if (!(fabs(_rigidity)>100.0)) _rigidity=sign*fkC->GetX(fabs(_rigidity));
883 
884  return _rigidity;
885 }
886 
887 Double_t getChisq(TrTrackR* _track, int par){
888 
889  static Double_t chisq;
890  chisq=0.0;
891 
892  if (!_track) return chisq;
893 
894  if (par>=0)
895  chisq=_track->GetChisq(par);
896 
897  return chisq;
898 }
899 
900 Double_t getNormChisqX(TrTrackR* _track, int par){
901 
902  static Double_t chisq;
903  chisq=0.0;
904 
905  if (!_track) return chisq;
906 
907  if (par>=0) {
908  // chisq=_track->GetChisqX(par);
909  // chisq/=_track->GetNdofX(par);
910  chisq=_track->GetNormChisqX(par);
911  }
912 
913  return chisq;
914 }
915 
916 Double_t getNormChisqY(TrTrackR* _track, int par){
917 
918  static Double_t chisq;
919  chisq=0.0;
920 
921  if (!_track) return chisq;
922 
923  if (par>=0) {
924  // chisq=_track->GetChisqY(par);
925  // chisq/=_track->GetNdofY(par);
926  chisq=_track->GetNormChisqY(par);
927  }
928 
929  return chisq;
930 }
931 
932 //Double_t getLogExtChisq(TrTrackR *trk, Int_t alg = 1) {
933 Double_t getLogExtChisq(TrTrackR *trk, int k89, int k8, int k9) {
934  // Double_t errx = TRFITFFKEY.ErrX; // 25e-4;
935  Double_t erry = TRFITFFKEY.ErrY; // 13e-4;
936  Double_t fitw8 = TRFITFFKEY.FitwMsc[7]; // 324;
937  Double_t fitw9 = TRFITFFKEY.FitwMsc[8]; // 243;
938 
939  // Tuned with muon MC ; play with them
940  Double_t sig8 = 3;
941  Double_t sig9 = 4;
942  Double_t ecor = 1.7;
943  // Tuned with muon MC ; play with them
944 
945  // Int_t mf8 = trk->iTrTrackPar(alg, 5, 0); // With L1N
946  // Int_t mf9 = trk->iTrTrackPar(alg, 6, 0); // With L9
947  // Int_t mff = trk->iTrTrackPar(alg, 7, 0); // Full span
948 
949  // Double_t rgt = trk->GetRigidity(mff);
950  Double_t rgt = trk->GetRigidity(k89);
951  Double_t fms8 = fitw8/rgt;
952  Double_t fms9 = fitw9/rgt;
953 
954  Double_t ery8 = sqrt(sig8*sig8+fms8*fms8)*erry*ecor;
955  Double_t ery9 = sqrt(sig9*sig9+fms9*fms9)*erry*ecor;
956 
957  // Double_t res8 = trk->GetResidual(7, mf9).y();
958  // Double_t res9 = trk->GetResidual(8, mf8).y();
959 #ifdef _B524_
960  Double_t res8 = trk->GetResidual(7, k9).y();
961  Double_t res9 = trk->GetResidual(8, k8).y();
962 #else
963  Double_t res8 = trk->GetResidualO(1+7, k9).y();
964  Double_t res9 = trk->GetResidualO(1+8, k8).y();
965 #endif
966 
967  return res8*res8/ery8/ery8+res9*res9/ery9/ery9;
968 }
969 
970 void GetMinDistanceNoiseHits(TrTrackR *trk, myPoint noisemindist[]){
971  // search nearest noise hit for each layer
972  int trk_layer_pattern[9];// = {0,0,0,0,0,0,0,0,0};
973  fill_n(trk_layer_pattern, 9, 0);
974  Double_t noise_y[9];// = {1.e9,1.e9,1.e9,1.e9,1.e9,1.e9,1.e9,1.e9,1.e9};
975 
976  fill_n(noise_y, 9, 1e9);
977  int laypa = 1;
978  AMSPoint pnoise[9];
979  int btp= trk->GetBitPatternJ();;
980 
981  // cout<<" @@@@@@@@@@@@@@@@ btp "<<btp<<endl;
982 
983 
984  for (int layer=1; layer<=9; layer++) {
985  int bitexam = int(pow(2.,layer-1));
986  // cout<<" bitexam "<<bitexam<<endl;
987  bitexam = (bitexam & btp);
988  if (bitexam<1){
989  laypa=laypa*10;
990  continue;
991  }
992  // cout<<" bitexam "<<bitexam<<" laypa "<<laypa<<endl;
993  trk_layer_pattern[layer-1] = 1; // flags hitted layers
994  pnoise[layer-1] = TrTrackSelection::GetMinDist(trk, laypa, 1);
995  noise_y[layer-1]=fabs(pnoise[layer-1].y());
996  laypa=laypa*10;
997  }
998  Double_t minoise[4];
999  int nlayer = 0;
1000  minoise[0] = 7.;
1001  int lminoise[4];
1002  lminoise[0] = -1;
1003  for (int layer=1; layer<=9; layer++) {
1004  if(trk_layer_pattern[layer-1]==1) nlayer=nlayer+1; // found number hitted layer
1005 
1006  if (noise_y[layer-1]<minoise[0] && noise_y[layer-1] >0.001 && trk_layer_pattern[layer-1] ==1) {
1007  noisemindist[0].CopyAMSPoint(pnoise[layer-1]);
1008  minoise[0] = noise_y[layer-1];
1009  lminoise[0]=layer;
1010  }
1011  }
1012  minoise[1] = 7.;
1013  lminoise[1] = -1;
1014  for (int layer=1; layer<=9; layer++) {
1015  if (noise_y[layer-1]<minoise[1] && noise_y[layer-1] >0.001 && trk_layer_pattern[layer-1] ==1 && layer != lminoise[0]) {
1016  noisemindist[1].CopyAMSPoint(pnoise[layer-1]);
1017  minoise[1] = noise_y[layer-1]; lminoise[1]=layer;
1018  }
1019  }
1020  minoise[2]=7.;
1021  lminoise[2] = -1;
1022  for (int layer=1; layer<=9; layer++) {
1023  if (noise_y[layer-1]<minoise[2] && noise_y[layer-1] >0.001 && trk_layer_pattern[layer-1] ==1 && layer != lminoise[0] && layer != lminoise[1]) {
1024  noisemindist[2].CopyAMSPoint(pnoise[layer-1]);
1025  minoise[2] = noise_y[layer-1]; lminoise[2]=layer;
1026  }
1027  }
1028  minoise[3]=7.;
1029  lminoise[3] = -1;
1030  for (int layer=1; layer<=9; layer++) {
1031  if (noise_y[layer-1]<minoise[3] && noise_y[layer-1] >0.001 && trk_layer_pattern[layer-1] ==1 && layer != lminoise[0] && layer != lminoise[1] && layer != lminoise[2]) {
1032  noisemindist[3].CopyAMSPoint(pnoise[layer-1]);
1033  minoise[3] = noise_y[layer-1]; lminoise[3]=layer;
1034  }
1035  }
1036 
1037  return;
1038 }
1039 
1040 Double_t HitCloseTrack(int layer, int side, const char* type, AMSEventR* fEvent, TrTrackR* track, bool kMC, int TrackId){
1041 
1042  static Double_t ret;
1043  ret=999999;
1044 
1045  static unsigned int currEvent=0;
1046  static unsigned int currRun=0;
1047 
1048  static std::vector<TrClusterR*> TrClusUsed;
1049 
1050  static map<int, Double_t**> mapHM;
1051  static map<int, Double_t**> mapNH;
1052 
1053  static map<int, Double_t**>::const_iterator pos;
1054 
1055  if (currEvent!=fEvent->Event() || currRun!=fEvent->Run()){
1056 
1057  //clear maps and delete array pointed by pointers in map
1058  for (pos=mapHM.begin(); pos!=mapHM.end(); ++pos) {
1059  for (int ii=0; ii<2; ii++) {
1060  delete [] pos->second[ii];
1061  }
1062  delete [] pos->second;
1063  }
1064  for (pos=mapNH.begin(); pos!=mapNH.end(); ++pos) {
1065  for (int ii=0; ii<2; ii++) {
1066  delete [] pos->second[ii];
1067  }
1068  delete [] pos->second;
1069  }
1070  mapHM.clear();
1071  mapNH.clear();
1072 
1073  //build the vector with the used clusters
1074  BuildVectorUsedClusters(fEvent, TrClusUsed);
1075  }
1076 
1077  if (mapHM.find(TrackId)==mapHM.end()) {// not found
1078 
1079  //create and allocate array where to store temporarely the results for this particular TrackId
1080  Double_t** distHM_temp=0;
1081  Double_t** distNH_temp=0;
1082  distHM_temp = new Double_t*[2];
1083  distNH_temp = new Double_t*[2];
1084  for (int ii=0; ii<2; ii++) {
1085  distHM_temp[ii] = new Double_t[9];
1086  distNH_temp[ii] = new Double_t[9];
1087  }
1088 
1089  //compute the residuals
1090  BuildHitCloseTrack(fEvent, distHM_temp, distNH_temp, track, TrackId, TrClusUsed, kMC);
1091 
1092  //insert the array pointers into the map
1093  mapHM.insert(pair<int, Double_t**>(TrackId, distHM_temp));
1094  mapNH.insert(pair<int, Double_t**>(TrackId, distNH_temp));
1095 
1096  distHM_temp=0;
1097  distNH_temp=0;
1098 
1099  }
1100 
1101  //now the pointers exist (if they were not found were created and computed)
1102  if (!strcmp(type, "HM")) ret=mapHM.find(TrackId)->second[side][layer-1];
1103  else if (!strcmp(type, "NH")) ret=mapNH.find(TrackId)->second[side][layer-1];
1104  else ret=-999999;
1105 
1106  //save the current event and run number and the
1107  currEvent=fEvent->Event();
1108  currRun=fEvent->Run();
1109 
1110  return ret;
1111 }
1112 
1113 void BuildVectorUsedClusters(AMSEventR* fEvent, std::vector<TrClusterR*>& TrClusUsed){
1114 
1115  TrClusUsed.clear();
1116 
1117  for (int ii=0; ii<fEvent->nTrTrack(); ii++) {
1118  TrTrackR* trk = fEvent->pTrTrack(ii);
1119  for (int jj=0; jj<trk->NTrRecHit(); jj++) {
1120  TrRecHitR* rechit = trk->pTrRecHit(jj);
1121  if (rechit) {
1122  if (rechit->GetXCluster()) TrClusUsed.push_back(rechit->GetXCluster());
1123  if (rechit->GetYCluster()) TrClusUsed.push_back(rechit->GetYCluster());
1124  }
1125  }
1126  }
1127 
1128  // //---------- only for debug---------------------------------------------
1129  // static int used;
1130  // used=0;
1131  // for (int ii=0; ii<fEvent->nTrCluster(); ii++) {
1132  // static TrClusterR* clu;
1133  // clu = fEvent->pTrCluster(ii);
1134  // if (clu->Used()) {
1135  // used++;
1136  // }
1137  // }
1138 
1139  // if (used!=TrClusUsed.size()) printf("The number of used clusters doesn't match: %d %d!\n", used, TrClusUsed.size());
1140  // //-----------------------------------------------------------------------
1141 
1142  return;
1143 }
1144 
1145 void BuildHitCloseTrack(AMSEventR* fEvent, Double_t** distHM, Double_t** distNH, TrTrackR* track, int TrackId, std::vector<TrClusterR*> TrClusUsed, bool kMC){
1146 
1147  static std::vector <TrClusterR*>::iterator TrClusUsedIt;
1148 
1149  static std::vector<Double_t> distHMx_temp[9];
1150  static std::vector<Double_t> distHMy_temp[9];
1151 
1152  static std::vector<Double_t> distNHx_temp[9];
1153  static std::vector<Double_t> distNHy_temp[9];
1154 
1155  for (int ii=0; ii<9; ii++) {
1156  distHMx_temp[ii].clear();
1157  distHMy_temp[ii].clear();
1158 
1159  distNHx_temp[ii].clear();
1160  distNHy_temp[ii].clear();
1161  }
1162 
1163  // printf("New event\n");//only for debug
1164  for (int ii=0; ii<fEvent->nTrCluster(); ii++) {
1165  static TrClusterR* clu;
1166  clu = fEvent->pTrCluster(ii);
1167  TrClusUsedIt = std::find(TrClusUsed.begin(), TrClusUsed.end(), clu);
1168  if (TrClusUsedIt==TrClusUsed.end()) {//This means that the cluster wasn't found in the used ones
1169  static int layer;
1170  layer = clu->GetLayer();
1171  // printf("Layer = %d\n", layer);//only for debug
1172  if (TrackId>=0) {
1173  static std::vector<Double_t>* vectorx;
1174  static std::vector<Double_t>* vectory;
1175 
1176  static TkSens* sensor;
1177  if (sensor) delete sensor;
1178  if (kMC) sensor=new TkSens(true);
1179  else sensor=new TkSens(false);
1180  static AMSDir direction;
1181  static AMSPoint global;
1182  static AMSPoint local;
1183  static int mult;
1184  static Double_t coord1;
1185  static TrRecHitR* rechit;
1186  static int rhtkid;
1187  rhtkid=-1234;
1188  static int clutkid;
1189  clutkid=clu->GetTkId();
1190 #ifdef _B524_
1191  rechit = track->GetHitL(layer-1);
1192 #else
1193  rechit = track->GetHitLO(1+layer-1);
1194 #endif
1195  if (rechit) {
1196  //------------ this is done with the distance cluster used - cluster around, instead is better to use track - cluster around ----
1197  // rhtkid=rechit->GetTkId();
1198  // if (clutkid!=rhtkid) continue;//the hit used and the cluster are on different ladders so no problem
1199  // mult=rechit->GetMultiplicity(); // for X side
1200  // local=rechit->GetLocalCoordinate(mult);
1201  //-------------------------------------------------------------------------------------------------------------------------------
1202 #ifdef _B524_
1203  track->InterpolateLayer(layer-1, global, direction, TrackId);
1204 #else
1205  track->InterpolateLayerO(1+layer-1, global, direction, TrackId);
1206 #endif
1207  sensor->SetGlobal(global, direction);
1208  rhtkid=sensor->GetLadTkID();
1209  if (clutkid!=rhtkid) continue;//the point intercepted by the track and the cluster are on different ladders so no problem
1210  mult=sensor->GetMultIndex(); // for X side
1211  local=sensor->GetLaddCoo();
1212 
1213  vectorx=&distNHx_temp[layer-1];
1214  vectory=&distNHy_temp[layer-1];
1215  }
1216  else {
1217 #ifdef _B524_
1218  track->InterpolateLayer(layer-1, global, direction, TrackId);
1219 #else
1220  track->InterpolateLayerO(1+layer-1, global, direction, TrackId);
1221 #endif
1222  sensor->SetGlobal(global, direction);
1223  rhtkid=sensor->GetLadTkID();
1224  if (clutkid!=rhtkid) continue;//the point intercepted by the track and the cluster are on different ladders so no problem
1225  mult=sensor->GetMultIndex(); // for X side
1226  local=sensor->GetLaddCoo();
1227 
1228  vectorx=&distHMx_temp[layer-1];
1229  vectory=&distHMy_temp[layer-1];
1230  }
1231 
1232  coord1=clu->GetCoord(mult);
1233  static Double_t coord2;
1234  static Double_t dist;
1235  if (clu->GetSide()==0) {
1236  coord2=local.x();
1237  dist=1.e4*(coord2-coord1); // micron
1238  (*vectorx).push_back(dist);
1239  // printf("SizeX = %d\n", (*vectorx).size());
1240  }
1241  else {
1242  coord2=local.y();
1243  dist=1.e4*(coord2-coord1); // micron
1244  (*vectory).push_back(dist);
1245  // printf("SizeY = %d\n", (*vectory).size());
1246  }
1247  // printf("Residual = %f\n", dist);//only for debug
1248  }
1249  }
1250  }
1251 
1252  //search for the closer clusters
1253  for (int ii=0; ii<9; ii++) {
1254  static Double_t minore;
1255 
1256  minore=999999999;
1257  for (int jj=0; jj<(int)(distHMx_temp[ii].size()); jj++) if (fabs(distHMx_temp[ii][jj])<minore) minore=distHMx_temp[ii][jj];
1258  distHM[0][ii]=minore;
1259 
1260  minore=999999999;
1261  for (int jj=0; jj<(int)(distHMy_temp[ii].size()); jj++) if (fabs(distHMy_temp[ii][jj])<minore) minore=distHMy_temp[ii][jj];
1262  distHM[1][ii]=minore;
1263 
1264  minore=999999999;
1265  for (int jj=0; jj<(int)(distNHx_temp[ii].size()); jj++) if (fabs(distNHx_temp[ii][jj])<minore) minore=distNHx_temp[ii][jj];
1266  distNH[0][ii]=minore;
1267 
1268  minore=999999999;
1269  for (int jj=0; jj<(int)(distNHy_temp[ii].size()); jj++) if (fabs(distNHy_temp[ii][jj])<minore) minore=distNHy_temp[ii][jj];
1270  distNH[1][ii]=minore;
1271  }
1272 
1273  return;
1274 }
1275 
1276 //----------------------------------------------------------------------------------------------------------
1277 
1278 void PathOnB(TrTrackR* tr, int kDef, MagField* magfield, Double_t& B, Double_t& SenThetaX, Double_t& SenThetaY, Double_t& L, Double_t& pMDR, Double_t& pMDR_known, Double_t step, bool modulus){
1279 
1280  // Z of tracker layers
1281  static Double_t zlay[9]={53.060000, 29.228000, 25.212000, 1.698000, -2.318000, -25.212000, -29.228000, 158.920000, -135.882000};//layer esterni alla fine così se looppo solo su 7 piani automaticamente è l'inner
1282 
1283  B=0;
1284  SenThetaX=0;
1285  SenThetaY=0;
1286  L=0;
1287  pMDR=0;
1288  pMDR_known=0;
1289 
1290  if (!tr) {
1291  printf("PathOnB:: You passed an empty TrTrackR pointer... Returning 0's!");
1292  return;
1293  }
1294 
1295  static int ntrhits;
1296  ntrhits = tr->GetNhits();
1297 
1298  static Double_t start;
1299  start=-999999;
1300 
1301  static Double_t stop;
1302  stop=999999;
1303 
1304  for (int ii=0; ii<ntrhits; ii++) {
1305  static TrRecHitR* rh;
1306  rh=tr->pTrRecHit(ii);
1307  if (rh) {
1308  if (rh->GetYCluster()) {//only if there's a measurement of Y
1309  static int lay;
1310  lay=rh->GetLayer();
1311  static Double_t zed;
1312  zed=zlay[lay-1];
1313  // printf("%d) layer = %d, zed = %f\n", ii, lay, zed);//only for debug
1314  if (zed>start) start=zed;//the greatest
1315  if (zed<stop) stop=zed;//the smallest
1316  }
1317  }
1318  }
1319  // printf("Start = %f, Stop = %f\n", start, stop);
1320 
1321  static float x[3]={0, 0, 0};
1322  static float _B[3]={0, 0, 0};
1323 
1324 #define ONLYY
1325 
1326 #ifdef ONLYY
1327  int howmany=0;
1328  static AMSPoint inmiddleout[3];
1329  static Double_t zinmiddleout[3];
1330  zinmiddleout[0]=start;
1331  zinmiddleout[1]=(start+stop)/2;//a Z/2
1332  zinmiddleout[2]=stop;
1333  tr->Interpolate(3, zinmiddleout, inmiddleout, 0, 0, kDef);
1334  static Double_t DZ;
1335  static Double_t DY;
1336  static Double_t DX;
1337  static Double_t D_ZY;
1338  static Double_t D_ZX;
1339  DZ=fabs(inmiddleout[2].z()-inmiddleout[0].z());
1340  DY=fabs(inmiddleout[2].y()-inmiddleout[0].y());
1341  DX=fabs(inmiddleout[2].x()-inmiddleout[0].x());
1342  D_ZY=sqrt(DZ*DZ+DY*DY);
1343  D_ZX=sqrt(DZ*DZ+DX*DX);
1344  L=D_ZY;
1345  SenThetaX=fabs(DZ/D_ZX);//vertical means ThetaX=90 -> SenThetaX=1
1346  SenThetaY=fabs(DZ/D_ZY);//vertical means ThetaY=90 -> SenThetaY=1
1347  static Double_t y1;
1348  static Double_t y2;
1349  static Double_t y3;
1350  y1=inmiddleout[0].y();
1351  y2=inmiddleout[1].y();
1352  y3=inmiddleout[2].y();
1353  static Double_t SigmaY;
1354  static TF1* sigmay = new TF1("sigmay", "13.75-3.437*(x-1.0)+31.86*((x-1.0)*(x-1.0))", 0.0, 1.0);
1355  static Double_t sigmay0=sigmay->Eval(1.0)/1000000.0;//in m
1356  SigmaY=sigmay->Eval(SenThetaY)/1000000.0;//in m
1357  // printf("SigmaY = %f\n", SigmaY);//only for debug
1358  static Double_t SigmaS;
1359  static TF1* sigmas = new TF1("sigmas", "18.003-139.973*x+599.883*x*x-1496.24*x*x*x+2235.34*x*x*x*x-1976.75*x*x*x*x*x+954.921*x*x*x*x*x*x-194.295*x*x*x*x*x*x*x", 0.0, 1.0);
1360  static Double_t sigmas0=sigmas->Eval(1.0);
1361  static Double_t kost=sqrt(3.0/2.0);//from Gluckstern simple model
1362  SigmaS=kost*SigmaY*sigmas->Eval(SenThetaY)/sigmas0;
1363  // printf("SigmaS = %f\n", SigmaS);//only for debug
1364 #endif
1365 
1366  static AMSPoint pvec[10000];
1367  static Double_t zpl[10000];
1368  static int nz;
1369  nz=0;
1370  for (Double_t ii=start; ii>=(stop-1); ii-=step) {
1371  zpl[nz]=ii;
1372  nz++;
1373  }
1374  tr->Interpolate(nz, zpl, pvec, 0, 0, kDef);
1375  static int index;
1376  index=0;
1377 
1378  for (Double_t ii=start+step; ii>=stop; ii-=step) {
1379  static AMSPoint global;
1380  static AMSPoint global_old;
1381  index++;
1382  global=pvec[index];
1383  global_old=pvec[index-1];
1384 #ifdef ONLYY
1385  x[0]=global.x();
1386  x[1]=global.y();
1387  x[2]=global.z();
1388  if (magfield) magfield->GuFld(x, _B);
1389  if (modulus) B+=fabs(_B[0]);
1390  else B+=_B[0];//va bene con il segno perchè se ci fosse una zona con campo nell'altro verso questo "diminuirebbe" la curvatura
1391  howmany++;
1392 #else
1393  x[0]=global.x();
1394  x[1]=global.y();
1395  x[2]=global.z();
1396  static float dx[3]={0, 0, 0};
1397  dx[0]=x[0]-global_old.x();
1398  dx[1]=x[1]-global_old.y();
1399  dx[2]=x[2]-global_old.z();
1400  if (magfield) magfield->GuFld(x, _B);
1401  static Double_t vec[3];
1402  //---------------------------sbagliato----------------------------
1403  // if (_B[0]) {
1404  // vec[1]=-dx[2];
1405  // vec[2]=dx[1];
1406  // }
1407  // if (_B[1]) {
1408  // vec[0]=dx[2];
1409  // vec[2]-=dx[0];
1410  // }
1411  // if (_B[2]) {
1412  // vec[0]-=dx[1];
1413  // vec[1]+=dx[0];
1414  // }
1415  // B+=vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2];
1416  //---------------------------sbagliato-----------------------------
1417  vec[0]=_B[1]*dx[2]*dx[2] - _B[2]*dx[1]*dx[1];
1418  vec[1]=_B[2]*dx[0]*dx[0] - _B[0]*dx[2]*dx[2];
1419  vec[2]=_B[0]*dx[1]*dx[1] - _B[1]*dx[0]*dx[0];
1420  BL2+=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
1421 #endif
1422 #ifndef ONLYY
1423  // printf("%f) %f %f %f, %f %f %f -> %f %f\n", ii, x[0], x[1], x[2], B[0], B[1], B[2], intB, intBL2);//only for debug
1424 #else
1425  // printf("%f) %f %f %f, %f\n", ii, x[0], x[1], x[2], _B[0]);//only for debug
1426 #endif
1427  }
1428 #ifdef ONLYY
1429  // printf("B total = %f 0.1*T\n", B);
1430  B/=howmany;
1431  // printf("B mean = %f 0.1*T\n", B);
1432  pMDR=0.3*B*L*L/(8.0*SenThetaX*SigmaS);
1433  //pMDR = L^2qB / 8 sigmas sinthetax sin^2thetay.
1434  //0.3 is the factor conversion from Kg m/s / C*T to GeV/c / T
1435  //L here is DZ/SenThetaY and so the sinthetay is already included.
1436  pMDR_known=0.3*B*DZ*DZ/(8.0*SenThetaX*sigmay0*sqrt(3.0/2.0));
1437  //pMDR = L^2qB / 8 sigmas sinthetax sin^2thetay.
1438  //0.3 is the factor conversion from Kg m/s / C*T to GeV/c / T
1439  //sigmay0 is sigmay at sinthetay=1.0
1440  //sigmas=sqrt(3/2)*sigmay from Gluckstern simple model
1441 #endif
1442 
1443  B/=10;//to have in T instead of 0.1*T
1444  // printf("B (T) = %f\n", B);
1445  L/=100;//to have in m instead of cm
1446  pMDR/=(10*100*100);//to have in m^2*T instead of 0.1*cm^2*T
1447  pMDR_known/=(10*100*100);//to have in m^2*T instead of 0.1*cm^2*T
1448 
1449  L/=(2*sqrt(SigmaS*10000.0));
1450  // printf("L = %f\n", L);//only for debug
1451 
1452  // sleep(3);
1453 
1454  return;
1455 }
1456 
1457 Double_t LocalResidual(int layer, int side, bool kMC, AMSEventR* fEvent, TrTrackR* track, int TrackId){
1458 
1459  static Double_t res;
1460  res=-999999999;
1461 
1462  static unsigned int currEvent=0;
1463  static unsigned int currRun=0;
1464 
1465  static map<int, Double_t**> mapp;
1466 
1467  static map<int, Double_t**>::const_iterator pos;
1468 
1469  if (currEvent!=fEvent->Event() || currRun!=fEvent->Run()){
1470  //clear maps and delete array pointed by pointers in map
1471  for (pos=mapp.begin(); pos!=mapp.end(); ++pos) {
1472  for (int ii=0; ii<2; ii++) {
1473  delete [] pos->second[ii];
1474  }
1475  delete [] pos->second;
1476  }
1477  mapp.clear();
1478  }
1479 
1480  if (mapp.find(TrackId)==mapp.end()) {// not found
1481 
1482  //create and allocate array where to store temporarely the results for this particular TrackId
1483  Double_t** res_temp=0;
1484  res_temp = new Double_t*[2];
1485  for (int ii=0; ii<2; ii++) {
1486  res_temp[ii] = new Double_t[9];
1487  }
1488 
1489  //compute the residuals
1490  BuildResiduals(kMC, res_temp, track, TrackId);
1491 
1492  //insert the array pointers into the map
1493  mapp.insert(pair<int, Double_t**>(TrackId, res_temp));
1494  res_temp=0;
1495  }
1496 
1497  //now the pointers exist (if they were not found were created and computed)
1498  res=mapp.find(TrackId)->second[side][layer-1];
1499 
1500  //save the current event and run number
1501  currEvent=fEvent->Event();
1502  currRun=fEvent->Run();
1503 
1504  return res;
1505 }
1506 
1507 void BuildResiduals(bool kMC, Double_t** res, TrTrackR* track, int TrackId){
1508 
1509  for (int ii=0; ii<9; ii++) {//all over the layers
1510 
1511  static Double_t resy;
1512  static Double_t resx;
1513  resx=-999999999;
1514  resy=-999999999;
1515 
1516  TrRecHitR* rechit;
1517 #ifdef _B524_
1518  rechit = track->GetHitL(ii);
1519 #else
1520  rechit = track->GetHitLO(1+ii);
1521 #endif
1522 
1523  static TrClusterR* clux;
1524  static TrClusterR* cluy;
1525  clux=0;
1526  cluy=0;
1527 
1528  if (rechit) {
1529  clux=rechit->GetXCluster();
1530  cluy=rechit->GetYCluster();
1531  }
1532 
1533  //----residuals made by hand-------------------------------
1534  static AMSDir direction;
1535  static AMSPoint global;
1536  static AMSPoint local;
1537  static Double_t coord1;
1538  static Double_t coord2;
1539  static int mult;
1540  mult=0;
1541 
1542  static TkSens* sensor;
1543  if (sensor) delete sensor;
1544  if (kMC) sensor=new TkSens(true);
1545  else sensor=new TkSens(false);
1546  static int iltkid;
1547 
1548 #ifdef _B524_
1549  track->InterpolateLayer(ii, global, direction, TrackId);
1550 #else
1551  track->InterpolateLayerO(1+ii, global, direction, TrackId);
1552 #endif
1553  // printf("Layer = %d, Z = %f\n", ii+1, global.z());
1554  sensor->SetGlobal(global, direction);
1555  local=sensor->GetLaddCoo();
1556  mult=sensor->GetMultIndex(); // for X side
1557  iltkid=sensor->GetLadTkID();
1558 
1559  static int clutkid;
1560 
1561  if (clux) {
1562  clutkid=clux->GetTkId();
1563  if (clutkid!=iltkid) {
1564  clux=0;
1565  resx=-999999999;
1566  }
1567  else {
1568  coord1=clux->GetCoord(mult);
1569  coord2=local.x();
1570  resx=10000.0*(coord2-coord1); // micron
1571  }
1572  }
1573  if (cluy) {
1574  clutkid=cluy->GetTkId();
1575  if (clutkid!=iltkid) {
1576  cluy=0;
1577  resx=-999999999;
1578  }
1579  else {
1580  coord1=cluy->GetCoord(mult);
1581  coord2=local.y();
1582  resy=10000.0*(coord2-coord1); // micron
1583  // printf("Coo1 = %f, Coo2 = %f, ResY = %f\n", coord1, coord2, resy);
1584  }
1585  }
1586 
1587  res[0][ii]=resx;
1588  res[1][ii]=resy;
1589  }
1590 
1591  return;
1592 }
1593 
1594 #endif //#ifdef _WITHGBATCH_
1595 
1596 //----------------------------------------------------------------------------------------------------
Generated on Tue Dec 29 2015 23:02:50 for AMSDST by   doxygen 1.8.3