EXP2003 » EXP2003: Express analysis

Int_t GetClusterMWPC(TClonesArray *data);

void MWPCprojection();

Float_t GetPosition(Float_t wire, Float_t wireStep,Float_t planeOffset);

void fillSi(TClonesArray *data,Float_t* amp,Float_t* time,Float_t *par1,Float_t *par2,Int_t *multiplicity, Float_t ampThreshold,Float_t timeThreshold=0);

void fillarrayCsI(TClonesArray *data,Float_t* amp,Float_t* time,Float_t *par1,Float_t *par2);

void processCsI_cal();

void zeroVars();

void fillF5(TClonesArray *data);

void fillF3(TClonesArray *data);

void fillMWPC(TClonesArray *data,Float_t *wire);

void readPar(TString fileName,Float_t *par1,Float_t *par2,Int_t size=16);

Int_t GetClusterSi(TClonesArray *data);

Int_t GetClusterMult(TClonesArray *data);

ERQTelescopeCsIDigi* processCsI(TClonesArray *data);

void fillCsI(ERQTelescopeCsIDigi *data);

//outtree vars

Int_t trigger; 

Float_t F5,tF5,F3,tF3;

Float_t wirex1,wirex2,wirey1,wirey2;

Float_t tMWPC[4];

Float_t fXt,fYt;

Float_t x1c, y1c, x2c, y2c;

Int_t nCsI;

Float_t aCsI,tCsI,aCsI_cal;

Float_t arCsI[16],trCsI[16];

Float_t DSD_X[32],DSD_Y[32];

Float_t tDSD_X[32],tDSD_Y[32];

Float_t SQ20_1[16],SSD1[16],SSD_V1[16];

Float_t tSQ20_1[16],tSSD1[16],tSSD_V1[16];

Float_t SQ20_2[16],SSD2[16],SSD_V2[16];

Float_t tSQ20_2[16],tSSD2[16],tSSD_V2[16];

Float_t SQ20_3[16],SSD3[16],SSD_V3[16];

Float_t tSQ20_3[16],tSSD3[16],tSSD_V3[16];

Float_t SQ20_4[16],SSD4[16],SSD_V4[16];

Float_t tSQ20_4[16],tSSD4[16],tSSD_V4[16];

// Sipars

Float_t pSQ201_1[16],pSQ201_2[16];

Float_t pSQ202_1[16],pSQ202_2[16];

Float_t pSQ203_1[16],pSQ203_2[16];

Float_t pSQ204_1[16],pSQ204_2[16];

Float_t pSSD1_1[16],pSSD1_2[16];

Float_t pSSD2_1[16],pSSD2_2[16];

Float_t pSSD3_1[16],pSSD3_2[16];

Float_t pSSD4_1[16],pSSD4_2[16];

Float_t pSSD_V1_1[16],pSSD_V1_2[16];

Float_t pSSD_V2_1[16],pSSD_V2_2[16];

Float_t pSSD_V3_1[16],pSSD_V3_2[16];

Float_t pSSD_V4_1[16],pSSD_V4_2[16];

Float_t pDSD_X1[32],pDSD_X2[32];

Float_t pDSD_Y1[32],pDSD_Y2[32];

Float_t pCsI_1[16],pCsI_2[16];

Float_t *nullPtr = NULL;

//

Int_t mult20_1,mult1_1,multv_1;

Int_t mult20_2,mult1_2,multv_2;

Int_t mult20_3,mult1_3,multv_3;

Int_t mult20_4,mult1_4,multv_4;

Int_t multc_x,multc_y,multCsI;

