[u/mrichter/AliRoot.git] / ITS / ITSSPDVertexDiamondda.cxx

/*
Contact: cvetan.cheshkov@cern.ch
Link: missing
Run Type: PHYSICS
DA Type: MON
Number of events needed: 10000
Input Files:
Output Files:
Trigger types used: PHYSICS
*/

#define X_LIMIT 2.0
#define Y_LIMIT 2.0
#define Z_LIMIT 50.0
#define X_DELTA 0.02
#define Y_DELTA 0.02
#define Z_DELTA 0.2
#define OUTPUT_FILE "SPDVertexDiamondDA.root"
#define CDB_STORAGE "local://$ALICE_ROOT"
#define N_EVENTS_AUTOSAVE 50
#define NFITPARAMS 5

extern "C" {
#include "daqDA.h"
}

#include "event.h"
#include "monitor.h"

#ifdef ALI_AMORE
#include <AmoreDA.h>
int amore::da::Updated(char const*) {}
#endif

#include <TTree.h>
#include <TH1.h>
#include <TH2.h>
#include <TF2.h>
#include <TFile.h>
#include <TPluginManager.h>
#include <TROOT.h>
#include <TFitter.h>

#include "AliRawReaderDate.h"
#include "AliGeomManager.h"
#include "AliCDBManager.h"
#include "AliESDVertex.h"
#include "AliITSDetTypeRec.h"
#include "AliITSInitGeometry.h"
#include "AliITSVertexer3DTapan.h"

TF2 *fitFcn = 0x0;

AliESDVertex* FitVertexDiamond(TH2F *hXY, TH1F *hZ);

int main(int argc, char **argv) {

  gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
					"*",
					"TStreamerInfo",
					"RIO",
					"TStreamerInfo()"); 
  TFitter *minuitFit = new TFitter(NFITPARAMS);
  TVirtualFitter::SetFitter(minuitFit);

  int status;
  if (argc<2) {
    printf("Wrong number of arguments\n");
    return -1;
  }

  /* define data source : this is argument 1 */  
  status=monitorSetDataSource( argv[1] );
  if (status!=0) {
    printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
    return -1;
  }

  /* declare monitoring program */
  status=monitorDeclareMp( __FILE__ );
  if (status!=0) {
    printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
    return -1;
  }

  /* define wait event timeout - 1s max */
  monitorSetNowait();
  monitorSetNoWaitNetworkTimeout(1000);
  
  /* log start of process */
  printf("Vertex-Diamond SPD DA started\n");  

  /* init some counters */
  int nevents_with_vertex = 0;
  int nevents_physics=0;
  int nevents_total=0;

  struct eventHeaderStruct *event;
  eventTypeType eventT;

  // Get run number
  if (getenv("DATE_RUN_NUMBER")==0) {
    printf("DATE_RUN_NUMBER not properly set.\n");
    return -1;
  }
  int runNr = atoi(getenv("DATE_RUN_NUMBER"));

  // Histograms initialization
  TH2F *hXY = new TH2F("hXY","Vertex Diamond (Y vs X)",
		       2*(Int_t)(X_LIMIT/X_DELTA),-X_LIMIT,X_LIMIT,
		       2*(Int_t)(Y_LIMIT/Y_DELTA),-Y_LIMIT,Y_LIMIT);
  TH1F *hZ  = new TH1F("hZ"," Longitudinal Vertex Profile",
		       2*(Int_t)(Z_LIMIT/Z_DELTA),-Z_LIMIT,Z_LIMIT);

  // Global initializations
  AliCDBManager *man = AliCDBManager::Instance();
  man->SetDefaultStorage(CDB_STORAGE);
  man->SetRun(runNr);
  AliGeomManager::LoadGeometry("geometry.root");
  AliGeomManager::ApplyAlignObjsFromCDB("ITS");

  // ITS initializations
  AliITSInitGeometry initgeom;
  AliITSgeom *geom = initgeom.CreateAliITSgeom();
  printf("Geometry name: %s\n",(initgeom.GetGeometryName()).Data());

  AliITSDetTypeRec *detTypeRec = new AliITSDetTypeRec();
  detTypeRec->SetITSgeom(geom);
  detTypeRec->SetDefaults();
  detTypeRec->SetDefaultClusterFindersV2(kTRUE);

  // Initialization of AMORE sender
