CaloPreSampler/src/CaloPreSampler.cpp

/**
 * \file CaloPreSampler.cpp
 * \author Emiliano Mocchiutti (2007/07/18)
 */
//
// headers
//
#include <CaloPreSampler.h>
//--------------------------------------
/**
 * Default constructor 
 */
CaloPreSampler::CaloPreSampler(){
  Clear();
}

void CaloPreSampler::SetCaloLevel2Pointer(CaloLevel2 *cl2p){
  if ( pcalo ) delete pcalo;
  pcalo = cl2p;
}

CaloPreSampler::CaloPreSampler(PamLevel2 *l2p){  
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  //
  // Default variables
  //
  event = new CaloLevel0();
  cstrip = new CaloStrip(false);
  //  c1 = new CaloLevel1();
  pcalo = new CaloLevel2();
  N = 4;
  NC = 22-N; 
  debug = false;
  sel = true;
  cont = false;
  emulate18 = true;
  usepl18x = false;
  simulation = false;
  withtrk = true;
  rigdefault = 50.;
  nox = false;
  noy = false;
  forcecalo = false;
  forcefitmode = -1;
  memset(mask,0,2*22*sizeof(Int_t));
  //
  Clear();
  //
  // loading magnetic field...
  //
  TrkLevel2 *trk = new TrkLevel2();
  GL_PARAM *q4 = new GL_PARAM();
  TSQLServer *dbc = 0;
  TString host = "mysql://localhost/pamelaprod";
  TString user = "anonymous";
  TString psw = "";
  const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
  const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
  const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
  if ( !pamdbhost ) pamdbhost = "";
  if ( !pamdbuser ) pamdbuser = "";
  if ( !pamdbpsw ) pamdbpsw = "";
  if ( strcmp(pamdbhost,"") ) host = pamdbhost;
  if ( strcmp(pamdbuser,"") ) user = pamdbuser;
  if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
  dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
  //
  q4->Query_GL_PARAM(1,1,dbc);
  printf(" Reading magnetic field maps at %s\n",(q4->PATH+q4->NAME).Data());
  trk->LoadField(q4->PATH+q4->NAME);
  //
}

CaloTrkVar* CaloPreSampler::GetCaloTrack(Int_t t){
  this->Process();
  if ( t < 0 ) return pcalo->GetCaloStoredTrack(t);
  if ( L2->GetTrack(t) ){
    return pcalo->GetCaloStoredTrack(L2->GetTrack(t)->GetTrkTrack()->GetSeqNo());
  };
  return NULL;
}

void CaloPreSampler::SetNoWpreSampler(Int_t n){
  if ( NC+n < 23 ){
    N = n;
  } else {
    printf(" ERROR! Calorimeter is made of 22 W planes\n");
    printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
    printf(" WARNING: using default values NWpre = 4, NWcalo = 18\n");
    NC = 18;
    N = 4;
  };
}

void CaloPreSampler::SetNoWcalo(Int_t n){
  if ( N+n < 23 ){
    NC = n;
  } else {
    printf(" ERROR! Calorimeter is made of 22 W planes\n");
    printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
    printf(" WARNING: using default values NWpre = 4, NWcalo = 18\n");
    NC = 18;
    N = 4;
  };
}

void CaloPreSampler::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
  this->SetNoWcalo(0);
  this->SetNoWpreSampler(0);
  if ( NoWpreSampler < NoWcalo ){
        this->SetNoWpreSampler(NoWpreSampler);
        this->SetNoWcalo(NoWcalo);
  } else {
        this->SetNoWcalo(NoWcalo);
        this->SetNoWpreSampler(NoWpreSampler);
  };
}

void CaloPreSampler::Clear(){
  //
  pcalo->Clear();
  //
}