void convert(TString InFile, TString OutFile) {

  TChain *ch = new TChain("er");

  ch->Add(InFile.Data());

  // ch->Add("/media/ivan/data/exp1904/digi/h7/h7_ct_40*.root");

  cout << "Found " << ch->GetEntries() << " entries" << endl;

  // readPar("SSD_20u_1_cal",pSQ201_1,pSQ201_2);

  // readPar("SSD_20u_2_cal",pSQ202_1,pSQ202_2);

  // readPar("SSD_20u_3_cal",pSQ203_1,pSQ203_2);

  // readPar("SSD_20u_4_cal",pSQ204_1,pSQ204_2);

  // readPar("SSD_1m_1_cal",pSSD1_1,pSSD1_2);

  // readPar("SSD_1m_2_cal",pSSD2_1,pSSD2_2);

  // readPar("SSD_1m_3_cal",pSSD3_1,pSSD3_2);

  // readPar("SSD_1m_4_cal",pSSD4_1,pSSD4_2);

  // readPar("SSD_1m_1_v_cal",pSSD_V1_1,pSSD_V1_2);

  // readPar("SSD_1m_2_v_cal",pSSD_V2_1,pSSD_V2_2);  

  // readPar("SSD_1m_3_v_cal",pSSD_V3_1,pSSD_V3_2);

  // readPar("SSD_1m_4_v_cal",pSSD_V4_1,pSSD_V4_2);

  // readPar("DSSD_X",pDSD_X1,pDSD_X2,32);

  // readPar("DSSD_Y",pDSD_Y1,pDSD_Y2,32);

  // readPar("empty16",pCsI_1,pCsI_2);

//--------------------------------------------------------------------------------

  ERBeamTimeEventHeader* header = new ERBeamTimeEventHeader();

  TClonesArray *v_DSDX_C = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_DSDY_C = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_CsI = new TClonesArray("ERQTelescopeCsIDigi");

  TClonesArray *v_F5 = new TClonesArray("ERBeamDetTOFDigi");

  TClonesArray *v_F3 = new TClonesArray("ERBeamDetTOFDigi");  

  TClonesArray *v_MWPCx1 = new TClonesArray("ERBeamDetMWPCDigi");

  TClonesArray *v_MWPCy1 = new TClonesArray("ERBeamDetMWPCDigi");  

  TClonesArray *v_MWPCx2 = new TClonesArray("ERBeamDetMWPCDigi");

  TClonesArray *v_MWPCy2 = new TClonesArray("ERBeamDetMWPCDigi"); 

  TClonesArray *v_SSD_V_1 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD20_1 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD_1 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD_V_2 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD20_2 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD_2 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD_V_3 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD20_3 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD_3 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD_V_4 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD20_4 = new TClonesArray("ERQTelescopeSiDigi");

  TClonesArray *v_SSD_4 = new TClonesArray("ERQTelescopeSiDigi");

  // setbranchAdress

  ch->SetBranchAddress("BeamDetToFDigi2",&v_F5);

  ch->SetBranchAddress("BeamDetToFDigi1",&v_F3);

  ch->SetBranchAddress("BeamDetMWPCDigiX1",&v_MWPCx1);

  ch->SetBranchAddress("BeamDetMWPCDigiX2",&v_MWPCx2);

  ch->SetBranchAddress("BeamDetMWPCDigiY1",&v_MWPCy1);

  ch->SetBranchAddress("BeamDetMWPCDigiY2",&v_MWPCy2);  

  ch->SetBranchAddress("EventHeader.",&header);

  TFile *f;

  TString cutName;

  // Creating outfile,outtree

  TFile *fw = new TFile(OutFile.Data(), "RECREATE");

  TTree *tw = new TTree("tree", "data");

  tw->Branch("trigger",&trigger,"trigger/I");

  tw->Branch("F5.",&F5,"F5/F");

  tw->Branch("tF5.",&tF5,"tF5/F");

  tw->Branch("F3.",&F3,"F3/F");

  tw->Branch("tF3.",&tF3,"tF3/F");

  tw->Branch("tMWPC.",&tMWPC,"tMWPC[4]/F");

  tw->Branch("wirex1.",&wirex1,"wirex1/F");

  tw->Branch("wirex2.",&wirex2,"wirex2/F");

  tw->Branch("wirey1.",&wirey1,"wirey1/F");

  tw->Branch("wirey2.",&wirey2,"wirey2/F");

  tw->Branch("x1c.",&x1c,"x1c/F"); // position in MWPC in mm

  tw->Branch("y1c.",&y1c,"y1c/F");

  tw->Branch("x2c.",&x2c,"x2c/F");

  tw->Branch("y2c.",&y2c,"y2c/F"); 

  tw->Branch("fXt.",&fXt,"fXt/F"); // beam profile at the target plane

  tw->Branch("fYt.",&fYt,"fYt/F"); 

  // for(Int_t nentry=0;nentry<ch->GetEntries();nentry++) {

  for(Int_t nentry=0;nentry<100;nentry++) {    

    if(nentry%1000000==0) cout << "#Event " << nentry << "#" << endl;

    ch->GetEntry(nentry);

    if(v_MWPCx1->GetEntriesFast()==0) continue;

    if(v_MWPCx2->GetEntriesFast()==0) continue;

    if(v_MWPCy1->GetEntriesFast()==0) continue;

    if(v_MWPCy2->GetEntriesFast()==0) continue;

    if (GetClusterMWPC(v_MWPCx1)!=1) continue;

    if (GetClusterMWPC(v_MWPCx2)!=1) continue;

    if (GetClusterMWPC(v_MWPCy1)!=1) continue;

    if (GetClusterMWPC(v_MWPCy2)!=1) continue;  

    trigger = header->GetTrigger();

    zeroVars();

  // beamdet

    fillF5(v_F5);

    fillF3(v_F3);

    tMWPC[0] = ((ERBeamDetMWPCDigi*)v_MWPCx1->At(0))->Time();

    tMWPC[1] = ((ERBeamDetMWPCDigi*)v_MWPCy1->At(0))->Time();

    tMWPC[2] = ((ERBeamDetMWPCDigi*)v_MWPCx2->At(0))->Time();

    tMWPC[3] = ((ERBeamDetMWPCDigi*)v_MWPCy2->At(0))->Time();

    fillMWPC(v_MWPCx1,&wirex1);

    fillMWPC(v_MWPCy1,&wirey1);

    fillMWPC(v_MWPCx2,&wirex2);

    fillMWPC(v_MWPCy2,&wirey2);

    MWPCprojection();

    tw->Fill();

  }

  fw->cd();

  tw->Write();

  fw->Close();

  return;

}