#ifdef ALI_AMORE
  amore::da::AmoreDA vtxAmore(amore::da::AmoreDA::kSender);
#endif
  /* main loop (infinite) */
  for(;;) {
    
    /* check shutdown condition */
    if (daqDA_checkShutdown()) {break;}
    
    /* get next event (blocking call until timeout) */
    status=monitorGetEventDynamic((void **)&event);
    if (status==MON_ERR_EOF) {
      printf ("End of File detected\n");
      break; /* end of monitoring file has been reached */
    }
    
    if (status!=0) {
      printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
      break;
    }
    
    /* retry if got no event */
    if (event==NULL) {
      continue;
    }

    nevents_total++;
    eventT=event->eventType;
    switch (event->eventType){
      
      /* START OF RUN */
    case START_OF_RUN:
      break;
      /* END START OF RUN */
      
    /* END OF RUN */
    case END_OF_RUN:
      break;
      
    case PHYSICS_EVENT:
      nevents_physics++;
      AliRawReader *rawReader = new AliRawReaderDate((void*)event);

      // Run SPD cluster finder
      TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
      detTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");

      // Run vertex-finder
      AliITSVertexer3DTapan *vertexer = new AliITSVertexer3DTapan(geom,1000);
      vertexer->LoadClusters(clustersTree);
      AliESDVertex *vtx = new AliESDVertex();
      vertexer->FindVertexForCurrentEvent(vtx);

      if (TMath::Abs(vtx->GetChi2()) < 0.1) {
	// Fill the vertex into the histos
	nevents_with_vertex++;
	hXY->Fill(vtx->GetXv(),vtx->GetYv());
	hZ->Fill(vtx->GetZv());

	// Auto save
	if ((nevents_with_vertex%N_EVENTS_AUTOSAVE) == 0) {
	  TFile outFile(OUTPUT_FILE, "update");
	  AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ);
	  if (fitVtx) {
	    fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite);
	    TH1 *fithXY = fitFcn->CreateHistogram();
	    fithXY->Write(fithXY->GetName(),TObject::kOverwrite);
	    delete fithXY;
	  }
	  hXY->Write(hXY->GetName(),TObject::kOverwrite);
	  hZ->Write(hZ->GetName(),TObject::kOverwrite);
	  outFile.Close();
	  delete fitVtx;
	}
      }

      delete vtx;
      delete vertexer;

      delete clustersTree;
      delete rawReader;

    }
    
    /* free resources */
    free(event);
    
    /* exit when last event received, no need to wait for TERM signal */
    if (eventT==END_OF_RUN) {
      printf("EOR event detected\n");
      break;
    }
  }
  
  if (detTypeRec) delete detTypeRec;

  // Store the final histograms
  TFile outFile(OUTPUT_FILE, "update");
  if (nevents_with_vertex > N_EVENTS_AUTOSAVE) { 
    // Fit XY & Z histograms
    AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ);
    if (fitVtx) {
      fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite);
      TH1 *fithXY = fitFcn->CreateHistogram();
      fithXY->Write(fithXY->GetName(),TObject::kOverwrite);
      delete fithXY;
    }
    delete fitVtx;
  }
  hXY->Write(hXY->GetName(),TObject::kOverwrite);
  hZ->Write(hZ->GetName(),TObject::kOverwrite);
  outFile.Close();

#ifdef ALI_AMORE
  // send the histos to AMORE pool
  printf("AMORE send status: %d",vtxAmore.Send(hXY->GetName(),hXY));
#endif

  delete minuitFit;
  TVirtualFitter::SetFitter(0);

  /* write report */
  printf("Run #%s, received %d events with vertex, out of %d physics and out of %d total events\n",getenv("DATE_RUN_NUMBER"),nevents_with_vertex,nevents_physics,nevents_total);

  status=0;

  /* export file to FXS */
  if (daqDA_FES_storeFile(OUTPUT_FILE, "VertexDiamond")) {
    status=-2;
  }
  
  return status;
}

Double_t fitFunction(Double_t *x, Double_t *par) {

  Double_t t1 =   x[0] - par[1];
  Double_t t2 =   x[1] - par[2];

  return par[0]*TMath::Exp(-0.5*(t1*t1/(par[3]*par[3])+t2*t2/(par[4]*par[4])));
}