void CaloPreSampler::Print(){
  //
  Process();
  //
  printf("========================================================================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
  printf(" debug                [debug flag]:.. %i\n",debug);
  printf(" simulation      [simulation flag]:.. %i\n",simulation);
  printf(" emulate18 [emulate dead plane 18]:.. %i\n",emulate18);
  printf(" selection mode                   :.. %i\n",sel);
  printf(" contamination mode               :.. %i\n",cont);
  printf(" pre-sampler planes               :.. %i\n",N);
  printf(" pcalo->qtot                      :.. %f\n",pcalo->qtot);
  printf(" pcalo->nstrip                    :.. %i\n",pcalo->nstrip);  
  if ( pcalo->ntrk() > 0 ){
    printf(" pcalo->track0->qtrack            :.. %f\n",pcalo->GetCaloTrkVar(0)->qtrack);  
    printf(" pcalo->track0->dX0l              :.. %f\n",pcalo->GetCaloTrkVar(0)->dX0l);  
  };
  printf("========================================================================\n");
  //
}

void CaloPreSampler::Delete(){
  Clear();
  delete pcalo;
  //delete this;
}


void CaloPreSampler::Process(){
  //  
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloPreSampler variables _NOT_ filled \n");
    return;
  };
  //
  // Clear structures used to communicate with fortran
  //
  event->ClearStructs();//ELENA
  if ( forcefitmode > 0 ){
    if ( forcefitmode != 1000 && forcefitmode != 1001 && forcefitmode != 1002 ){
      printf(" ERROR! forcefitmode=%i \n Use forcefitmode = 1000 for fit mode 0, 1001 fit mode 1, 1002 fit mode 3\n",forcefitmode);
      forcefitmode = -1;
    } else {
      event->clevel2->fmode[0] = forcefitmode;
      event->clevel2->fmode[1] = forcefitmode;
    };
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( debug ) printf(" I am here, we have orbital infos \n");
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || sel != ssel ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
      ssel = sel;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  // Some variables
  //
  Int_t S3 = 0;
  Int_t S2 = 0;
  Int_t S12 = 0;
  Int_t S11 = 0;
  Float_t tmptrigty = -1.; 
  Bool_t trackanyway = true;
  //  Float_t rigdefault = 50.;
  Bool_t hZn = true;
  //  Bool_t withtrk = true;
  Bool_t st = true;
  Int_t ntrkentry = 0;
  TrkLevel2 *trk = L2->GetTrkLevel2();
  //  Bool_t filled = false;
  //
  if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  // find out if we have trkseqno = -1, -2 or -3
  //
  Bool_t m1 = false;
  Bool_t m2 = false;
  Bool_t m3 = false;
  for (Int_t mm=0; mm < L2->GetCaloLevel2()->ntrk(); mm++ ){
    if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -1 ) m1 = true;
    if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -2 ) m2 = true;
    if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -3 ) m3 = true;
  };
  if ( !withtrk ) m3 = true;
  //
  // copy variables calculated during calibration process which is skipped here... this has been moved before this->Clear() for reprocessing purpose
  //
  event->clevel2->good = L2->GetCaloLevel2()->good;
  memcpy(event->clevel2->perr,L2->GetCaloLevel2()->perr,sizeof(L2->GetCaloLevel2()->perr));
  memcpy(event->clevel2->swerr,L2->GetCaloLevel2()->swerr,sizeof(L2->GetCaloLevel2()->swerr));
  memcpy(event->clevel2->crc,L2->GetCaloLevel2()->crc,sizeof(L2->GetCaloLevel2()->crc));
  event->clevel2->selftrigger = L2->GetCaloLevel2()->selftrigger;
  //
  //
  this->Clear(); // moved after all calls to L2->GetCaloLevel2() since pcalo could be the same as L2->GetCaloLevel2() !!!
  //
  //
  if ( debug ) printf(" Fill estrip matrix needed to calculate variables \n");
  //
  // Fill the estrip matrix
  //
  memset(event->clevel1->estrip, 0, 2*22*96*sizeof(Float_t));
  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;
  Float_t mip = 0.;
  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    //
    if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
    //
    // Mask x or y view if nox and/or noy are true (default false)
    //
    if ( mask[view][plane] ) mip = 0.;
    if ( nox && view == 0 ) mip = 0.; 
    if ( noy && view == 1 ) mip = 0.; 
    //
    if ( emulate18 && view == 0 && plane == 18 ) mip = 0.;
    //
    // Selection mode: fill the matrix only for plane < (22 - N) REMEMBER N = number of W planes to be used as presampler, ie if N = 2 then we want to use planes from 0 to 19
    //                 included so plane < (22 - 2)
    //
    if ( sel ){
      if ( plane < (22 - N) ){
        //
        if ( emulate18 && plane == (18 - N) && view == 0 ) mip = 0.;
        if ( plane >= NC ) mip = 0.;
        event->clevel1->estrip[strip][plane][view] = mip;
        //
      };
    };
    //
    // Contamination mode: fill the matrix only for planes from N to 22 but shift all planes up to the first one
    //
    if ( cont ){
      if ( plane >= N ){
        //
        if ( emulate18 && plane == (18 + N) && view == 0 ) mip = 0.;
        if ( (plane-N) >= NC ) mip = 0.;
        event->clevel1->estrip[strip][(plane-N)][view] = mip;
        //
      };
    };    
    //
  };
  //
  // if data comes from the simulation we must use mechanical alignment parameters (default is flight parameters)
  //
  if ( simulation ){
    cstrip->UseMechanicalAlig();
  };
  //
  // Set alignment parameter 
  //
  event->clevel1->xalig = cstrip->GetXalig();
  event->clevel1->yalig = cstrip->GetYalig();
  event->clevel1->zalig = cstrip->GetZalig();
  //
  event->clevel1->emin = 0.7;
  //
  // in case of the contamination mode we must play with the Z alignment in order to have the correct track in the calo since we have moved the planes up...
  //
  if ( cont ){
    if ( !(N%2) ){
      event->clevel1->reverse = 0; // if the number of planes is even we have taken away a full module no need to do anything strange...
      event->clevel1->zalig -= (N/2) * (8.09 + 10.09);
    } else {
      event->clevel1->reverse = 1; // if the number of planes is odd we have taken away half a module, we need to reverse silicon planes shifting
      event->clevel1->zalig -= ((N+1)/2) * 8.09 + ((N-1)/2) * 10.09;
    };
  };
  if ( debug ) printf(" xalig = %f \n",event->clevel1->xalig);
  if ( debug ) printf(" yalig = %f \n",event->clevel1->yalig);
  if ( debug ) printf(" zalig = %f \n",event->clevel1->zalig);
  //
  if ( debug ) printf(" Calculate variables as done in CaloCore, N = %i \n",N);
  //
  // Calculate variables
  //
  //
  // use only N W planes
  //
  //  event->clevel1->npla = 22-N;
  event->clevel1->npla = NC;
  //
  S3 = 0;
  S2 = 0;
  S12 = 0;
  S11 = 0;
  S3 = L2->GetTrigLevel2()->patterntrig[2];
  S2 = L2->GetTrigLevel2()->patterntrig[3];
  S12 = L2->GetTrigLevel2()->patterntrig[4];
  S11 = L2->GetTrigLevel2()->patterntrig[5];
  if ( L2->GetTrigLevel2()->patterntrig[1] & (1<<0) ) tmptrigty = 1.;
  if ( L2->GetTrigLevel2()->patterntrig[0] ) tmptrigty = 2.;
  if ( S3 || S2 || S12 || S11 )  tmptrigty = 0.;    
  if ( !(L2->GetTrigLevel2()->patterntrig[1] & (1<<0)) && !L2->GetTrigLevel2()->patterntrig[0] && !S3 && !S2 && !S12 && !S11 ) tmptrigty = 1.;
  event->clevel2->trigty = tmptrigty;
  //
  // do we have at least one track from the tracker? this check has been disabled
  //
  event->clevel1->good2 = 1;
  //
  // copy variables calculated during calibration process which is skipped here... this has to be moved before this->Clear() for reprocessing purpose
  //