Int_t GetClusterMWPC(TClonesArray *data) {

  if (!data) return 0;

  Int_t entries = data->GetEntriesFast();

  if (entries<2) return entries;

  Int_t wire1, wire2;

  Int_t noclusters = 1;

  for (Int_t i = 1; i < entries; i++) {

    //check if entries are in specific order

    wire1 = ((ERBeamDetMWPCDigi*)data->At(i))->Channel();

    wire2 = ((ERBeamDetMWPCDigi*)data->At(i-1))->Channel();

    //todo number 32 is related to number of wires

    // and should be taken from Parameters

    if ( abs(wire1 - wire2) > 1 && abs(wire1 - wire2) < 32) {

      noclusters++;

    }

  }

  return noclusters;

}

void zeroVars() {

  // trigger = -1;

  F5 = 0;

  tF5 = 0;

  F3 = 0;

  tF3 = 0;

  wirex1 = -1;

  wirex2 = -1;

  wirey1 = -1;

  wirey2 = -1;

  for(Int_t i = 0;i<4;i++) {

    tMWPC[i] = 0;

  }

  x1c = -50;

  y1c = -50;

  x2c = -50;

  y2c = -50;

  fXt = -100;

  fYt = -100;

  nCsI = -1;

  aCsI = 0;

  tCsI = 0;

  for(Int_t i=0;i<16;i++) {

    arCsI[i] = 0;

    trCsI[i] = 0;

    SQ20_1[i] = 0;

    SSD1[i] = 0;

    SSD_V1[i] = -10;

    tSQ20_1[i] = 0;

    tSSD1[i] = 0;

    tSSD_V1[i] = -10;

    SQ20_2[i] = 0;

    SSD2[i] = 0;

    SSD_V2[i] = 0;

    tSQ20_2[i] = 0;

    tSSD2[i] = -10;

    tSSD_V2[i] = -10;

    SQ20_3[i] = 0;

    SSD3[i] = 0;

    SSD_V3[i] = 0;

    tSQ20_3[i] = 0;

    tSSD3[i] = -10;

    tSSD_V3[i] = -10;

    SQ20_4[i] = 0;

    SSD4[i] = 0;

    SSD_V4[i] = 0;

    tSQ20_4[i] = 0;

    tSSD4[i] = -10;

    tSSD_V4[i] = -10;

  }

  for(Int_t i=0;i<32;i++) {

    DSD_X[i] = 0;

    DSD_Y[i] = 0;

    tDSD_X[i] = 0;

    tDSD_Y[i] = 0;

  }

}

void fillF5(TClonesArray *data){

  if (!data) return;

  ERBeamDetTOFDigi *temp_F5 = ((ERBeamDetTOFDigi*)data->At(0));

  if(!temp_F5) return;

  F5 = temp_F5->Edep();

  tF5 = temp_F5->Time();

  return; 

}

void fillF3(TClonesArray *data){

  if (!data) return;

  ERBeamDetTOFDigi *temp_F3 = ((ERBeamDetTOFDigi*)data->At(0));

  if(!temp_F3) return;

  F3 = temp_F3->Edep();

  tF3 = temp_F3->Time();

  return; 

}

void fillMWPC(TClonesArray *data,Float_t *wire) {

  if (!data) return;

  *(wire) = ((ERBeamDetMWPCDigi*)data->At(0))->Channel() + ((ERBeamDetMWPCDigi*)data->At(data->GetEntriesFast()-1))->Channel();

  *(wire) = *(wire)/2; 

  return;

}