AliESDVertex* FitVertexDiamond(TH2F *hXY, TH1F *hZ)
{

  if (!fitFcn) {
    fitFcn = new TF2("fitFcn",fitFunction,
		     -X_LIMIT,X_LIMIT,
		     -Y_LIMIT,Y_LIMIT,NFITPARAMS);
    fitFcn->SetNpx(2*(Int_t)(X_LIMIT/X_DELTA));
    fitFcn->SetNpy(2*(Int_t)(Y_LIMIT/Y_DELTA));
    fitFcn->SetParameters(hXY->GetMaximum(),0,0,hXY->GetRMS(1),hXY->GetRMS(2));
    fitFcn->Update();
  }

  if (hXY->Fit("fitFcn","L0V+") != 0) {
    printf("XY fit failed!");
    return 0x0;
  }
  
  Double_t pos[3],poserr[3];
  pos[0] = fitFcn->GetParameter(1);
  pos[1] = fitFcn->GetParameter(2);
  poserr[0] = fitFcn->GetParameter(3);
  poserr[1] = fitFcn->GetParameter(4);

  // Could be replaced by something more robust...
  pos[2] = hZ->GetMean();
  poserr[2] = hZ->GetRMS();
 
  return new AliESDVertex(pos,poserr);
}
Commit	Line	Data
eb35e591	1	/*
	2	Contact: cvetan.cheshkov@cern.ch
	3	Link: missing
	4	Run Type: PHYSICS
	5	DA Type: MON
	6	Number of events needed: 10000
	7	Input Files:
	8	Output Files:
	9	Trigger types used: PHYSICS
	10	*/
	11
	12	#define X_LIMIT 2.0
	13	#define Y_LIMIT 2.0
	14	#define Z_LIMIT 50.0
	15	#define X_DELTA 0.02
	16	#define Y_DELTA 0.02
	17	#define Z_DELTA 0.2
	18	#define OUTPUT_FILE "SPDVertexDiamondDA.root"
	19	#define CDB_STORAGE "local://$ALICE_ROOT"
	20	#define N_EVENTS_AUTOSAVE 50
	21	#define NFITPARAMS 5
	22
	23	extern "C" {
	24	#include "daqDA.h"
	25	}
	26
	27	#include "event.h"
	28	#include "monitor.h"
	29
55c5e86d	30	#ifdef ALI_AMORE
	31	#include <AmoreDA.h>
	32	int amore::da::Updated(char const*) {}
	33	#endif
	34
eb35e591	35	#include <TTree.h>
	36	#include <TH1.h>
	37	#include <TH2.h>
	38	#include <TF2.h>
	39	#include <TFile.h>
	40	#include <TPluginManager.h>
	41	#include <TROOT.h>
	42	#include <TFitter.h>
	43
	44	#include "AliRawReaderDate.h"
	45	#include "AliGeomManager.h"
	46	#include "AliCDBManager.h"
	47	#include "AliESDVertex.h"
	48	#include "AliITSDetTypeRec.h"
	49	#include "AliITSInitGeometry.h"
	50	#include "AliITSVertexer3DTapan.h"
	51
	52	TF2 *fitFcn = 0x0;
	53
	54	AliESDVertex* FitVertexDiamond(TH2F hXY, TH1F hZ);
	55
	56	int main(int argc, char **argv) {
	57
	58	gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
	59	"*",
	60	"TStreamerInfo",
	61	"RIO",
	62	"TStreamerInfo()");
	63	TFitter *minuitFit = new TFitter(NFITPARAMS);
	64	TVirtualFitter::SetFitter(minuitFit);
	65
	66	int status;
	67	if (argc<2) {
	68	printf("Wrong number of arguments\n");
	69	return -1;
	70	}
	71
	72	/* define data source : this is argument 1 */
	73	status=monitorSetDataSource( argv[1] );
	74	if (status!=0) {
	75	printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
	76	return -1;
	77	}
	78
	79	/* declare monitoring program */
	80	status=monitorDeclareMp( __FILE__ );
	81	if (status!=0) {
	82	printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
	83	return -1;
	84	}
	85
	86	/* define wait event timeout - 1s max */
	87	monitorSetNowait();
	88	monitorSetNoWaitNetworkTimeout(1000);
	89
	90	/* log start of process */
	91	printf("Vertex-Diamond SPD DA started\n");
	92
	93	/* init some counters */
	94	int nevents_with_vertex = 0;
	95	int nevents_physics=0;
	96	int nevents_total=0;
	97
	98	struct eventHeaderStruct *event;
