ITS/ITSSPDVertexDiamondda.cxx

   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
  30 #ifdef ALI_AMORE
  31 #include <AmoreDA.h>
  32 int amore::da::Updated(char const*) {}
  33 #endif
  34
  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
 132   // Initialization of AMORE sender
 133 #ifdef ALI_AMORE
 134   amore::da::AmoreDA vtxAmore(amore::da::AmoreDA::kSender);
 135 #endif
 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
 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
 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*(t1*t1/(par[3]*par[3])+t2*t2/(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 }