void MWPCprojection() {

  const Float_t fMWPCwireStepX1 = -1.25;

  const Float_t fMWPCwireStepY1 = 1.25;   //step between two wires

  const Float_t fMWPCwireStepX2 = -1.25;    //step between two wires

  const Float_t fMWPCwireStepY2 = 1.25;   //step between two wires

  const Float_t fMWPC1_X_offset = 0;

  const Float_t fMWPC1_Y_offset = 0;

  const Float_t fMWPC2_X_offset = 0;

  const Float_t fMWPC2_Y_offset = 0;

  const Float_t fMWPCz1 = -815.;  //z coordinate of the center of MWPC1

  const Float_t fMWPCz2 = -270.;  //z coordinate of the center of MWPC2

  Float_t xtc, ytc;

  //cluster multiplicity equal to 1

  x1c = GetPosition(wirex1, fMWPCwireStepX1, fMWPC1_X_offset);

  y1c = GetPosition(wirey1, fMWPCwireStepY1, fMWPC1_Y_offset);

  x2c = GetPosition(wirex2, fMWPCwireStepX2, fMWPC2_X_offset);

  y2c = GetPosition(wirey2, fMWPCwireStepY2, fMWPC2_Y_offset);

  xtc = x1c - (x2c -x1c)*fMWPCz1/(fMWPCz2-fMWPCz1);

  ytc = y1c + (xtc - x1c)*(y2c-y1c)/(x2c-x1c);

  fXt = xtc;

  fYt = ytc;

}

Float_t GetPosition(Float_t wire, Float_t wireStep,

    Float_t planeOffset) {

  //TODO: number of wires (16) as parameter

  //TODO: omit gRandom

  return (wire - 16.5)*wireStep + planeOffset;

}

void fillSi(TClonesArray *data,Float_t* amp,Float_t* time,Float_t *par1,Float_t *par2,Int_t *multiplicity, Float_t ampThreshold,Float_t timeThreshold=0) {

  if(!data) return;

  if(data->GetEntriesFast()==0) return;

  Bool_t isCal = kTRUE;

  if(par1 == NULL || par2 == NULL) isCal = kFALSE;

  Int_t nCh;

  for(Int_t i=0;i<data->GetEntriesFast();i++) { 

    nCh = ((ERQTelescopeSiDigi*)data->At(i))->Channel();

    // cout << 1000*((ERQTelescopeSiDigi*)data->At(i))->Edep()*(*(par2+nCh)) + (*(par1+nCh)) << " " << nCh << " " << *(par2+nCh) << " " << *(par1+nCh) << endl;

    if (isCal) *(amp+nCh) = ((ERQTelescopeSiDigi*)data->At(i))->Edep();

    if (!isCal) *(amp+nCh) = ((ERQTelescopeSiDigi*)data->At(i))->Edep();

    *(time+nCh) = ((ERQTelescopeSiDigi*)data->At(i))->Time();

    if ( isCal && ( *(amp+nCh)<=ampThreshold || *(time+nCh)<=timeThreshold ) ) {

      *(amp+nCh) = 0;

      *(time+nCh) = 0;

      *multiplicity = *multiplicity - 1;

    } 

  } 

}

void fillarrayCsI(TClonesArray *data,Float_t* amp,Float_t* time,Float_t *par1,Float_t *par2) {

  if(!data) return;

  if(data->GetEntriesFast()==0) return;

  Int_t nCh;

  for(Int_t i=0;i<data->GetEntriesFast();i++) {

    nCh = ((ERQTelescopeCsIDigi*)data->At(i))->Channel();

    *(amp+nCh) = ((ERQTelescopeCsIDigi*)data->At(i))->Edep();

    *(time+nCh) = ((ERQTelescopeCsIDigi*)data->At(i))->Time();

  }

  return;

}

void processCsI_cal() {

  for(Int_t i=0;i<16;i++){

    if (arCsI[i]>aCsI) {

      tCsI = trCsI[i];

      nCsI = i;

      aCsI = arCsI[i];

    }

  }

  Int_t nMax = 0;

  for(Int_t i=0;i<16;i++){

    if (arCsI[i]==aCsI) {

      nMax++;

    }

  }

  if (nMax!=1) {

     cout << " ##ENTRY## " << endl;

    for(Int_t i=0;i<16;i++) {

      if (arCsI[i]>0 && trCsI[i]>0) cout << i << " " << arCsI[i] << " " << trCsI[i] << endl;

    }

    cout << " maximum choice " << endl;

    cout << nCsI << " " << aCsI << " " << tCsI << " " << endl;

    cout << " nMAX= " << nMax << endl;

    nCsI = -1;

    tCsI = 0;

    aCsI = 0;

  } 

  return;

}

void readPar(TString fileName,Float_t *par1,Float_t *par2,Int_t size=16){

  TString line;

  ifstream myfile;

  Int_t count=0;

  TString file = "/home/ivan/work/macro/h7_1904/parameters/" + fileName + ".cal";

  myfile.open(file.Data());

  while (! myfile.eof() ){

    line.ReadLine(myfile);

    // if(count < 0){

    //   count++;

    //   continue;

    // }

    if(line.IsNull()) break;

    sscanf(line.Data(),"%g %g", par1+count,par2+count);

    count++;

  }

  cout << endl << fileName.Data() << endl;

  for(Int_t i=0;i<size;i++) cout << par1[i] << " " << par2[i] << endl;

  return;

}