99	eventTypeType eventT;
100
101	// Get run number
102	if (getenv("DATE_RUN_NUMBER")==0) {
103	printf("DATE_RUN_NUMBER not properly set.\n");
104	return -1;
105	}
106	int runNr = atoi(getenv("DATE_RUN_NUMBER"));
107
108	// Histograms initialization
109	TH2F *hXY = new TH2F("hXY","Vertex Diamond (Y vs X)",
110	2*(Int_t)(X_LIMIT/X_DELTA),-X_LIMIT,X_LIMIT,
111	2*(Int_t)(Y_LIMIT/Y_DELTA),-Y_LIMIT,Y_LIMIT);
112	TH1F *hZ = new TH1F("hZ"," Longitudinal Vertex Profile",
113	2*(Int_t)(Z_LIMIT/Z_DELTA),-Z_LIMIT,Z_LIMIT);
114
115	// Global initializations
116	AliCDBManager *man = AliCDBManager::Instance();
117	man->SetDefaultStorage(CDB_STORAGE);
118	man->SetRun(runNr);
119	AliGeomManager::LoadGeometry("geometry.root");
120	AliGeomManager::ApplyAlignObjsFromCDB("ITS");
121
122	// ITS initializations
123	AliITSInitGeometry initgeom;
124	AliITSgeom *geom = initgeom.CreateAliITSgeom();
125	printf("Geometry name: %s\n",(initgeom.GetGeometryName()).Data());
126
127	AliITSDetTypeRec *detTypeRec = new AliITSDetTypeRec();
128	detTypeRec->SetITSgeom(geom);
129	detTypeRec->SetDefaults();
130	detTypeRec->SetDefaultClusterFindersV2(kTRUE);
131
55c5e86d	132	// Initialization of AMORE sender
	133	#ifdef ALI_AMORE
	134	amore::da::AmoreDA vtxAmore(amore::da::AmoreDA::kSender);
	135	#endif
eb35e591	136	/* main loop (infinite) */
	137	for(;;) {
	138
	139	/* check shutdown condition */
	140	if (daqDA_checkShutdown()) {break;}
	141
	142	/* get next event (blocking call until timeout) */
	143	status=monitorGetEventDynamic((void **)&event);
	144	if (status==MON_ERR_EOF) {
	145	printf ("End of File detected\n");
	146	break; /* end of monitoring file has been reached */
	147	}
	148
	149	if (status!=0) {
	150	printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
	151	break;
	152	}
	153
	154	/* retry if got no event */
	155	if (event==NULL) {
	156	continue;
	157	}
	158
	159	nevents_total++;
	160	eventT=event->eventType;
	161	switch (event->eventType){
	162
	163	/* START OF RUN */
	164	case START_OF_RUN:
	165	break;
	166	/* END START OF RUN */
	167
	168	/* END OF RUN */
	169	case END_OF_RUN:
	170	break;
	171
	172	case PHYSICS_EVENT:
	173	nevents_physics++;
	174	AliRawReader rawReader = new AliRawReaderDate((void)event);
	175
	176	// Run SPD cluster finder
	177	TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
	178	detTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
	179
	180	// Run vertex-finder
	181	AliITSVertexer3DTapan *vertexer = new AliITSVertexer3DTapan(geom,1000);
	182	vertexer->LoadClusters(clustersTree);
	183	AliESDVertex *vtx = new AliESDVertex();
	184	vertexer->FindVertexForCurrentEvent(vtx);
	185
	186	if (TMath::Abs(vtx->GetChi2()) < 0.1) {
	187	// Fill the vertex into the histos
	188	nevents_with_vertex++;
	189	hXY->Fill(vtx->GetXv(),vtx->GetYv());
	190	hZ->Fill(vtx->GetZv());
	191
	192	// Auto save
	193	if ((nevents_with_vertex%N_EVENTS_AUTOSAVE) == 0) {
	194	TFile outFile(OUTPUT_FILE, "update");
	195	AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ);
	196	if (fitVtx) {
	197	fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite);
	198	TH1 *fithXY = fitFcn->CreateHistogram();
	199	fithXY->Write(fithXY->GetName(),TObject::kOverwrite);
200	delete fithXY;
201	}
202	hXY->Write(hXY->GetName(),TObject::kOverwrite);
203	hZ->Write(hZ->GetName(),TObject::kOverwrite);
204	outFile.Close();
205	delete fitVtx;
206	}
207	}
208
209	delete vtx;
210	delete vertexer;
211
212	delete clustersTree;
213	delete rawReader;
214
215	}
216
217	/* free resources */
218	free(event);
219
220	/* exit when last event received, no need to wait for TERM signal */
221	if (eventT==END_OF_RUN) {
222	printf("EOR event detected\n");
223	break;
224	}
225	}
226
227	if (detTypeRec) delete detTypeRec;
228
229	// Store the final histograms
230	TFile outFile(OUTPUT_FILE, "update");
231	if (nevents_with_vertex > N_EVENTS_AUTOSAVE) {
232	// Fit XY & Z histograms
233	AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ);
234	if (fitVtx) {
235	fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite);
236	TH1 *fithXY = fitFcn->CreateHistogram();
237	fithXY->Write(fithXY->GetName(),TObject::kOverwrite);
238	delete fithXY;
239	}
240	delete fitVtx;
241	}
242	hXY->Write(hXY->GetName(),TObject::kOverwrite);
243	hZ->Write(hZ->GetName(),TObject::kOverwrite);
244	outFile.Close();
245
55c5e86d	246	#ifdef ALI_AMORE
	247	// send the histos to AMORE pool
	248	printf("AMORE send status: %d",vtxAmore.Send(hXY->GetName(),hXY));
	249	#endif
	250