//  event->clevel2->good = L2->GetCaloLevel2()->good;
//  memcpy(event->clevel2->perr,L2->GetCaloLevel2()->perr,sizeof(L2->GetCaloLevel2()->perr));
//  memcpy(event->clevel2->swerr,L2->GetCaloLevel2()->swerr,sizeof(L2->GetCaloLevel2()->swerr));
//  memcpy(event->clevel2->crc,L2->GetCaloLevel2()->crc,sizeof(L2->GetCaloLevel2()->crc));
//  event->clevel2->selftrigger = L2->GetCaloLevel2()->selftrigger;
  //
  // Calculate variables common to all tracks (qtot, nstrip, etc.)
  //
  if ( debug ) printf("1 Call GetCommonVar() \n");
  event->GetCommonVar();
  //
  // Fill common variables
  //
  if ( debug ) printf("1 Call FillCommonVar() \n");
  event->FillCommonVar(NULL,pcalo);
  //
  // Calculate variables related to tracks only if we have at least one track (from selftrigger and/or tracker)
  //
  ntrkentry = 0;
  //
  //  filled = false;
  //
  // Run over tracks (tracker or calorimeter )
  //
  if ( withtrk ){
    //
    for (Int_t nt=0; nt < trk->ntrk(); nt++){  
      //
      event->clevel1->good2 = 1;
      //
      TrkTrack *ptt = trk->GetStoredTrack(nt); 
      //
      event->clevel1->trkchi2 = 0;
      //
      // Copy the alpha vector in the input structure
      //
      for (Int_t e = 0; e < 5 ; e++){
        event->clevel1->al_p[e][0] = ptt->al[e];
      };          
      //
      // Get tracker related variables for this track
      //
      if ( debug ) printf("track %i Call GetTrkVar() \n",nt);
      event->GetTrkVar();
      if ( debug ) printf(" event->clevel2->dX0l %f \n",event->clevel2->dX0l);
      //
      // Save tracker track sequence number
      //        
      event->trkseqno = nt;
      //
      // Copy values in the class ca from the structure clevel2
      //
      if ( debug ) printf("track %i Call FillTrkVar() \n",nt);
      event->FillTrkVar(pcalo,ntrkentry);


      ntrkentry++;      
      //      filled = true;
      //
    }; // loop on all the tracks
  };
  // 
  // if no tracks found but there is the possibility to have a good track we should try to calculate anyway the track related variables using the calorimeter 
  // fit of the track (to be used for example when TRK is off due to any reason like IPM3/5 off).
  // here we make an event selection so it must be done very carefully...
  //
  // conditions are: 0) no track from the tracker 1) we have a track fit both in x and y 2) no problems with calo for this event 3) no selftrigger event
  //
  //  if ( trackanyway && !filled && event->clevel2->npcfit[0] >= 2 && event->clevel2->npcfit[1] >= 2 && event->clevel2->good != 0 && event->clevel2->trigty < 2. ){
  if ( (trackanyway && m3) || forcecalo ){
    if ( debug ) printf(" Event with a track not fitted by the tracker \n");
    //
    // Disable "track mode" in the fortran routine
    //
    event->clevel1->good2 = 0;
    event->clevel1->riginput = rigdefault;
    if ( debug ) printf(" Using as default rigidity: %f \n",event->clevel1->riginput);
    //
    // We have a selftrigger event to analyze.
    //
    for (Int_t e = 0; e < 5 ; e++){
      event->clevel1->al_p[e][0] = 0.;
      event->clevel1->al_p[e][1] = 0.;
    };
    event->clevel1->trkchi2 = 0;
    //
    if ( debug ) printf("-3 a Call GetTrkVar() \n");
    event->GetTrkVar();
    //
    // if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
    //
    if ( event->clevel1->good2 == 0 ) {
      //
      // In selftrigger mode the trkentry variable is set to -1
      //
      event->trkseqno = -3;
      //
      // Copy values in the class ca from the structure clevel2
      //
      if ( debug ) printf("-3 a Call FillTrkVar() \n");
      event->FillTrkVar(pcalo,ntrkentry);
      ntrkentry++;
      //      filled = true;
      //
    } else {
      if ( debug ) printf(" Selftrigger: problems with event \n");
    };
    //
  };
  //
  // Call high energy nuclei routine
  //
  //  if ( hZn && event->clevel2->trigty >= 2. ){
  if ( hZn && m2 ){
    if ( debug ) printf(" Calling selftrigger high energy nuclei routine \n");
    //
    // Disable "track mode" in the fortran routine
    //
    event->clevel1->good2 = 0;
    //
    // Set high energy nuclei flag to one
    // 
    event->clevel1->hzn = 1; 
    event->clevel1->riginput = rigdefault;
    //
    // We have a selftrigger event to analyze.
    //
    for (Int_t e = 0; e < 5 ; e++){
      event->clevel1->al_p[e][0] = 0.;
      event->clevel1->al_p[e][1] = 0.;
    };
    event->clevel1->trkchi2 = 0;
    //
    if ( debug ) printf("-2 a Call GetTrkVar() \n");
    event->GetTrkVar();
    //
    // if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
    //
    if ( event->clevel1->good2 == 0 ) {
      //
      // In selftrigger mode the trkentry variable is set to -1
      //
      event->trkseqno = -2;
      //
      // Copy values in the class ca from the structure clevel2
      //
      if ( debug ) printf("-2 a Call FillTrkVar() \n");
      event->FillTrkVar(pcalo,ntrkentry);
      ntrkentry++;
      //      filled = true;
      //
    } else {
      if ( debug ) printf(" Selftrigger: problems with event \n");
    };
    //
  };
  //
  // self trigger event
  //
  //  if ( st && event->clevel2->trigty >= 2. ){
  if ( st && m1 ){
    if ( debug ) printf(" Selftrigger event  \n");
    //
    // Disable "track mode" in the fortran routine
    //
    event->clevel1->good2 = 0;
    //
    // disable high enery nuclei flag;
    //
    event->clevel1->hzn = 0; 
    //
    // We have a selftrigger event to analyze.
    //
    for (Int_t e = 0; e < 5 ; e++){
      event->clevel1->al_p[e][0] = 0.;
      event->clevel1->al_p[e][1] = 0.;
    };
    event->clevel1->trkchi2 = 0;
    //
    if ( debug ) printf("-1 a Call GetTrkVar() \n");
    event->GetTrkVar();
    //
    // if we had no problem (clevel2->good = 0, NOTICE zero, not one in selftrigger mode!), fill and go on
    //
    if ( event->clevel1->good2 == 0 ) {
      //
      // In selftrigger mode the trkentry variable is set to -1
      //
      event->trkseqno = -1;
      //
      // Copy values in the class ca from the structure clevel2
      //
      if ( debug ) printf("-1 a Call FillTrkVar() \n");
      event->FillTrkVar(pcalo,ntrkentry);
      ntrkentry++;
      //      filled = true;
      //
    } else {
      if ( debug ) printf(" Selftrigger: problems with event \n");
    };
  };
//   //
//   // Clear structures used to communicate with fortran
//   //
//   event->ClearStructs();
// ELENA: moved @ beginning
  //
  //
  //
  if ( debug ) this->Print();
  if ( debug ) printf(" exit \n");
  //
}

//
// Method to add a calorimeter track, evaluated around a tracker track defined by a status vector.
// (can be used to evaluate the calorimeter track around an arbitrary axis, by setting the status vector with zero deflection )
// 
//
CaloTrkVar* CaloPreSampler::AddCaloTrkVar(float *al,int trktag){

    int ntrkentry = pcalo->ntrk();
    //
    for (Int_t nt=0; nt < ntrkentry; nt++){  
        if( pcalo->GetCaloTrkVar(nt)->trkseqno == trktag){
            cout << " CaloTrkVar* CaloPreSampler::AddCaloTrkVar(float *al,int trktag)"<<endl;
            cout << " --> trktag = "<<trktag<<" already defined "<<endl;
            return NULL;
        }
    }
    //
    event->clevel1->good2 = 1; //is a trk track
    event->clevel1->trkchi2 = 0;
    event->clevel1->hzn = 0; 
    //
    // Copy the alpha vector in the input structure
    //
    for (Int_t e = 0; e < 5 ; e++){
        event->clevel1->al_p[e][0] = al[e];
    };    
    //
    // Get tracker related variables for this track
    //
    if ( debug ) printf("track %i Call GetTrkVar() \n",trktag);
    event->GetTrkVar();
    if ( debug ) printf(" event->clevel2->dX0l %f \n",event->clevel2->dX0l);
    //
    // Save tracker track sequence number
    //  
    event->trkseqno = trktag;
    //
    // Copy values in the class ca from the structure clevel2
    //
    if ( debug ) printf("track %i Call FillTrkVar() \n",trktag);
    event->FillTrkVar(pcalo,ntrkentry);

    return pcalo->GetCaloTrkVar(ntrkentry);
    

};//ELENA
1	mocchiut	1.1	/**
2			* \file CaloPreSampler.cpp
3			* \author Emiliano Mocchiutti (2007/07/18)
4			*/
5			//
6			// headers
7			//
8			#include <CaloPreSampler.h>
9			//--------------------------------------
10			/**
11			* Default constructor
12			*/
13			CaloPreSampler::CaloPreSampler(){
14			Clear();
15	mocchiut	1.6	}
16	mocchiut	1.1
17	mocchiut	1.17	void CaloPreSampler::SetCaloLevel2Pointer(CaloLevel2 *cl2p){
18			if ( pcalo ) delete pcalo;
19			pcalo = cl2p;
20			}
21
22	mocchiut	1.1	CaloPreSampler::CaloPreSampler(PamLevel2 *l2p){
23			//
24			L2 = l2p;
25			//
26			if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
27			//
28			OBT = 0;
29			PKT = 0;
30			atime = 0;
31			//
32			// Default variables
33			//
34			event = new CaloLevel0();
35			cstrip = new CaloStrip(false);
36			// c1 = new CaloLevel1();
37			pcalo = new CaloLevel2();
38	mocchiut	1.6	N = 4;
39			NC = 22-N;
40	mocchiut	1.1	debug = false;
41			sel = true;
42			cont = false;
43			emulate18 = true;
44	mocchiut	1.16	usepl18x = false;
45	mocchiut	1.1	simulation = false;
46	mocchiut	1.4	withtrk = true;
47			rigdefault = 50.;
48	mocchiut	1.5	nox = false;
49			noy = false;
50	mocchiut	1.8	forcecalo = false;
51	mocchiut	1.13	forcefitmode = -1;
52	mocchiut	1.9	memset(mask,0,222sizeof(Int_t));
53	mocchiut	1.1	//
54			Clear();
55			//
56			// loading magnetic field...
57			//
58			TrkLevel2 *trk = new TrkLevel2();
59			GL_PARAM *q4 = new GL_PARAM();
60			TSQLServer *dbc = 0;
61			TString host = "mysql://localhost/pamelaprod";
62			TString user = "anonymous";
63			TString psw = "";
64			const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
65			const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
66			const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
67			if ( !pamdbhost ) pamdbhost = "";
68			if ( !pamdbuser ) pamdbuser = "";
69			if ( !pamdbpsw ) pamdbpsw = "";
70			if ( strcmp(pamdbhost,"") ) host = pamdbhost;
71			if ( strcmp(pamdbuser,"") ) user = pamdbuser;
72			if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
73			dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
74			//
75			q4->Query_GL_PARAM(1,1,dbc);
76			printf(" Reading magnetic field maps at %s\n",(q4->PATH+q4->NAME).Data());
77			trk->LoadField(q4->PATH+q4->NAME);
78			//
79	mocchiut	1.6	}
80
81	mocchiut	1.15	CaloTrkVar* CaloPreSampler::GetCaloTrack(Int_t t){
82			this->Process();
83			if ( t < 0 ) return pcalo->GetCaloStoredTrack(t);
84			if ( L2->GetTrack(t) ){
85			return pcalo->GetCaloStoredTrack(L2->GetTrack(t)->GetTrkTrack()->GetSeqNo());
86			};
87			return NULL;
88			}
89
90	mocchiut	1.6	void CaloPreSampler::SetNoWpreSampler(Int_t n){
91			if ( NC+n < 23 ){
92			N = n;
93			} else {
94			printf(" ERROR! Calorimeter is made of 22 W planes\n");
95			printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
96			printf(" WARNING: using default values NWpre = 4, NWcalo = 18\n");
97			NC = 18;
98			N = 4;
99			};
100			}
101
102			void CaloPreSampler::SetNoWcalo(Int_t n){
103			if ( N+n < 23 ){
104			NC = n;
105			} else {
106			printf(" ERROR! Calorimeter is made of 22 W planes\n");
107			printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
108			printf(" WARNING: using default values NWpre = 4, NWcalo = 18\n");
109			NC = 18;
110			N = 4;
111			};
112			}
113	mocchiut	1.1
114	mocchiut	1.11	void CaloPreSampler::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
115			this->SetNoWcalo(0);
116			this->SetNoWpreSampler(0);
117			if ( NoWpreSampler < NoWcalo ){
118			this->SetNoWpreSampler(NoWpreSampler);
119			this->SetNoWcalo(NoWcalo);
120			} else {
121			this->SetNoWcalo(NoWcalo);
122			this->SetNoWpreSampler(NoWpreSampler);
123			};
124			}
125
126	mocchiut	1.1	void CaloPreSampler::Clear(){
127			//
128			pcalo->Clear();
129			//
130	mocchiut	1.6	}
131	mocchiut	1.1
132			void CaloPreSampler::Print(){
133			//
134			Process();
135			//
136			printf("========================================================================\n");
137			printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
138			printf(" debug [debug flag]:.. %i\n",debug);
139			printf(" simulation [simulation flag]:.. %i\n",simulation);
140			printf(" emulate18 [emulate dead plane 18]:.. %i\n",emulate18);
141			printf(" selection mode :.. %i\n",sel);
142			printf(" contamination mode :.. %i\n",cont);
143			printf(" pre-sampler planes :.. %i\n",N);
144			printf(" pcalo->qtot :.. %f\n",pcalo->qtot);
145			printf(" pcalo->nstrip :.. %i\n",pcalo->nstrip);
146	mocchiut	1.2	if ( pcalo->ntrk() > 0 ){
147			printf(" pcalo->track0->qtrack :.. %f\n",pcalo->GetCaloTrkVar(0)->qtrack);
148			printf(" pcalo->track0->dX0l :.. %f\n",pcalo->GetCaloTrkVar(0)->dX0l);
149			};
150	mocchiut	1.1	printf("========================================================================\n");
151			//
152	mocchiut	1.6	}
153	mocchiut	1.1
154			void CaloPreSampler::Delete(){
155			Clear();
156			delete pcalo;
157			//delete this;
158	mocchiut	1.6	}
159	mocchiut	1.1
160
161			void CaloPreSampler::Process(){
162			//
163			if ( !L2 ){
164			printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
165			printf(" ERROR: CaloPreSampler variables _NOT_ filled \n");
166			return;
167			};
168			//
169	pam-fi	1.12	// Clear structures used to communicate with fortran
170			//
171			event->ClearStructs();//ELENA
172	mocchiut	1.13	if ( forcefitmode > 0 ){
173			if ( forcefitmode != 1000 && forcefitmode != 1001 && forcefitmode != 1002 ){
174			printf(" ERROR! forcefitmode=%i \n Use forcefitmode = 1000 for fit mode 0, 1001 fit mode 1, 1002 fit mode 3\n",forcefitmode);
175			forcefitmode = -1;
176			} else {
177			event->clevel2->fmode[0] = forcefitmode;
178			event->clevel2->fmode[1] = forcefitmode;
179			};
180			};
181	pam-fi	1.12	//
182	mocchiut	1.1	Bool_t newentry = false;
183			//
184			if ( L2->IsORB() ){
185	mocchiut	1.14	if ( debug ) printf(" I am here, we have orbital infos \n");
186	mocchiut	1.1	if ( L2->GetOrbitalInfo()->pkt_num != PKT \|\| L2->GetOrbitalInfo()->OBT != OBT \|\| L2->GetOrbitalInfo()->absTime != atime \|\| sel != ssel ){
187			newentry = true;
188			OBT = L2->GetOrbitalInfo()->OBT;
189			PKT = L2->GetOrbitalInfo()->pkt_num;
190			atime = L2->GetOrbitalInfo()->absTime;
191			ssel = sel;
192			};
193			} else {
194			newentry = true;
195			};
196			//
197			if ( !newentry ) return;
198			//
199			// Some variables
200			//
201			Int_t S3 = 0;
202			Int_t S2 = 0;
203			Int_t S12 = 0;
204			Int_t S11 = 0;
205			Float_t tmptrigty = -1.;
206			Bool_t trackanyway = true;
207	mocchiut	1.4	// Float_t rigdefault = 50.;
208	mocchiut	1.1	Bool_t hZn = true;
209	mocchiut	1.4	// Bool_t withtrk = true;
210	mocchiut	1.1	Bool_t st = true;
211			Int_t ntrkentry = 0;
212			TrkLevel2 *trk = L2->GetTrkLevel2();
213	mocchiut	1.19	// Bool_t filled = false;
214	mocchiut	1.1	//
215			if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
216			//
217			// find out if we have trkseqno = -1, -2 or -3
218			//
219			Bool_t m1 = false;
220			Bool_t m2 = false;
221			Bool_t m3 = false;
222			for (Int_t mm=0; mm < L2->GetCaloLevel2()->ntrk(); mm++ ){
223			if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -1 ) m1 = true;
224			if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -2 ) m2 = true;
225			if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -3 ) m3 = true;
226			};
227	mocchiut	1.4	if ( !withtrk ) m3 = true;
228	mocchiut	1.1	//
229	mocchiut	1.18	// copy variables calculated during calibration process which is skipped here... this has been moved before this->Clear() for reprocessing purpose
230			//
231			event->clevel2->good = L2->GetCaloLevel2()->good;
232			memcpy(event->clevel2->perr,L2->GetCaloLevel2()->perr,sizeof(L2->GetCaloLevel2()->perr));
233			memcpy(event->clevel2->swerr,L2->GetCaloLevel2()->swerr,sizeof(L2->GetCaloLevel2()->swerr));
234			memcpy(event->clevel2->crc,L2->GetCaloLevel2()->crc,sizeof(L2->GetCaloLevel2()->crc));
235			event->clevel2->selftrigger = L2->GetCaloLevel2()->selftrigger;
236			//
237			//
238			this->Clear(); // moved after all calls to L2->GetCaloLevel2() since pcalo could be the same as L2->GetCaloLevel2() !!!
239			//
240			//
241	mocchiut	1.1	if ( debug ) printf(" Fill estrip matrix needed to calculate variables \n");
242			//
243			// Fill the estrip matrix
244			//
245			memset(event->clevel1->estrip, 0, 22296*sizeof(Float_t));
246			Int_t view = 0;
247			Int_t plane = 0;
248			Int_t strip = 0;
249			Float_t mip = 0.;
250			for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
251			//
252			mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
253			//
254	mocchiut	1.16	if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
255			//
256	mocchiut	1.5	// Mask x or y view if nox and/or noy are true (default false)
257			//
258	mocchiut	1.9	if ( mask[view][plane] ) mip = 0.;
259	mocchiut	1.5	if ( nox && view == 0 ) mip = 0.;
260			if ( noy && view == 1 ) mip = 0.;
261			//
262	mocchiut	1.10	if ( emulate18 && view == 0 && plane == 18 ) mip = 0.;
263			//
264	mocchiut	1.1	// Selection mode: fill the matrix only for plane < (22 - N) REMEMBER N = number of W planes to be used as presampler, ie if N = 2 then we want to use planes from 0 to 19
265			// included so plane < (22 - 2)
266			//
267			if ( sel ){
268			if ( plane < (22 - N) ){
269			//
270	mocchiut	1.8	if ( emulate18 && plane == (18 - N) && view == 0 ) mip = 0.;
271	mocchiut	1.6	if ( plane >= NC ) mip = 0.;
272	mocchiut	1.1	event->clevel1->estrip[strip][plane][view] = mip;
273			//
274			};
275			};
276			//
277			// Contamination mode: fill the matrix only for planes from N to 22 but shift all planes up to the first one
278			//
279			if ( cont ){
280			if ( plane >= N ){
281			//
282	mocchiut	1.8	if ( emulate18 && plane == (18 + N) && view == 0 ) mip = 0.;
283	mocchiut	1.6	if ( (plane-N) >= NC ) mip = 0.;
284	mocchiut	1.1	event->clevel1->estrip[strip][(plane-N)][view] = mip;
285			//
286			};
287			};
288			//
289			};
290			//
291			// if data comes from the simulation we must use mechanical alignment parameters (default is flight parameters)
292			//
293			if ( simulation ){
294			cstrip->UseMechanicalAlig();
295			};
296			//
297			// Set alignment parameter
298			//
299			event->clevel1->xalig = cstrip->GetXalig();
300			event->clevel1->yalig = cstrip->GetYalig();
301			event->clevel1->zalig = cstrip->GetZalig();
302			//
303			event->clevel1->emin = 0.7;
304			//
305			// in case of the contamination mode we must play with the Z alignment in order to have the correct track in the calo since we have moved the planes up...
306			//
307			if ( cont ){
308			if ( !(N%2) ){
309			event->clevel1->reverse = 0; // if the number of planes is even we have taken away a full module no need to do anything strange...
310			event->clevel1->zalig -= (N/2) * (8.09 + 10.09);
311			} else {
312			event->clevel1->reverse = 1; // if the number of planes is odd we have taken away half a module, we need to reverse silicon planes shifting
313			event->clevel1->zalig -= ((N+1)/2) * 8.09 + ((N-1)/2) * 10.09;
314			};
315			};
316			if ( debug ) printf(" xalig = %f \n",event->clevel1->xalig);
317			if ( debug ) printf(" yalig = %f \n",event->clevel1->yalig);
318			if ( debug ) printf(" zalig = %f \n",event->clevel1->zalig);
319			//
320			if ( debug ) printf(" Calculate variables as done in CaloCore, N = %i \n",N);
321			//
322			// Calculate variables
323			//
324			//
325			// use only N W planes
326			//
327	mocchiut	1.7	// event->clevel1->npla = 22-N;
328			event->clevel1->npla = NC;
329	mocchiut	1.1	//
330			S3 = 0;
331			S2 = 0;
332			S12 = 0;
333			S11 = 0;
334			S3 = L2->GetTrigLevel2()->patterntrig[2];
335			S2 = L2->GetTrigLevel2()->patterntrig[3];
336			S12 = L2->GetTrigLevel2()->patterntrig[4];
337			S11 = L2->GetTrigLevel2()->patterntrig[5];
338			if ( L2->GetTrigLevel2()->patterntrig[1] & (1<<0) ) tmptrigty = 1.;
339			if ( L2->GetTrigLevel2()->patterntrig[0] ) tmptrigty = 2.;
340			if ( S3 \|\| S2 \|\| S12 \|\| S11 ) tmptrigty = 0.;
341			if ( !(L2->GetTrigLevel2()->patterntrig[1] & (1<<0)) && !L2->GetTrigLevel2()->patterntrig[0] && !S3 && !S2 && !S12 && !S11 ) tmptrigty = 1.;
342			event->clevel2->trigty = tmptrigty;
343			//
344			// do we have at least one track from the tracker? this check has been disabled
345			//
346			event->clevel1->good2 = 1;
347			//
348	mocchiut	1.18	// copy variables calculated during calibration process which is skipped here... this has to be moved before this->Clear() for reprocessing purpose
349	mocchiut	1.3	//
350	mocchiut	1.18	// event->clevel2->good = L2->GetCaloLevel2()->good;
351			// memcpy(event->clevel2->perr,L2->GetCaloLevel2()->perr,sizeof(L2->GetCaloLevel2()->perr));
352			// memcpy(event->clevel2->swerr,L2->GetCaloLevel2()->swerr,sizeof(L2->GetCaloLevel2()->swerr));
353			// memcpy(event->clevel2->crc,L2->GetCaloLevel2()->crc,sizeof(L2->GetCaloLevel2()->crc));
354			// event->clevel2->selftrigger = L2->GetCaloLevel2()->selftrigger;
355	mocchiut	1.3	//
356	mocchiut	1.1	// Calculate variables common to all tracks (qtot, nstrip, etc.)
357			//
358			if ( debug ) printf("1 Call GetCommonVar() \n");
359			event->GetCommonVar();
360			//
361			// Fill common variables
362			//
363			if ( debug ) printf("1 Call FillCommonVar() \n");
364			event->FillCommonVar(NULL,pcalo);
365			//
366			// Calculate variables related to tracks only if we have at least one track (from selftrigger and/or tracker)
367			//
368			ntrkentry = 0;
369			//
370	mocchiut	1.19	// filled = false;
371	mocchiut	1.1	//
372			// Run over tracks (tracker or calorimeter )
373			//
374			if ( withtrk ){
375			//
376			for (Int_t nt=0; nt < trk->ntrk(); nt++){
377			//
378			event->clevel1->good2 = 1;
379			//
380			TrkTrack *ptt = trk->GetStoredTrack(nt);
381			//
382			event->clevel1->trkchi2 = 0;
383			//
384			// Copy the alpha vector in the input structure
385			//
386			for (Int_t e = 0; e < 5 ; e++){
387			event->clevel1->al_p[e][0] = ptt->al[e];
388			};
389			//
390			// Get tracker related variables for this track
391			//
392			if ( debug ) printf("track %i Call GetTrkVar() \n",nt);
393			event->GetTrkVar();
394			if ( debug ) printf(" event->clevel2->dX0l %f \n",event->clevel2->dX0l);
395			//
396			// Save tracker track sequence number
397			//
398			event->trkseqno = nt;
399			//
400			// Copy values in the class ca from the structure clevel2
401			//
402			if ( debug ) printf("track %i Call FillTrkVar() \n",nt);
403			event->FillTrkVar(pcalo,ntrkentry);
404
405
406			ntrkentry++;
407	mocchiut	1.19	// filled = true;
408	mocchiut	1.1	//
409			}; // loop on all the tracks
410			};
411			//
412			// if no tracks found but there is the possibility to have a good track we should try to calculate anyway the track related variables using the calorimeter
413			// fit of the track (to be used for example when TRK is off due to any reason like IPM3/5 off).
414			// here we make an event selection so it must be done very carefully...
415			//
416			// conditions are: 0) no track from the tracker 1) we have a track fit both in x and y 2) no problems with calo for this event 3) no selftrigger event
417			//
418			// if ( trackanyway && !filled && event->clevel2->npcfit[0] >= 2 && event->clevel2->npcfit[1] >= 2 && event->clevel2->good != 0 && event->clevel2->trigty < 2. ){
419	mocchiut	1.8	if ( (trackanyway && m3) \|\| forcecalo ){
420	mocchiut	1.1	if ( debug ) printf(" Event with a track not fitted by the tracker \n");
421			//
422			// Disable "track mode" in the fortran routine
423			//
424			event->clevel1->good2 = 0;
425			event->clevel1->riginput = rigdefault;
426			if ( debug ) printf(" Using as default rigidity: %f \n",event->clevel1->riginput);
427			//
428			// We have a selftrigger event to analyze.
429			//
430			for (Int_t e = 0; e < 5 ; e++){
431			event->clevel1->al_p[e][0] = 0.;
432			event->clevel1->al_p[e][1] = 0.;
433			};
434			event->clevel1->trkchi2 = 0;
435			//
436			if ( debug ) printf("-3 a Call GetTrkVar() \n");
437			event->GetTrkVar();
438			//
439			// if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
440			//
441			if ( event->clevel1->good2 == 0 ) {
442			//
443			// In selftrigger mode the trkentry variable is set to -1
444			//
445			event->trkseqno = -3;
446			//
447			// Copy values in the class ca from the structure clevel2
448			//
449			if ( debug ) printf("-3 a Call FillTrkVar() \n");
450			event->FillTrkVar(pcalo,ntrkentry);
451			ntrkentry++;
452	mocchiut	1.19	// filled = true;
453	mocchiut	1.1	//
454			} else {
455			if ( debug ) printf(" Selftrigger: problems with event \n");
456			};
457			//
458			};
459			//
460			// Call high energy nuclei routine
461			//
462			// if ( hZn && event->clevel2->trigty >= 2. ){
463			if ( hZn && m2 ){
464			if ( debug ) printf(" Calling selftrigger high energy nuclei routine \n");
465			//
466			// Disable "track mode" in the fortran routine
467			//
468			event->clevel1->good2 = 0;
469			//
470			// Set high energy nuclei flag to one
471			//
472			event->clevel1->hzn = 1;
473			event->clevel1->riginput = rigdefault;
474			//
475			// We have a selftrigger event to analyze.
476			//
477			for (Int_t e = 0; e < 5 ; e++){
478			event->clevel1->al_p[e][0] = 0.;
479			event->clevel1->al_p[e][1] = 0.;
480			};
481			event->clevel1->trkchi2 = 0;
482			//
483			if ( debug ) printf("-2 a Call GetTrkVar() \n");
484			event->GetTrkVar();
485			//
486			// if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
487			//
488			if ( event->clevel1->good2 == 0 ) {
489			//
490			// In selftrigger mode the trkentry variable is set to -1
491			//
492			event->trkseqno = -2;
493			//
494			// Copy values in the class ca from the structure clevel2
495			//
496			if ( debug ) printf("-2 a Call FillTrkVar() \n");
497			event->FillTrkVar(pcalo,ntrkentry);
498			ntrkentry++;
499	mocchiut	1.19	// filled = true;
500	mocchiut	1.1	//
501			} else {
502			if ( debug ) printf(" Selftrigger: problems with event \n");
503			};
504			//
505			};
506			//
507			// self trigger event
508			//
509			// if ( st && event->clevel2->trigty >= 2. ){
510			if ( st && m1 ){
511			if ( debug ) printf(" Selftrigger event \n");
512			//
513			// Disable "track mode" in the fortran routine
514			//
515			event->clevel1->good2 = 0;
516			//
517			// disable high enery nuclei flag;
518			//
519			event->clevel1->hzn = 0;
520			//
521			// We have a selftrigger event to analyze.
522			//
523			for (Int_t e = 0; e < 5 ; e++){
524			event->clevel1->al_p[e][0] = 0.;
525			event->clevel1->al_p[e][1] = 0.;
526			};
527			event->clevel1->trkchi2 = 0;
528			//
529			if ( debug ) printf("-1 a Call GetTrkVar() \n");
530			event->GetTrkVar();
531			//
532			// if we had no problem (clevel2->good = 0, NOTICE zero, not one in selftrigger mode!), fill and go on
533			//
534			if ( event->clevel1->good2 == 0 ) {
535			//
536			// In selftrigger mode the trkentry variable is set to -1
537			//
538			event->trkseqno = -1;
539			//
540			// Copy values in the class ca from the structure clevel2
541			//
542			if ( debug ) printf("-1 a Call FillTrkVar() \n");
543			event->FillTrkVar(pcalo,ntrkentry);
544			ntrkentry++;
545	mocchiut	1.19	// filled = true;
546	mocchiut	1.1	//
547			} else {
548			if ( debug ) printf(" Selftrigger: problems with event \n");
549			};
550			};
551	pam-fi	1.12	// //
552			// // Clear structures used to communicate with fortran
553			// //
554			// event->ClearStructs();
555			// ELENA: moved @ beginning
556	mocchiut	1.1	//
557			//
558			//
559			if ( debug ) this->Print();
560			if ( debug ) printf(" exit \n");
561			//
562	mocchiut	1.6	}
563	pam-fi	1.12
564			//
565			// Method to add a calorimeter track, evaluated around a tracker track defined by a status vector.
566			// (can be used to evaluate the calorimeter track around an arbitrary axis, by setting the status vector with zero deflection )
567			//
568			//
569			CaloTrkVar* CaloPreSampler::AddCaloTrkVar(float *al,int trktag){
570
571			int ntrkentry = pcalo->ntrk();
572			//
573			for (Int_t nt=0; nt < ntrkentry; nt++){
574			if( pcalo->GetCaloTrkVar(nt)->trkseqno == trktag){
575			cout << " CaloTrkVar* CaloPreSampler::AddCaloTrkVar(float *al,int trktag)"<<endl;
576			cout << " --> trktag = "<<trktag<<" already defined "<<endl;
577			return NULL;
578			}
579			}
580			//
581			event->clevel1->good2 = 1; //is a trk track
582			event->clevel1->trkchi2 = 0;
583			event->clevel1->hzn = 0;
584			//
585			// Copy the alpha vector in the input structure
586			//
587			for (Int_t e = 0; e < 5 ; e++){
588			event->clevel1->al_p[e][0] = al[e];
589			};
590			//
591			// Get tracker related variables for this track
592			//
593			if ( debug ) printf("track %i Call GetTrkVar() \n",trktag);
594			event->GetTrkVar();
595			if ( debug ) printf(" event->clevel2->dX0l %f \n",event->clevel2->dX0l);
596			//
597			// Save tracker track sequence number
598			//
599			event->trkseqno = trktag;
600			//
601			// Copy values in the class ca from the structure clevel2
602			//
603			if ( debug ) printf("track %i Call FillTrkVar() \n",trktag);
604			event->FillTrkVar(pcalo,ntrkentry);
605
606			return pcalo->GetCaloTrkVar(ntrkentry);
607
608
609			};//ELENA