eb35e591	251	delete minuitFit;
	252	TVirtualFitter::SetFitter(0);
	253
	254	/* write report */
	255	printf("Run #%s, received %d events with vertex, out of %d physics and out of %d total events\n",getenv("DATE_RUN_NUMBER"),nevents_with_vertex,nevents_physics,nevents_total);
	256
	257	status=0;
	258
	259	/* export file to FXS */
	260	if (daqDA_FES_storeFile(OUTPUT_FILE, "VertexDiamond")) {
	261	status=-2;
	262	}
	263
	264	return status;
	265	}
	266
	267	Double_t fitFunction(Double_t x, Double_t par) {
	268
	269	Double_t t1 = x[0] - par[1];
	270	Double_t t2 = x[1] - par[2];
	271
	272	return par[0]TMath::Exp(-0.5(t1t1/(par[3]par[3])+t2t2/(par[4]par[4])));
	273	}
	274
	275	AliESDVertex* FitVertexDiamond(TH2F hXY, TH1F hZ)
	276	{
	277
	278	if (!fitFcn) {
	279	fitFcn = new TF2("fitFcn",fitFunction,
	280	-X_LIMIT,X_LIMIT,
	281	-Y_LIMIT,Y_LIMIT,NFITPARAMS);
	282	fitFcn->SetNpx(2*(Int_t)(X_LIMIT/X_DELTA));
	283	fitFcn->SetNpy(2*(Int_t)(Y_LIMIT/Y_DELTA));
	284	fitFcn->SetParameters(hXY->GetMaximum(),0,0,hXY->GetRMS(1),hXY->GetRMS(2));
	285	fitFcn->Update();
	286	}
	287
	288	if (hXY->Fit("fitFcn","L0V+") != 0) {
	289	printf("XY fit failed!");
	290	return 0x0;
	291	}
	292
	293	Double_t pos[3],poserr[3];
	294	pos[0] = fitFcn->GetParameter(1);
	295	pos[1] = fitFcn->GetParameter(2);
	296	poserr[0] = fitFcn->GetParameter(3);
	297	poserr[1] = fitFcn->GetParameter(4);
	298
	299	// Could be replaced by something more robust...
	300	pos[2] = hZ->GetMean();
	301	poserr[2] = hZ->GetRMS();
	302
	303	return new AliESDVertex(pos,poserr);
	304	}